libquad question

[Steve Kargl]

On Mon, Sep 06, 2010 at 06:11:18PM +0200, FX wrote:
> >> The problem is that on FreeBSD (and any other system)
> >> where precision(1.0_10) and precision(1.0_8) are
> >> the same
> 
> On i686-freebsd, you have precision(0._8) == precision(0._10) == 15? How come?
> 
> FX

It goes back to the start of the *BSD.  In src/sys/i386/include/npx.h,
one finds the comment

/*
 * The hardware default control word for i387's and later coprocessors is
 * 0x37F, giving:
 *
 *      round to nearest
 *      64-bit precision
 *      all exceptions masked.
 *
 * We modify the affine mode bit and precision bits in this to give:
 *
 *      affine mode for 287's (if they work at all) (1 in bitfield 1<<12)
 *      53-bit precision (2 in bitfield 3<<8)
 *
 * 64-bit precision often gives bad results with high level languages
 * because it makes the results of calculations depend on whether
 * intermediate values are stored in memory or in FPU registers.
 */

It is my understanding that the issue comes up when one thinks
about the implementation of various libm routines.  Many of
the routines use minimax polynomial approximations on some small
interval.  For example, sin(x) reduces x into the interval [0,pi/4).
The minimax polynomials were developed with the precision of
the type, so allowing the 80-bit precision of the FPU may cause
a problem.  While for the float functions, it is possible to test
all possible inputs for problems and adjust the 'minimax' polynomial,
for double and especially long double it simply is not possible to
test for problems (except over a very limited subset of possible
values).

Code spelunking, finds the 17 year old version of the file

http://www.freebsd.org/cgi/cvsweb.cgi/src/sys/i386/include/npx.h?rev=1.1;content-type=text%2Fplain

Notice the 17 year old comment.

-- 
Steve