David Chase wrote:
> 
> At 11:59 PM 8/4/00 +0200, Stephan Rudlof wrote:
> >Before 'bugfix':
> >        1.000000e+20
> >and
> >        -7.469219e-01;
> >after bugfix by using compile flag -D__NO_MATH_INLINES
> >in both cases
> >        -7.469219e-01
> >.
> >Reading specs from /usr/lib/gcc-lib/i486-linux/egcs-2.91.66/specs
> >gcc version egcs-2.91.66 19990314/Linux (egcs-1.1.2 release)
> 
> >> I got identical behavior, in both cases printing:
> >> sin(arg1): -6.452513e-01

Reading the informations below, this should be the *correct* one! And my
without optimizations is wrong! :-(

> >>
> >> Gcc version:
> >> gcc version egcs-2.91.57 19980901 (egcs-1.1 release)
> >
> >Which machine/CPU?
> 
> Pentium Pro, 200 Mhz, Windows NT Service Pack 6a.
> 
> >> Checking against two Java implementations:
> >> -0.6452512852657808442058
> >> (they both agreed, and if rounded to 7 decimal
> >> digits, they agree with the C program)
> >>
> >> Checking against contructive reals calculator:
> >> -0.6452512852657808442058117113125230074069...
> >> (agrees with Java, if rounded)
> >
> >To illuminate this point:
> >- Does Java use a C lib to implement double arithmetics?
> 
> Java (1.0, 1.1, 1.2) is specified to give results "as if"
> produced by the fdlibm (Freely Distributable Lib Math)
> C library where all the operations (*, /, +, -) had
> strictfp Java semantics.  The specification of java.lang.Math
> in Java 1.3 was "relaxed" (made to fit Sun's implementation,
> sort of), but the old no-wiggle-room specification was
> renamed "StrictMath".  At least one of the two implementations
> I tried conforms exactly to the Java specification on
> strictfp floating point, transcendental functions, and
> floating point formatting.  You have to be very careful
> with C floating point code; the behavior is specified
> badly, and often implemented worse.
> 
> If you want some possible description of what could be
> going wrong, as well as pointers to fdlibm, Hans Boehm's
> Constructive Reals calculator, etc, see
> 
>   http://www.naturalbridge.com/floatingpoint/
> 
> [ This is a plug, this is my employer, and I did most of
>   the work to ensure that all things floating-point work
>   correctly in the product. ]
> 
> I have heard, but never seen an example demonstrating this
> myself, that Java's floating point specification of sin,
> cos, and tan is sometimes a few bits off of the best
> approximation.  I would love to know this for sure, and
> would really like to have access to a library that provided
> the "best" results if fdlibm is not already that library.
> 
> >- What kind of calculator have you used: A calculator emulated by your
> >computer or a separate one (e.g. a Casio)?
> 
> I used Hans Boehm's Constructive Reals Calculator.
> It will calculate as many digits of the exact answer
> as you demand of it.  If you want to know the answer,
> this is the answer, as much of it as you care to scroll,
> wait for, and write down.
> 
> http://oss.sgi.com/projects/crcalc/CRCalc.html

Thank you for this good hint!
So we are able to prove results of libraries...

> You need to view this with a Java-applet-enabled browser;
> I just checked it with Opera 4.02 and JRE 1.2.2, and it
> worked fine for me.

It works for me, too: Netscape® Communicator 4.61.

> In terms of whether sin(1.0e20) is meaningful, if you
> believe that you know the input, and it happens to be
> 1.0e20, then sin(1.0e20) is *a* number (*), not a pile of
> fuzz, and a computer implementation of sin should
> return the nearest approximation of that number.

There is also an ANSI definition of errno set to EDOM, if the argument is
outside the correct domain; this doesn't work for my implementation (see
also other mail of mine), too!

> (*) 1.0e20 is exactly representable as a double, since
> it is equal to 2**20 * 5**20, and log-base-2 of 5**20 is
> less than 47, so it can be represented exactly in the 53
> bits available to an IEEE double.

Thank you for this information; I didn't know this before.
So the lib implementors don't have the excuse, that the precision isn't
sufficient here.

Independently of this point, at least EDOM should be set correctly...

> It's also not just a problem with large numbers; inputs
> near zeros of sin and cos (N * PI for sine,
> PI/2 + N * PI for cosine) lose quite a few bits of
> precision.

:-(


Greetings,

Stephan

> 
> David Chase
> chase at world.std.com
> drchase at alumni.rice.edu
> chase at naturalbridge.com

-- 
Stephan Rudlof (sr at evolgo.de)
   "Genius doesn't work on an assembly line basis.
    You can't simply say, 'Today I will be brilliant.'"
    -- Kirk, "The Ultimate Computer", stardate 4731.3