Inverting a power series
As a student, one of the things that seemed curious to me about power series was that a function might have a simple power series, but its inverse could be much more complicated. For example, the coefficients in the series for arctangent have a very simple pattern
but the coefficients in the series for tangent are mysterious.
There’s no obvious pattern to the coefficients. It is possible to write the sum in closed form, but this requires introducing the Bernoulli numbers which are only slightly less mysterious than the power series for tangent.
To a calculus student, this is bad news: a simple, familiar function has a complicated power series. But this is good news for combinatorics. Reading the equation from right to left, it says a complicated sequence has a simple generating function!
Computing the coefficients
The example above suggests that inverting a power series, i.e. starting with the power series for a function and computing the power series for its inverse, is going to be complicated. We introduce exponential Bell polynomials to encapsulate some of the complexity, analogous to introducing Bernoulli numbers above.
Assume the function we want to invert, f(x), satisfies f(0) = 0 and its derivative satisfies f‘(0) ≠ 0. Assume f and its inverse g have the series representations
Note that this isn’t the power series per se. The coefficients here are k! times the ordinary power series coefficients. (Compare with ordinary and exponential generating functions explained here.)
Then you can compute the series coefficients for g from the coefficients for f as follows. We have g1 = 1/f1 (things start out easy!) and for n ≥ 2,
where the B‘s are exponential Bell polynomials,
is the kth rising power of n. (More on rising and falling powers here.)
Let’s do an example in Python. Suppose we want a series for the inverse of the gamma function near 2. The equations above assume we’re working in a neighborhood of 0 and that our function is 0 at 0. So we will invert the series for f(x) = Γ(x+2) – 1 and then adjust the result to find the inverse of Γ(x).
import numpy as np from scipy.special import gamma from sympy import factorial, bell def rising_power(a, k): return gamma(a+k)/gamma(a) # Taylor series coefficients for gamma centered at 2 gamma_taylor_coeff = [ 1.00000000000000, 0.42278433509846, 0.41184033042643, 0.08157691924708, 0.07424901075351, -0.0002669820687, 0.01115404571813, -0.0028526458211, 0.00210393334069, -0.0009195738388, 0.00049038845082, ] N = len(gamma_taylor_coeff) f_exp_coeff = np.zeros(N) for i in range(1, N): f_exp_coeff[i] = gamma_taylor_coeff[i]*factorial(i) # Verify the theoretical requirements on f assert( f_exp_coeff == 0 ) assert( f_exp_coeff != 0 ) f_hat = np.zeros(N-1) for k in range(1, N-1): f_hat[k] = f_exp_coeff[k+1] / ((k+1)*f_exp_coeff) g_exp_coeff = np.zeros(N) g_exp_coeff = 1/f_exp_coeff for n in range(2, N): s = 0 for k in range(1, n): # Note that the last argument to bell must be a list, # and not a NumPy array. b = bell(n-1, k, list(f_hat[1:(n-k+1)])) s += (-1)**k * rising_power(n, k) * b g_exp_coeff[n] = s / f_exp_coeff**n def g(x): return sum([ x**i * g_exp_coeff[i] / factorial(i) for i in range(1, N) ]) def gamma_inverse(x): return 2. + g(x-1.) # Verify you end up where you started when you do a round trip for x in [1.9, 2.0, 2.1]: print( gamma_inverse(gamma(x)) )
1.90000003424331 2.00000000000000 2.09999984671755