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

Backward compatibility.

>If I'd know for an algorithm just to work in a specific range (e.g.
>from -1.0e10 to +1.0e10), 

You do understand, that if you want the result to be accurate to
the full 53 bits that it could be accurate, then your approximation
to PI must have about 106 bits of accuracy, even in the range of
0 to 2PI?  Supposing you have the machine approximation to PI, all
53 bits of it, and you want to take the sine of THAT APPROXIMATION,
WHICH IS NOT PI.

sin(x) = -sin(x-PI), so let's subtract PI

Approx pi       =  4 * 0.c90fdaa22168c
Actual pi       =  4 * 0.c90fdaa22168c234c4c6628b80dc1cd129024e088
approx - actual = -4 * 0.0000000000000234c4c6628b80dc1cd129024e088...
-sin(app-actual)=  4 * 0.0000000000000234c4c6628b80dc1cd129024e088...
                =  4 * 2 ** (-52) * 0.234c4c6628b80dc1cd129024e088...
                =      2 ** (-50) * 0.234c4c6628b80dc1cd129024e088...
round to 53     =      2 ** (-50) * 0.234c4c6628b80dc1cd129024e088...
                             +------->  1  2 3  4 5^^ 
        bits in mantissa ----+------->2604826048260|+--+
                             |                     |   |
                             |                     11011100
                             +-------------------->123xy, xy > 1/2
round to 53     =      2 ** (-50) * 0.234c4c6628b80e
                =      2 ** (-53) * 1.1a62633145c07

That's the answer, NOT zero, and it apparently takes 109 bits of PI
to get it.

On the other hand, this means that you must get comfortable with the
fact that except for zero itself, there is no multiple N of machine_PI
such that N*machine_PI has a sine of zero.

If you are working with IEEE "float" inputs, you can work
with a double approximation to PI and generally do ok in
a restricted range near zero, because a 53-bit PI will get you
the correct 24-bit results given a 24-bit input.  This still
isn't good enough for sin of really large inputs.

This problem is not unique to sin/cos/tan, though they're pretty
much the worst case (*).  If you decide to implement complex arithmetic
and just transcribe the book formulas, that'll lose precision as
well.  Consider (a + bi) * (c + di) =  (ac - bd) + (ad + bc)i;
for float-complex you need to put ac, bd, ad, and bc into double
fp to avoid loss of accuracy from cancellation, and you need
"128-bit" FP to accurately compute the result of a double-complex
multiplication.

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