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 *x*^{2m − 1} then we need

*x*^{2m + 1} / (2*m* + 1)! < ε.

Since we’re interested in *x* = 1 and Γ(*n* + 1) = *n*!, we need

1/ε < Γ(2*m* + 2).

That means we need

2*m* + 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 2*m* + 2 to be larger than Γ^{−1}(10^{100}), and the code above tells us that Γ^{−1}(10^{100}) 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 10^{1000} 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`

= 10^{100} 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).