long double broken on i386?
Steve Kargl
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
#include <stdio.h>
int
main (void)
{
int i;
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);
return 0;
}
amd64:sgk[206] ./z
64 1.084202e-19
i386:kargl[205] ./z
54 5.551115e-17
The above results for i386 are consistent with DBL_MANT_DIG
and DBL_EPSILON. Is this intentional, and should float.h be
fixed?
--
Steve
More information about the freebsd-standards
mailing list