How do you evaluate the sine of a large number in floating point arithmetic? What does the result even mean?
Sine of a trillion
Let’s start by finding the sine of a trillion (1012) using floating point arithmetic. There are a couple ways to think about this. The floating point number t = 1.0e12 can only be represented to 15 or 16 significant decimal figures (to be precise, 53 bits [1]), and so you could think of it as a representative of the interval of numbers that all have the same floating point representation. Any result that is the sine of a number in that interval should be considered correct.
Another way to look at this is to say that t can be represented exactly—its binary representation requires 42 bits and we have 53 bits of significand precision available—and so we expect sin(t) to return the sine of exactly one trillion, correct to full precision.
It turns out that IEEE arithmetic does the latter, computing sin(1e12) correctly to 16 digits.
Here’s the result in Python
>>> sin(1.0e12)
-0.6112387023768895
and verified by Mathematica by computing the value to 20 decimal places
In:= N[Sin[10^12], 20]
Out= -0.61123870237688949819
Range reduction
Note that the result above is not what you’d get if you were first to take the remainder when dividing by 2π and then taking the sine.
>>> sin(1.0e12 % (2*pi))
-0.6112078499756778
This makes sense. The result of dividing a trillion by the floating point representation of 2π, 159154943091.89536, is correct to full floating point precision. But most of that precision is in the integer part. The fractional part is only has five digits of precision, and so we should expect the result above to be correct to at most five digits. In fact it’s accurate to four digits.
When your processor computes sin(1e12) it does not naively take the remainder by 2π and then compute the sine. If it did, you’d get four significant figures rather than 16.
We started out by saying that there were two ways of looking at the problem, and according to the first one, returning only four significant figures would be quite reasonable. If we think of the value t as a measured quantity, measured to the precision of floating point arithmetic (though hardly anything can be measured that accurately), then four significant figures would be all we should expect. But the people who designed the sine function implementation on your processor did more than they might be expected to do, finding the sine of a trillion to full precision. They do this using a range reduction algorithm that retains far more precision than naively doing a division. (I worked on this problem a long time ago.)
Sine of a googol?
What if we try to take the sine of a ridiculously large number like a googol (10100)? This won’t work because a googol cannot be represented exactly as a floating point number (i.e. as a IEEE 754 double). It’s not too big; floating point numbers can be as large as around 10308. The problem is that a number that big cannot be represented to full precision. But we can represent 2333 exactly, and a googol is between 2332 and 2333. And amazingly, Python’s sine function (or rather the sine function built into the AMD processor on my computer) returns the correct result to full precision.
>>> sin(2**333)
0.9731320373846353
How do we know this is right? I verified it in Mathematica:
In:= Sin[2.0^333]
Out= 0.9731320373846353
How do we know Mathematica is right? We can do naive range reduction using extended precision, say 200 decimal places.
In:= p = N[Pi, 200]
In:= theta = x - IntegerPart[x/ (2 p)] 2 p
Out= 1.8031286178334699559384196689…
and verify that the sine of the range reduced value is correct.
In:= Sin[theta]
Out= 0.97313203738463526…
Interval arithmetic interpretation
Because floating point numbers have 53 bits of precision, all real numbers between 256 − 22 and 256 + 22 are represented as the floating point number 256. This is a range of width 8, wider than 2π, and so the sines of the numbers in this range cover the possible values of sine, i.e. [−1, 1]. So if you think of a floating point number as a sample from the set of all real numbers with the same binary representation, every possible value of sine is a valid return value for numbers bigger than 256. But if you think of a floating point number as representing itself exactly, then it makes sense to ask for the sine of numbers like 2333 and larger, up to the limits of floating point representation [2].
More numerical math posts
[1] An IEEE 754 double has 52 significand bits, but these bits can represent 53 bits of precision because the first bit of the fractional part is always 1, and so it does not need to be represented. See more here.
[2] The largest exponent of an IEEE double is 1023, and the largest significand is 2 − 2−53 (i.e. all 1’s), so the largest possible value of a double is
(253 − 1)21024−53
and in fact the Python expression sin((2**53 - 1)*2**(1024-53)) returns the correct value to 16 significant figures.
