#196 Range reduction for periodic trig functions

[Jeronimo Pellegrini]

Currently,

float(cos(2^10000));
Maxima encountered a Lisp error:
 arithmetic error FLOATING-POINT-OVERFLOW signalled

But it would be nice if Maxima had done range reduction on the number and answered the correct value. For example, pari/gp does this:

? cos(2^10000)
%1 = -0.79208722661996953073200680951110467514

[Jeronimo Pellegrini]

That would likely require using rational arithmetic for large values (so as to be precise enough). Using Taylor expansion, perhaps.

[Stavros Macrakis]

Maxima's default floating-point system is machine floating point (typically 64-bit IEEE 754, ~16 decimal digits of precision, limited to ~2^1024), and uses conventional one-pass bottom-up expression evaluation, so float(cos(2^10000)) is in effect t1 <- 2^10000 (a very large integer); t2 <- float(t1); t3 <- sin(t2). The overflow happens with float(t1), and the cos function never even "sees" the argument.

Maxima also provides an arbitrary fixed-precision floating point type called bfloat or bigfloat. Pari/GP uses a similar system as the default, with 38 digits of precision. To reproduce that in Maxima:

fpprec: 38$
bfloat(cos(2^10000)) =>
      -7.9208722661996953073200680951110467514b-1

and in fact Maxima (and I assume Pari/GP) uses enough precision for π to get the right result even when running at lower precision:

fpprec: 5$
bfloat(cos(2^10000));
  =>  -7.9209b-1

Of course, this still doesn't handle cases like atan(%e^100)-atan(%e^100-1), which requires at least 87 digits of precision to get a non-zero result.

Apparently some systems (Mathematica?) have adaptive precision for cases like that. Maxima doesn't, and I don't believe Pari/GP does either.