Strange results about the sine function of a big integer

Question

I want to calculate the sine of a big integer, e.g., Sin[10^50], but found that N[Sin[10^50]], N[Sin[10^50],10], and N[Sin[10^50],20] give different results:

N[Sin[10^50]]
-0.4805

N[Sin[10^50], 10]
-0.7896724934

N[Sin[10^50], 20]
-0.78967249342931008271

It seems that N[Sin[10^50]] did not give a correct result as N[Sin[10^50], 20] uses higher precision and is more likely to give a correct one. But I was curious why N[Sin[10^50], 10] can give a result close to N[Sin[10^50], 20] with a lower precision than the machine precision while N[Sin[10^50]] cannot?

I did some further research by calculating the remainder of 10^50/(2 Pi), in the following way:

remainder := 10^50 - 2 Pi (10^50/(2 Pi) // Floor);

N[remainder]
2.07692*10^34

N[remainder, 10]
4.051867659

N[remainder, 20]
4.0518676593164822448

It is obvious that N[remainder] has a very large error, while N[remainder,10] and N[remainder,20] seem to give the correct result.

I also tried the built-in Mod function, and the results are similar.

So I was curious why the function N with the default machine precision can have such a large error while N with a specified precision even lower than the machine precision can give a good result?

Can anyone tell the reason? Thanks.

Answer 1

To get a feeling for the inaccuracies in IEEE 754 double-precision floating-point numbers, have a look at how the number $10^{50}$ is represented in this form:

x = N[10^50]
(*    1.*10^50    *)

SetPrecision[x, ∞]
(*    100000000000000007629769841091887003294964970946560    *)

So the number $10^{50}$ has been modified by the finite internal representation of machine numbers to be about $10^{50}+7.62977\times10^{33}$. There simply aren't more bits available to do a better job.

Naturally, calculating the sine of this number has nothing to do with calculating the sine of $10^{50}$.

This is approximately what is going on when doing this kind of machine-precision calculations.

Let's do a bit of stochastic math (with a physics spin) on top of this. Converting an exact number $x$ to a machine-precision floating point number $\tilde{x}$ always brings a relative error of about $\pm2^{-54}$. In the absence of more details, let's assume that this relative error has a Gaussian distribution of mean $x$ and standard deviation $\sigma=x\delta$ (with $\delta\approx 2^{-54}$):

$$ p(\tilde{x})=\frac{1}{x\delta\sqrt{2\pi}}e^{-\frac{(\tilde{x}/x-1)^2}{2\delta^2}} $$

Calculating $\sin(x)$ can then proceed over this distribution:

$$ \langle\sin(\tilde{x})\rangle=\int_{-\infty}^{\infty}\sin(\tilde{x})p(\tilde{x})d\tilde{x} = \sin(x)\cdot e^{-\frac12x^2\delta^2} $$

This means that the expected result $\sin(x)$ is attenuated (smeared out) by the uncertainty factor $e^{-\frac12x^2\delta^2}$. In your case, with $x=10^{50}$ and $\delta\approx2^{-54}$, we get

E^(-1/2 x^2 δ^2) /. {x -> 10^50, δ -> 2^-54}
(*    1/E^(7888609052210118054117285652827862296732064351090230047702789306640625/512)    *)

Log[10, %] // N
(*    -6.69137*10^66    *)

This is a very small number, about $10^{-6.69\times10^{66}}$. So there is absolutely no chance that you'll be getting the right answer $\sin(x)$ from this calculation.

Answer 2

To add to @Roman's answer.

1. Sin[x] is accurate when x is a machine precision number. When x comes from rounding an exact number x0 to the nearest machine-precision number x, the error Sin[x] - Sin[x0] may be as large as x * $MachineEpsilon/2, which is approximately 10^-16 * x. (For very large x, the error is at most to 2 due to the range of Sin[].) Once the roundoff error gets large enough, it's not really worth evaluating Sin[x], it seems to me.

2. N[Sin[10^50]] first converts the argument of Sin[] to machine precision, which TracePrint[ N[Sin[10^50]], _Sin, TraceInternal -> True ] shows. Sin[1.*^50] evaluates to -0.4805. As @Roman pointed out, the roundoff error is well over 10^33, and that is why the value is far from the exact Sin[10^50].

3. N[expr, p] uses adaptive precision to achieve a numerical value that is accurate to a precision of p, when possible. This may require internal computations at a much higher precision than p. So the output N[Sin[10^50], 5] is accurate to 5 digits, but internally, it uses much higher precision. According to the error message from Block[{$MaxExtraPrecision = 49}, N[Sin[10^50], 5]], it needs more than 49 extra digits to obtain five-digit precision. For very large integers x, N[Sin[x], p], for whatever precision p you desire, is the best way to go in Mathematica. And N@N[Sin[x], $MachinePrecision] will give a machine real that is accurate to machine precision.