Welcome back! I am experimenting with out writing these posts in Jupyter and then exporting to HTML, in hopes that this will be more convenient in the future.
In this post, I want to cover a mind blowing way (at least to me) of doing combinatorics: generating functions. Generating functions are a established combinatorial tool, and frankly the only way an analyst such as myself can stomach combinatorics. For those of you who are unfamiliar with generating functions, you can probably infer a lot about them from the contents of this post, and if you hunger to more of learn these dark arts you can read the legendary generatingfunctionology book (link redirects to the author's official webpage).
While I have a good working knowledge of generating functions, I recently learned another one of their utilities while completing my preliminary examination at UC Berkeley. The problem which inspired this blog post is fairly simple, and anyone with some generating functions experience should recognize how to solve this problem.
Let $c_n$ be the number of ways to make $n$ cents with pennies, nickels, dimes, and quarters. Compute $$\lim_{n\to \infty} \frac{c_n}{n^3}.$$
My solution, unsurprisingly, uses generating functions:
Let $$f(z) = \sum_{n\geq 0} c_n z^n.$$ Then $$f(z) = (1 + z + z^2 + \cdots)(1 + z^5 + z^{10} + \cdots)(1 + z^{10} + z^{20} + \cdots)(1 + z^{25} + z^{50} + \cdots)= \frac{1}{(1-z)(1-z^5)(1-z^{10})(1-z^{25})}.$$ So $f$ is a meromorphic function with a pole of order $4$ at $z = 1$, and all the other poles have order at most $3$. So $$f(z) = \frac{c}{(1 - z)^4} + \sum_p \frac{c_p}{(p - z)^{k_p}} + g(z),$$ where the sum is over all the poles of $f$ of order at most $3$, with $1 \leq k_p \leq 3$ being the multiplicities of each pole. The function $g$ is an entire function, and $c, c_p$ are constants. Taking derivatives yields $$c_n = \frac{f^{(n)}(0)}{n!} = c\binom{n+3}{3} + \sum_p c_p\binom{n + k_p -1}{k_p - 1} + \frac{g^{(n)}(0)}{n!}.$$ Thus $$\frac{c_n}{n^3} \sim \frac{c}{n^3}\binom{n+3}{3}$$ since $1 \leq k_p \leq 3$ implies $\binom{n + k_p -1}{k_p - 1}$ grows at most quadratically in $n$ and since $g$ is entire $\frac{g^{(n)}(0)}{n^3n!} \to 0$. So we get $$\lim_{n\to \infty} \frac{c_n}{n^3} = \frac{c}{3!}.$$ To compute $c$, we observe that $c = \lim_{z \to 1}(1 - z)^4f(z)$, and as $$(1 - z)^4f(z) = \frac{1}{(z^4 + \cdots + z + 1)(z^9 + \cdots + z + 1)(z^{24} + \cdots + z + 1)}$$ we get $c = \frac{1}{5 \cdot 10 \cdot 25}$. Thus $$\lim_{n\to \infty} \frac{c_n}{n^3} = \frac{1}{3!\cdot 5 \cdot 10 \cdot 25} = \frac{1}{7500}.$$
This technique is nothing special, it is the standard way to do such problems using generating functions. That being said, I wanted to numerically verify this answer. In the interest of time, I tried the rather stupid but easy to implement brute force method for calculating the $c_n$, which takes ages to converge:
def c(k):
total = 0
for p in range(k+1):
for n in range((k)//5 +1):
for d in range((k)//10 + 1):
for q in range((k )//25 + 1):
if p + 5 * n + 10 * d + 25 * q == k:
total += 1
return total
result = %timeit -o c(100)
print(result.best) # best time in seconds
16.6 ms ± 108 µs per loop (mean ± std. dev. of 7 runs, 100 loops each) 0.01653641791
Since the complexity of $c(k)$ scales quartically in $k$, this method is too inefficient to be of any practical use for large $k$. Alas, the generating function presents us a way of computing $c_n$ in a much more efficient manner. We know from complex analysis that if $f(z) = \sum_{n\geq 0} c_n z^n,$ then $$c_n = \frac{1}{2\pi i}\int_{|z|=r}\frac{f(z)}{z^{n+1}}\,dz,$$ where $0 < r < 1$ since $f$ has poles on $|z| = 1$. Parameterizing this integral gives $$c_n = \frac{r^{-n}}{2\pi}\int_0^{2\pi} f(re^{i\theta})e^{-in\theta}\,d\theta.$$ The astute reader should recognize this as a Fourier transform. In particular, $c_n$ is just the $n$-th Fourier coefficient of the function $\theta \mapsto r^{-n}f(re^{i\theta})$. Using the approximation $$\frac{r^{-n}}{2\pi}\int_0^{2\pi} f(re^{i\theta})e^{-in\theta}\,d\theta \approx \frac{r^{-n}}{2\pi N}\sum_{k=0}^N f(re^{2\pi ij/N})e^{-2\pi inj/N}$$ for large enough $N$, we can approximate $c_n$ using the fast fourier transform:
import numpy as np
f = lambda z : 1./((1 - z) * (1 - z**5) * (1 - z**10) * (1 - z**25))
def fft_c(f, k, N=4096, R=0.99):
'''
Computes the k-th coefficient of the power series expansion for an analytic function f, provided said power
series has radius of convergence > R.
f: the function
k: the k-th coefficient
N: number of partitions for the numerical integration
R: a positive real number smaller than the radius of convergence.
'''
j = np.arange(N)
theta = 2.0 * np.pi * j/N
z = R * np.exp(1j * theta)
S = np.fft.fft(f(z))
# Extract the k-th DFT coefficient (mod N)
sum_term = S[k % N]
a_k = (1.0 / N) * (R ** (-k)) * sum_term
return a_k.real
# sanity check on small k
Nval = lambda k : max(4096, 2**(int(np.log2(2*k+1)) + 1)) # A function for computing a large N:
k = 25
N = Nval(k)
# we round the values of the FFT method since they will be some stuff after the decimal place
print(f'Actual values of cn: {[c(j) for j in range(1,k+1)]}\nFFT computed values: {[round(fft_c(f, j, N=N)) for j in range(1,k+1)]}')
Actual values of cn: [1, 1, 1, 1, 2, 2, 2, 2, 2, 4, 4, 4, 4, 4, 6, 6, 6, 6, 6, 9, 9, 9, 9, 9, 13] FFT computed values: [1, 1, 1, 1, 2, 2, 2, 2, 2, 4, 4, 4, 4, 4, 6, 6, 6, 6, 6, 9, 9, 9, 9, 9, 13]
## Just to prove I am not hiding anything by rounding the FFT method, here are the unrounded values:
print(f'Unrounded values: {[fft_c(f, j, N=N) for j in range(1,k+1)]}')
Unrounded values: [1.0000000000082496, 1.0000000000081788, 1.000000000008175, 1.0000000000081695, 2.000000000008132, 2.0000000000080687, 2.000000000007995, 2.0000000000079723, 2.000000000007962, 4.00000000000792, 4.000000000007874, 4.000000000007791, 4.000000000007789, 4.000000000007739, 6.000000000007687, 6.000000000007629, 6.0000000000075495, 6.000000000007532, 6.0000000000075175, 9.000000000007462, 9.000000000007406, 9.000000000007285, 9.000000000007281, 9.000000000007224, 13.000000000007182]
The results speak for themselves, though of course there will be runtime differences depending on implementation, hardware, etc.
## Time computing c_100 using both methods, and print the value each method computes for c_100
k=100
N=Nval(k)
result = %timeit -o fft_c(f, k, N=N)
print(f'FFT Best Time : {result.best}s, Computed Value = {round(fft_c(f, k, N=N))}')
result = %timeit -o c(k)
print(f'Naive Best Time : {result.best}s, Computed Value = {c(k)}')
317 µs ± 1.78 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) FFT Best Time : 0.0003160092910000003s, Computed Value = 242 16.3 ms ± 54.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each) Naive Best Time : 0.01625790417000001s, Computed Value = 242
In theory, we should have $c_n/n^3 \to 1/7500$ as $n\to \infty$. Is this what the numerics suggest? Indeed, it appears to be so, but the convergence is very slow.
## approximate c_n/n^3
k=15000
N=Nval(k)
theory = 1/7500
practice = fft_c(f,k,N=N,R=0.999)/(k**3)
print(f'Theoretical value={theory}, experimental={practice}')
Theoretical value=0.00013333333333333334, experimental=0.00013393411875637439
The thing that blows my mind, and the reason I am even writing this post, is that we have used the Fourier transform to quickly and accurately compute a combinatorial identity. It seems totally absurd to me that the Fourier transform would ever appear in an area of mathematics such as this. Moreover, this technique generalizes beyond just the sequence $\{c_n\}$. Another interesting generating function corresponds to the sequence of even zetas, i.e., $\zeta(2n) := \sum_{k\geq 1} k^{-2n}$: $$\frac{1 - \pi z \cot(\pi z)}{2} = \sum_{n\geq 1} \zeta(2n) z^{2n}$$ for $|z| < 1$. From this we can rapidly numerically compute the first several even zetas:
## Even zetas.
# check if it works
g = lambda z : (1 - np.pi * z / np.tan(np.pi * z)) / 2
k = 20
N = Nval(k)
# We'll zero out the odd values of k because they are 0 in the generating function
print(f'First {k} even zetas: {[fft_c(g, j, N=N) if j % 2 == 0 else 0 for j in range(1,k+1)]}')
First 20 even zetas: [0, 1.6449340668482246, 0, 1.0823232337111361, 0, 1.0173430619844472, 0, 1.0040773561979421, 0, 1.000994575127816, 0, 1.0002460865533058, 0, 1.0000612481350564, 0, 1.0000152822594062, 0, 1.0000038172932624, 0, 1.000000953962031]
Note that $\zeta(2) = \pi^2/6 \approx 1.6449340668482264$ and $\zeta(4) = \pi^4/90 \approx 1.082323233711138$, so our FFT method correctly predicts the value of $\zeta(2)$ and $\zeta(4)$ to at least the first 16 and 13 decimal points, respectively. The possible choices of generating function are essentially endless. Have fun playing with different generating functions. On your own, you can try calculating the Bernoulli numbers $B_n$ using the generating function $$\frac{z}{e^z - 1} = \sum_{n\geq 0} \frac{B_n}{n!}z^n$$ for $|z| < 2\pi$.
In short, we saw how we could use a generating function to compute the asymptotics of a sequence $\{c_n\}$ which is difficult to explicitly compute without tedious combinatorics. We also saw, mind blowingly, how Fourier transforms can rapidly be efficiently be used to compute the values of the sequence. One of my questions still is how to guarantee rounding accuracy for the FFT method. Values for $N$ were essentially guesses based off of cursory internet searches, but I personally have not put in the effort yet to make sure the choices of $N$ were enough to guarantee adequate numerical precision. Perhaps that will be the topic of a future blog post...