9

This Mathematica tutorial notes some strange behaviour with calculations such as Sin[10^50]:

N[Sin[10^50]]

> -0.4805

Sin[10.0^50]

> 0.92455

Wait, what?

I know of the infamous $0.1+0.2\not=0.3$ in floating point arithmetic - but I didn't think that would apply for a number like 10.0, surely?

For example, Python doesn't have such a discrepancy:

>>> import math
>>> math.sin(10**50)
-0.4805001434937588
>>> math.sin(10.0**50)
-0.4805001434937588

I then noticed that, if I type in 10.0^50 and hit shift-enter, then copy-paste the result, I see:

9.999999999999997`*^49

Also:

10^50 - 10.0^50
> 4.15384*10^34

$10^{34}$? That's quite a difference. Roughly 41 decillion according to IntegerName[Round[10^50-10.0^50]].

So my main question is: Why does this inaccuracy happen for an integer expressed in machine-precision, rather than it only happening for decimal numbers? Usually issues like $0.2+0.1=0.30000000000000004$ only happen because you can't express $0.1$ in binary - but $10.0$ can be expressed accurately in binary, right? So how does 10.0^50 go so wrong?

numbermaniac
  • 607
  • 7
  • 14
  • 1
    In this case, it's 10.^50 that can't be represented exactly as machine number (although it's indeed an integer). 2^50 or 4^50, however, can. – user202729 Dec 02 '17 at 04:56
  • To add to the remark by @user202729, check the base 2 representation of 10^50:IntegerDigits[10^50, 2] will show that one cannot get the "exact" representation in base 2 using only the 53 bits of a machine double mantissa. The error is, not surprisingly, on the order of machine epsilon, which explains the approx 10^34 difference in that subtraction. – Daniel Lichtblau Dec 02 '17 at 15:51
  • 1
    It is also worth mentioning that the first step in computing the sine of a number far away from zero is to reduce it modulo $2\pi$. Since $\pi$ is an irrational number, that means you need arbitrary-precision arithmetic (and an arbitrarily long approximation to $\pi$) just to do the division and not get a garbage remainder. – zwol Dec 11 '17 at 22:45

1 Answers1

22

The result given by python is completely wrong, as are the results given by Mathematica for machine numbers. To get a correct result, you need to use extended precision numbers:

N[Sin[10^50], 10]
Precision @ %

-0.7896724934

10.

This result is correct to 10 digits. In order to get a correct result, Mathematica needs to use extra precision beyond the 10 that is requested above. One can see this in 2 ways. First, use a high precision argument to Sin:

Sin[10`60^50]
Precision @ %

-0.78967249

8.41064

We start with 60 digits of precision, but the result is only good to 8 digits.

Second, limit the amount of extra precision Mathematica will use:

Block[{$MaxExtraPrecision = 20}, N[Sin[10^50], 10]]

N::meprec: Internal precision limit $MaxExtraPrecision = 20.` reached while evaluating Sin[100000000000000000000000000000000000000000000000000].

0.

With only 20 digits of extra precision, Mathematica is unable to compute any significant digits. With the default (50) Mathematica is able to produce 10 digits as requested.

Block[{$MaxExraPrecision = 50}, N[Sin[10^50], 10]]
Precision @ %

-0.7896724934

10.

Finally, a machine number can only have 16 digits of precision. This is why you get a result on the order of 10^34 when you do the subtraction:

10^50 - 10.^50

4.15384*10^34

Carl Woll
  • 130,679
  • 6
  • 243
  • 355