long double broken on i386?
sgk at troutmask.apl.washington.edu
Fri Sep 28 08:24:51 PDT 2007
So, in my never ending crusade to implement missing
C99 math functions, I decided to tackle sinl, cosl,
and tanl. I have implementations that use argument
reduction from k_rem_pio2.c and Chebyshev interpolation
for the trig functions over the range [0,pi/4].
On amd64, these routines give <= 1/2 ULP for the various
values I've tested in the reduced range. Now, if I test
these on i386, I see between 400 and 1200 ULP. I spent
a week or so trying to understand the descrepancy and "fix"
the problem using a double-double precision arithemetic.
I finally found the problem! /usr/include/float.h on i386
is lying about its numerical representation of long double.
In particular, at least these 3 values appear to be wrong
#define LDBL_MANT_DIG 64
#define LDBL_EPSILON 1.0842021724855044340E-19L
#define LDBL_DIG 18
long double a;
a = 1.L;
for (i = 1; i < 64; i++)
a /= 2;
if (1.L - a == 1.L) break;
printf("%d %Le\n", i, a);
The above results for i386 are consistent with DBL_MANT_DIG
and DBL_EPSILON. Is this intentional, and should float.h be
More information about the freebsd-standards