[Larceny-users] A Floating Hole

[David Rush]

I am currently moderately annoyed that:

[David Rush]

I am currently moderately annoyed that

(cos (acos 0))

6.123233995736766e-17

on Larceny. Now it appears that this is an inescapable feature of the
C runtime; however, given that FLT_MIN on the x86 platform is 1E-37,
mapping this result to #i0 is somewhat problematic. Since I'm finding
that neither R5RS nor R6RS are really offering much help on this score
(the requirement appear to be that the implementation not lose
precision from that present in the arguments) I'm wondering if there
is a Larceny-specific way to either

1) get a more accurate cos function (e.g. POSIX cosl which returns a
long double)
2) extract some kind of underlying precision information about the
Larceny's flonums (and maybe individual operations)

I suspect the answers to both of the above are probably "no" (modulo
an excursion into the FFI). But I am trying to use Larceny for
numerics these days and I would be interested as to what the state of
the art best practices are...

david rush

[William D Clinger]

Post by David Rush
I am currently moderately annoyed that

(cos (acos 0))

6.123233995736766e-17
on Larceny. Now it appears that this is an inescapable feature of the
C runtime;

Indeed, it illustrates the inescapable misfeature known as
roundoff error. Most real numbers, including pi/2, are not
representable in binary floating point of any finite precision.
The consequent small errors are more likely to accumulate
during subsequent computations than to cancel themselves out.

Post by David Rush
however, given that FLT_MIN on the x86 platform is 1E-37,
mapping this result to #i0 is somewhat problematic.

FLT_MIN isn't relevant here, since Larceny uses IEEE double
precision. Were Larceny to use single precision, as you are
assuming, the roundoff errors would be much worse.

Post by David Rush
I'm wondering if there
is a Larceny-specific way to either
1) get a more accurate cos function (e.g. POSIX cosl which returns a
long double)

No, although you could use Larceny's FFI to do that.

Post by David Rush
2) extract some kind of underlying precision information about the
Larceny's flonums (and maybe individual operations)

(system-features) will give you that information. Currently,
however, all supported varieties of Larceny use IEEE double
precision on all supported architectures.

Perhaps you will feel better about this once you realize that
C programs experience the same kinds of roundoff error. For
example:

#include <stdlib.h>
#include <math.h>
#include <stdio.h>

int main( int argc, char **argv ) {
printf("%g\n", cos(acos(0.0)));
}

On the two machines I tried, the result was 6.12323e-17,
which is consistent with Larceny's result. Note that C
uses double or extended precision for the above program.

Other data points for (cos (acos 0)):

Bigloo 6.1230317691119e-17
Chicken 6.12303176911189e-17
Gambit-C 6.123031769111886e-17
Ikarus 6.123031769111886e-17
Kawa 6.123233995736766E-17
Larceny 6.123031769111886e-17
MIT 6.123031769111886e-17
PLT 6.123031769111886e-17
Petite 6.123031769111886e-17
Scheme48 6.123031769111886e-17
scm 6.123031769111886e-17
Ypsilon 6.123031769111886e-17

The outputs shown for Bigloo and Chicken illustrate
their violation of R5RS 6.2.6. Kawa's result is the
same as what I get from C, which may mean Kawa uses
an extended precision for the intermediate result
(of acos) instead of rounding to double precision;
that usually improves accuracy or has no effect, but
occasionally gives a less accurate result.

Will

[David Rush]

Post by William D Clinger

Post by David Rush

(cos (acos 0))

6.123233995736766e-17
on Larceny. Now it appears that this is an inescapable feature of the
C runtime;

Indeed, it illustrates the inescapable misfeature known as
roundoff error.

Yup. I actually got here by using an exact rational approximation to
pi, derived via the Machin formula and accurate to significantly more
places than those returned by acos, which produces the same result. I
just didn't really think that 266 digits of ratnum was going to make
the post any clearer.

Post by William D Clinger

Post by David Rush
however, given that FLT_MIN on the x86 platform is 1E-37,
mapping this result to #i0 is somewhat problematic.

FLT_MIN isn't relevant here

Well, the man page on my Debian Etch says that FLT_MIN == DBL_MIN ==
LDBL_MIN == 1E-37, which seems silly but there you go. In any case,
the cos approximation is rather less accurate than any of those
minima.

Post by William D Clinger

Post by David Rush
2) extract some kind of underlying precision information about the
Larceny's flonums (and maybe individual operations)

(system-features) will give you that information.  Currently,
however, all supported varieties of Larceny use IEEE double
precision on all supported architectures.

That's a nice feature. It would be nice to have it include a
flonum-epsilon entry or something similar so that programs could avoid
building a mechanism to convert the symbolic information into
computable data.

Post by William D Clinger
Perhaps you will feel better about this once you realize that
C programs experience the same kinds of roundoff error.

Naah, I figured that part out already. What I'm really looking for is
a way to minimize the roundoff error in trig functions (so I can
actually find the zeros). It's looking like reducing all trig
evaluations to the +x+y quadrant and using symmetry to return to the
actual value might be the best answer. _Numerical Recipes_ also
recommends using a mutant recurrence relation for iterating sines and
cosines over equally spaced angular samples, and I'm thinking of
trying that.

david rush

[William D Clinger]

Post by David Rush
Well, the man page on my Debian Etch says that FLT_MIN == DBL_MIN ==
LDBL_MIN == 1E-37, which seems silly but there you go. In any case,
the cos approximation is rather less accurate than any of those
minima.

You need to understand that Larceny's value for (acos 0)
is the best IEEE double precision approximation to pi/2.
It is not pi/2, however, so the cos procedure correctly
returns a nonzero value for (acos 0). If cos were to
return 0.0 for (acos 0), then cos would be *less* accurate
than is currently the case.

Post by David Rush
Naah, I figured that part out already. What I'm really looking for is
a way to minimize the roundoff error in trig functions (so I can
actually find the zeros). It's looking like reducing all trig
evaluations to the +x+y quadrant and using symmetry to return to the
actual value might be the best answer. _Numerical Recipes_ also
recommends using a mutant recurrence relation for iterating sines and
cosines over equally spaced angular samples, and I'm thinking of
trying that.

As explained above, the only way to do better is to use
more precision throughout, or to use something other than
binary floating point.

Larceny will continue to represent inexact reals in IEEE
floating point. Given that design decision, which is
regarded as best current practice, this isn't something
Larceny can fix.

Will