David,

David Chase wrote:
> 
> At 05:33 PM 8/6/00 +0200, Stephan Rudlof wrote:
> >> >> 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! :-(
> 
> > > -7.469219e-01;
> 
> >> (*) 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.
> 
> It is actually a problem of precision,

Thank you for the clearification.

> even when the input can
> be exactly represented as a machine floating point number.  In order
> to get the correct approximation for a number as large as 1e20
> (about 2**66.4) you need to take the remainder (modulo 2*PI --
> the algorithm in fdlibm is actually mod PI/4 plus quadrant
> information, I think) and to get the remainder accurately you
> need roughly 53+66.4 = ~129 bits (more or less, I am not an
> expert here) of accuracy in your "PI".  Pentiums use 68 bits,
> so you've got *roughly* 1.5 bits of accuracy in your
> input to sin.

I have seen such an algorithm in the ST implementation for failing the sin
primitive, which didn't work (because it calls double routines again), too.

> For me, there are three inputs that should trigger EDOM:
> +Inf, -Inf, and NaN.  People knew how to do range reduction
> properly at least ten years ago, and fdlibm is about
> that old (fdlibm came from the Sun floating point group, and
> I worked at Sun from 1990-1993, and talked to them often).
> 
> As to whether a programmer has any business asking for the
> sine of 1e20, I have no idea, that is not my department.
> It seems to me that returning a wrong answer (not the right
> answer, not NaN, not throwing an exception) is a very bad idea
> indeed.

Throwing an exception by signals could result into a SIGFPE for ANSI-C, so
it seems to be reasonable to give a SIGFPE, an EDOM or NaN here.

If I'd know for an algorithm just to work in a specific range (e.g. from
-1.0e10 to +1.0e10), I'd like to get an EDOM out of this range: It'd be a
practical, not a mathematical EDOM, of course (from my K&R ANSI-C book it
isn't clear, if EDOM has a pure mathematical or a numerical definition). NaN
would have the advantage of giving *no* answer - you'd see that something
goes wrong here - in the common case of not checking errno though. But
probably an EDOM doesn't forbid to give also a NaN as result! Drawback: I
haven't seen any NaN in my ANSI-C book... 

I'm somewhat surprised, that these *basic* problems exist after decades of
numerical computations...

Stephan

> David Chase                    --  chase at naturalbridge.com
> NaturalBridge Inc              --  http://www.naturalbridge.com
> BulletTrain bytecode compiler  --  when you can't wait for performance

-- 
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