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