Unit 19: Multiplying Fast

Every pair, added into its slot

Multiply 1 + 2x by 1 + 3x and you get 1 + 5x + 6x². Each term of the first meets each term of the second, and the powers add.

So the slot k of the answer is the sum of a[i] * b[k - i]. That operation is called a convolution, and it costs n * n.

Python

def schoolbook(a, b):
    result = [0] * (len(a) + len(b) - 1)
    for i in range(len(a)):
        for j in range(len(b)):
            result[i + j] += a[i] * b[j]
    return result

print(schoolbook([1, 2], [1, 3]))
print(schoolbook([1, 2, 3], [4, 5, 6]))

It prints

[1, 5, 6]
[4, 13, 28, 27, 18]

Why a contest cares

Count the pairs of cards that add up to each total. Count the ways to make each score with two spinners. Multiply two enormous numbers.

All three are the same sum — a[i] * b[k - i] — so anything that multiplies polynomials quickly does all of them quickly.

Python

first = [0, 1, 1, 0, 1]
second = [0, 1, 0, 1, 1]

result = [0] * 9
for i in range(5):
    for j in range(5):
        result[i + j] += first[i] * second[j]

for total in range(2, 9):
    print(total, result[total])

It prints

2 1
3 1
4 1
5 3
6 1
7 1
8 1

And why it has to be faster

Two polynomials of a hundred thousand terms is ten thousand million multiplications. Far too slow.

The fast way gets it to about n log n — under two million — and the rest of the unit is how.

A polynomial can be written two ways

As a list of coefficients — 1, 5, 6 — or as its value at enough points. Two points pin down a straight line; three pin down a quadratic; n pin down anything of degree below n.

Both descriptions hold exactly the same information.

Python

def value_at(poly, x):
    total = 0
    power = 1
    for c in poly:
        total += c * power
        power *= x
    return total

poly = [1, 5, 6]
for x in [0, 1, 2, 3]:
    print(x, value_at(poly, x))

It prints

0 1
1 12
2 35
3 70

And in point form, multiplying is trivial

If you know both polynomials at the same x, the product at that x is just the two values multiplied. One multiplication per point instead of n of them.

Here (1 + 2x) and (1 + 3x) are worked out at four points, multiplied point by point, and the answers match the product 1 + 5x + 6x² exactly.

Python

def value_at(poly, x):
    total = 0
    power = 1
    for c in poly:
        total += c * power
        power *= x
    return total

a = [1, 2]
b = [1, 3]
product = [1, 5, 6]

for x in [0, 1, 2, 3]:
    print(x, value_at(a, x) * value_at(b, x), value_at(product, x))

It prints

0 1 1
1 12 12
2 35 35
3 70 70

So the plan is three steps

Evaluate both polynomials at n points. Multiply the values point by point — that part is n. Interpolate: turn the points back into coefficients.

And it saves nothing at all, because evaluating at n points the obvious way costs n * n. The saving has to come from choosing the points cleverly.

Choose points that repeat when squared

Evaluating at x and at -x shares work: the even-power terms give the same answer both times, and the odd ones only change sign.

So pick points that come in plus-and-minus pairs, and whose squares are again such a set. On the number line that runs out immediately. In the complex numbers it goes on for ever.

The roots of unity

The n numbers spread evenly round a circle of radius 1, each one e to the power 2πik/n. Squaring them gives the n/2 roots of unity, twice over — exactly the folding needed.

For n = 4 they are 1, i, -1 and -i, and yes, they really are the four numbers whose fourth power is 1.

Python

import cmath

n = 4
for k in range(n):
    z = cmath.exp(2j * cmath.pi * k / n)
    print(k, int(round(z.real)), int(round(z.imag)), int(round((z ** 4).real)))

It prints

0 1 0 1
1 0 1 1
2 -1 0 1
3 0 -1 1

Evaluating at all of them: the transform

