One of the really useful properties of the Fourier transform is that convolution in one domain becomes multiplication in the other and vice-versa.
But it’s also interesting that the way we multiply regular base-10 numbers together is really a form of convolution. The decimal number ABC multiplied by DEF can be calculated by looking at C*(DEF) + B0*(DEF) + A00*(DEF), which is essentially convolving ADC with DEF (as this paper points out). If multiplication is a convolution, then we can convert to the Fourier domain to use multiplication to perform our convolution to perform our multiplication!
Here’s an example in python that takes two lists of digits and returns a new list representing the two numbers multiplied together:
import numpy as np
fft = np.fft.fft
ifft = np.fft.ifft
def mult_lists(a, b):
"""
Multiply the arrays of digits a,b by convolving a with b via a
Fourier transform.
A and b should have units in the first element, then base**1, base**2, etc.
Eg 3,412 (in base 10) would be [2, 1, 4, 3]
"""
out_size = len(a) + len(b)
a = a + [0] * (out_size - len(a))
b = b + [0] * (out_size - len(b))
f = ifft(fft(a)*fft(b)).real
return f
To use this we can write some helper functions to convert numbers to and from lists and we then have a bigint multiplication routine:
import random
def to_digit_list(num, base=10):
ret = []
while num > 0:
ret.append(num % base)
num //= base
return ret
def from_digit_list(nums, base=10):
return sum(int(round(num)) * base**i for i,num in enumerate(nums))
def fourier_mult(a, b):
return from_digit_list(mult_lists(to_digit_list(a), to_digit_list(b)))
a=random.getrandbits(256)
b=random.getrandbits(256)
ab = a*b
abf = fourier_mult(a, b)
if ab - abf == 0:
print("Numbers are the same!")
There’s something else we can do with this. We can think of each row of Pascal’s Triangle as the previous row convoluted with (1, 1). This provides us with an easy way to calculate a row of Pascal’s triangle using Fourier transforms to do the convolution for us:
def pascal_row(n):
x = [1,1] + [0] * (n-1)
fx = fft(x)
return np.round(ifft(fx ** n).real)
Rounding errors mean that only rows up to 50 or so get calculated accurately.
It’s interesting that Wikipedia talks about obtaining an alternating-sign row of Pascal’s triangle by taking the Fourier transform of sin(x)^n / x. I expect finding the right combinations of convolutions and multiplications might be able to convert between the two expressions.
[I spent some time trying this, but I found analytically finding something whose inverse Fourier transform was [1,1,0…] not that straight forward - a sinc function corresponding to a top-hat width of two shifted by half a pixel didn’t give me a particularly clean result and started giving errors for the very first rows of Pascal’s triangle. I put my attempt on github just in case anyone’s slightly interested…]
.png)
