long double broken on i386?
sgk at troutmask.apl.washington.edu
Mon Dec 3 06:54:35 PST 2007
On Mon, Dec 03, 2007 at 02:44:08AM -0500, David Schultz wrote:
> Is the latest version of the patch the one you sent to the list?
> If cosl and friends produce accurate answers within a reasonable
> domain, it's worth committing; the whole business about how
> accurate cosl(1000000000) is can be dealt with later.
No, I've rewritten most of the code (numerous times in fact :-).
First, I followed Bruce's recommendation to compute the
ULP in the higher precision available in GMP/MPFR instead
of the long double type. The old method appears to
have been using a chop rounding because I would get
only 0.5, 1.0, 1.5, ... ULPs (ie, I never saw anything like
0.6312 ULP). With the new method, I see ULPs in the 0.6
to 0.75 range. I haven't had time to implement a full
blown Remes algorithm to see if I can reduce the ULP.
Second, I've implemented the algorithms found in
s_cos.c, s_sin.c, and s_tan.c that use the extra
precision from the argument reduction. If I use
these algorithms, then the ULPs go up! For example,
the ULP for sinl() goes up to about 1.5. Bruce
had suggested that the 396 hex digits in the table
for 2/pi may need to be checked. I found a paper
by K.C. Ng at David Hough's website that (if I read it
correctly) suggests that we need a much larger
table. So, I need to check how well argument reduction
More information about the freebsd-standards