Suppose you wanted to calculate sin(x) to a million decimal places using a Taylor series. How many terms of the series would you need?
You can use trig identities to reduce the problem to finding sin(x) for |x| ≤ 1. Let’s take the worst case and assume we want to calculate sin(1).
The series for sine alternates, and so by the alternating series theorem we need the first term left out of our Taylor approximation to be less than our error tolerate ε. If we keep the terms of the Taylor series up to x2m – 1 then we need
x2m + 1 / (2m + 1)! < ε.
Since we’re interested in x = 1 and Γ(n + 1) = n!, we need
1/ε < Γ(2m + 2).
That means we need
2m + 2 > Γ-1(1/ε).
But how do you compute the inverse of the gamma function? This is something I wrote about a few years ago.
Here is a function for approximately computing the inverse of the gamma function. See the earlier post for details.
from numpy import pi, e, sqrt, log
from scipy.special import lambertw
def inverse_gamma(x):
c = 0.03653381448490056
L = log((x+c)/sqrt(2*pi))
return L / (lambertw(L/e)) + 0.5
Suppose we want to compute sin(1) to 100 decimal places. We need 2m + 2 to be larger than Γ-1(10100), and the code above tells us that Γ-1(10100) is something like 70.9. This tells us we can choose m = 35.
If we want to compute sin(1) to thousands of digits, the code above will fail because we cannot represent 101000 as a floating point number. I will assume for this post that we will use an extended precision library for summing the series for sin(1), but we’re going to use ordinary precision to plan this calculation, i.e. to decide how many terms to sum.
If we look closely at the function inverse_gamma above we see that it only depends on x via log(x + c). Since we’re interested in large x, we can ignore the difference between log(x + c) and log(x). This lets us write a new version of inverse_gamma that takes the log of x rather than x as an argument.
def inverse_gamma2(logx):
L = logx - log(sqrt(2*pi))
return L/lambertw(L/e) + 0.5
Calling inverse_gamma with x = 10100 gives the same result, down to the last decimal place, as calling inverse_gamma2 with 100 log(10).
We asked at the top of the post about computing sine to a million decimal places. If we call inverse_gamma2(1e6*log(10)) we get 205023.17… and so m = 102,511 would be large enough.
***
If you enjoyed reading this post, you may like reading this post on planning a world record calculation for computing ζ(3).