Suppose you need to write software to compute the sine function. You were told in a calculus class that this is done using Taylor series—it’s not, but that’s another story—and so you start writing code to implement Taylor series.
How many terms do you need? Uh, …, 20? Let’s go with that.
from math import sin, factorial
def sine(x):
s = 0
N = 20
for n in range(N//2):
m = 2*n + 1
s += (-1)**n * x**m / factorial(m)
return s
You test your code against the built-in sin function and things look good.
>>> sin(1)
0.8414709848078965
>>> sine(1)
0.8414709848078965
Larger arguments
Next, just for grins, you try x = 100.
>>> sin(100)
-0.5063656411097588
>>> sine(100)
-7.94697857233433e+20
Hmm. Off by 20 orders of magnitude.
So you do a little pencil and paper calculation and determine that the code should work if you set N = 300 in the code above. You try again, and now you’re off by 25 orders of magnitude. And when you try entering 100.0 rather than 100, the program crashes.
What’s going on? In theory, setting N = 300 should work, according to the remainder term in the Taylor series for sine. But the largest representable floating point number is on the order of 10308 and so x**m terms in the code overflow. The reason the program didn’t crash when you typed in 100 with no decimal point is that Python executed x**m in integer arithmetic.
If you were patient and willing to use a huge number of terms in your Taylor series, you couldn’t include enough terms in your series to get accurate results for even moderately large numbers because overflow would kill you before you could get the series truncation error down. Unless you used extended precision floating point. But this would be using heroic measures to rescue a doomed approach.
Taylor series are good local approximations. They’re only valid within a radius of convergence; outside of that radius they’re worthless even in theory. Inside the radius they may be valid in theory but not in practice, as is the case above.
Range reduction
Taylor series isn’t the best approach for evaluating sine numerically, but you could make it work if you limited yourself to a small range of angles, say 0 to π/4, and used trig identities to reduce the general problem to computing sine in the range [0, π/4]. So how do you do that?
Let’s look at reducing arguments to live in the interval [0, 2π]. In practice you’d want to reduce the range further, but that requires keeping up with what quadrant you’re in etc. We’ll just consider [0, 2π] for simplicity.
The natural thing to do is subtract off the largest multiple of 2π you can.
>>> from math import pi
>>> k = int(100/(2*pi))
>>> x = 100 - k*2*pi
>>> sine(x)
-0.5065317636696883
>>> sin(100)
-0.5063656411097588
That’s better. Now we don’t overflow, and we get three correct decimal places. We could increase N and get more. However, there’s a limit to how well we could do, as the calculation below shows.
>>> sin(x)
-0.5063656411097495
>>> sin(100)
-0.5063656411097588
When we reduce 100 mod 2π and call sin(x), we get a result that differs from sin(100) in the last three decimal places. That’s not terrible, depending on the application, but the loss of precision increases with the size of x and would be a complete loss for large enough values of x.
Let’s look back at where we divided 100 by 2π:
>>> 100/(2*pi)
15.915494309189533
When we throw away the integer part (15) and keep the rest (.915494309189533) we lose two figures of precision. The larger the integer part, the more precision we lose.
Our naive attempt at range reduction doesn’t work so well, but there are clever ways to make range reduction work. It is possible, for example, to compute the sine of a googol (i.e. sin 10100).
At first it seems simple to reduce the range. Then when you become aware of the precision problems and think about it a while, it seems precise range reduction is impossible. But if you keep thinking about it long enough you might find a way to make it work. A lot of effort went into range reduction algorithms years ago, and now they’re part of the computational infrastructure that we take for granted.
![D\bigl(f^{-1}\bigr)(y) = \left[Df\bigl(f^{-1}(y)\bigr)\right]^{-1}](https://www.johndcook.com/ift3.svg){kind=small}
