Error in Sin(x)

There seems to be a problem with the small precision calculations:

PariGP 2.4.3 gives sin(10^22) = -0.8522008497671888017727058938
with default precision of 28 decimal digits.

TTCalc 0.9.1 gives sin(10^22) = -0.852200868269788305456900652
with small precision of 26 decimal digits.

This indicates this result has only about 7 digits accuracy.

The accuracy of the PariGP result is confirmed by TTCalc at higher precisions. Indeed, Matlab 2007 gets a correctly-rounded 16-digit result using just IEEE double precision.

Derek O'Connor

Added by: tomek, 2010 III 23, Last modified: 2010 III 23

Unfortunately at the moment there is a very poor implementation of the Mod() function for floating point numbers. The Mod() is calculated in this way:
Mod(x; y) = x - int(x/y) * y

This formula gives inaccurate results when the first operand x is much larger than y. If you calculate sin(10^22) then first the Sin has to reduce 2*PI period. It is done by using Mod(x;2*PI) and the result from Mod is then pass back to the Sin. And here is the problem.

Try to calculate:
mod(10^22; 2*pi)
ttcalc with medium precision gives:
5.2630079146204995036070847812778413131302932913465264190190 15937[...]
and then change precision to small and try to calculate Sin with this argument:
sin(5.2630......) = -0.852200849767188801772705894
which is a good aproximation of sin(10^22)

The Mod function will be the first thing to change in the next release of ttmath.

Added by: ~Derek O'Connor, 2010 III 24

The real error here is that your argument(x) reduction is inaccurate. This means that all trigonometric functions will be inaccurate.

The Ng's paper below shows how difficult it is to do this correctly:

Added by: tomek, 2010 III 30, Last modified: 2010 III 30

Thanks for the link, your site has a lot interesting materials.
In this paper there is shown algorithm for double precision (64 bit floating point numbers). Im not sure if this is applicable to ttmath as we can use arbitrary large numbers. What about sin(10^2200) or even sin(10^20000000000)?

But now I know that the problem is called 'argument reduction' and I make some research in this subject.