The discrete Fourier transform of a list is its polynomial's value at every root of unity. Done directly it is n * n — this is only here to say what the answer should be.

For 1, 2, 3, 4 it comes out as 10, then -2+2i, -2, -2-2i. The 10 is the plain sum, because the first root of unity is 1.

Python

import cmath

def slow_dft(a):
    n = len(a)
    out = []
    for k in range(n):
        total = 0
        for j in range(n):
            total += a[j] * cmath.exp(-2j * cmath.pi * j * k / n)
        out.append(total)
    return out

result = slow_dft([1, 2, 3, 4])
print([(int(round(z.real)), int(round(z.imag))) for z in result])

It prints

[(10, 0), (-2, 2), (-2, 0), (-2, -2)]

Now split odds from evens

Write the polynomial as its even-numbered terms plus x times its odd-numbered ones. Each half is a polynomial of half the size — and it needs evaluating at the squares of the points, which are the half-size roots of unity.

Two problems of size n/2, joined in n steps. That is n log n, and it is the fast Fourier transform.

Halve, halve, and join back up

Split the list into even and odd positions, transform each half, then combine: entry k is the even half's k plus a twiddle times the odd half's k, and entry k + n/2 is the same with a minus.

The inverse is the identical function with the angle turned the other way, and everything divided by n at the end.

Python

import cmath

def fft(a, invert):
    n = len(a)
    if n == 1:
        return a[:]
    even = fft(a[0::2], invert)
    odd = fft(a[1::2], invert)
    angle = 2 * cmath.pi / n * (1 if invert else -1)
    out = [0] * n
    for k in range(n // 2):
        t = cmath.exp(1j * angle * k) * odd[k]
        out[k] = even[k] + t
        out[k + n // 2] = even[k] - t
    return out

result = fft([complex(1), complex(2), complex(3), complex(4)], False)
print([(int(round(z.real)), int(round(z.imag))) for z in result])

It prints

[(10, 0), (-2, 2), (-2, 0), (-2, -2)]

The same four numbers as the slow one

That is the check that matters: the fast transform agrees exactly with the direct sums of the previous lesson.

It only works when the length is a power of two, because the halving has to keep working all the way down. So pad with zeroes up to one.

Multiplying, at last

Pad both lists to a power of two big enough for the answer, transform both, multiply point by point, transform back, divide by the size and round.

And the answers match schoolbook exactly — on lists this small it is slower, and on a hundred thousand terms it is the only thing that finishes.

Python

import cmath

def fft(a, invert):
    n = len(a)
    if n == 1:
        return a[:]
    even = fft(a[0::2], invert)
    odd = fft(a[1::2], invert)
    angle = 2 * cmath.pi / n * (1 if invert else -1)
    out = [0] * n
    for k in range(n // 2):
        t = cmath.exp(1j * angle * k) * odd[k]
        out[k] = even[k] + t
        out[k + n // 2] = even[k] - t
    return out

def multiply(p, q):
    size = 1
    while size < len(p) + len(q):
        size *= 2
    a = [complex(x) for x in p] + [0] * (size - len(p))
    b = [complex(x) for x in q] + [0] * (size - len(q))
    fa = fft(a, False)
    fb = fft(b, False)
    fc = [fa[i] * fb[i] for i in range(size)]
    back = fft(fc, True)
    return [int(round(z.real / size)) for z in back][:len(p) + len(q) - 1]

print(multiply([1, 2], [1, 3]))
print(multiply([1, 2, 3], [4, 5, 6]))

It prints

[1, 5, 6]
[4, 13, 28, 27, 18]

The catch, and it is a real one

Everything came back as a float and had to be rounded. With big enough numbers the rounding stops being safe, and the answer is quietly wrong.

Real contest code for exact answers uses a number-theoretic transform instead: the same algorithm with roots of unity taken under a prime modulus, so every value is a whole number and nothing is ever rounded.