Fast polynomial multiplication
- Motivation
- Problem definition
- Karatsuba algorithm
- Multiplication using polynomial interpolation
- Conclusion
- Few practical tips
- Test your understanding
Motivation
Many computer science problems can be reduced to a polynomial multiplication. The most familiar example is the integer multiplication. Have you ever wondered why multiplying 100000-digit numbers in Python takes seconds and not minutes or hours?
The school algorithm is very simple, but also very slow. Its time complexity is O(n2). There are much more efficient algorithms and the goal of this article is to explore some of them.
After formalizing the problem, we’ll have a look at the Karatsuba algorithm for multiplying polynomials in O(n1.58) time. Next, we’ll try to multiply polynomials using polynomial interpolation which will be the basis of O(nlogn) algorithms that utilize Fast Fourier Transform (FFT) and Number Theoretic Transform (NTT).
All these algorithms use the Divide and Conquer paradigm, so this article can also serve as an advanced introduction (through examples) to that problem-solving approach.
Problem definition
In the rest of the article we’ll denote a polynomial of degree n using one of the following:
p=p(x)=n∑j=0pjxj=p0+p1x+p2x2+..+pnxnThe input consists of two polynomials a and b of the same degree n (if that’s not the case we can always pad with zeroes), and the output is a polynomial c of degree 2n:
c(x)=a(x)b(x)=n∑i=0n∑j=0aibjxi+jThe school algorithm written in C is simple and unambiguously illustrates the problem:
void mul(const int* a, const int* b, int* c, int n) {
for (int i = 0; i <= 2*n; ++i) {
c[i] = 0;
}
for (int i = 0; i <= n; ++i) {
for (int j = 0; j <= n; ++j) {
c[i+j] += a[i] * b[j];
}
}
}
Karatsuba algorithm
We will solve the problem recursively:
- the base case is n=0, then simply c=a0b0
- otherwise (n>0):
- split each of the input polynomials into two smaller polynomials:
- m=⌈n/2⌉
- a=pxm+q
- b=rxm+s
- for example, if a=x3−2x2+x+4, then m=2, q=x+4 and p=x−2
- write c using new polynomials p,q,r,s:
- c=ab=(pxm+q)(rxm+s)=prx2m+(ps+qr)xm+qs
- introduce simple substitutions:
- z2=pr
- z1=ps+qr
- z0=qs
- now we have:
- c=z2x2m+z1xm+z0
- it’s clear we can compute z0,z1 and z2 using four multiplications of half-sized polynomials, but can we do better?
- if we compute z0 and z2 first, it turns out we can compute z1 as follows:
- z1=(p+q)(r+s)−pr−qs=(p+q)(r+s)−z0−z2
- split each of the input polynomials into two smaller polynomials:
One multiplication is traded for several additions (which are cheap) so now we need only three multiplications of half-degreed polynomials. We can multiply those by invoking the algorithm recursively.
The time complexity of this algorithm is O(nlog23) which is close to, and is usually written as, O(n1.58).
If you’d like to know how to compute the complexity of this and similar algorithms in an elegant way, see the Master theorem (Wikipedia).
Multiplication using polynomial interpolation
It is well known that n+1 pairs (x0,y0),..,(xn,yn) where xi≠xj for i≠j uniquely determine a polynomial p of degree n such that p(xi)=yi. We say that from evaluations of a polynomial in n+1 points we can reconstruct (interpolate) the polynomial. A simple intuition for this is that the polynomial has n+1 unknown coefficients, and each evaluation is an equation.
In the polynomial multiplication problem, we are looking for the polynomial c of degree 2n. What if we find values of c in 2n+1 points? Then we should be able to reconstruct the solution!
How to find the value of c in a point x? Since c(x)=a(x)b(x), it suffices to evaluate both polynomials in x and multiply the results.
Now we are ready to sketch out the high-level, three-step algorithm:
- Choose 2n+1 different numbers (points) x0,..,x2n and evaluate polynomials a and b in them,
- Multiply evaluations for each point: yi=a(xi)b(xi),
- Interpolate c from pairs (x0,y0),..,(x2n,y2n).
This is great, but already the first step is difficult to perform in better than O(n2) time. Lucky for us, someone has noticed that evaluation of a polynomial in many points (multipoint evaluation) can be done more efficiently if the set of points is carefully chosen. Furthermore, the third step can then also be reduced to the multipoint evaluation of a similar set of points. The second step is clearly linear, so this sounds promising.
Let’s first see how to do the multipoint evaluation efficiently.
Multipoint evaluation
The algorithm we are going to see takes a polynomial p, a natural number N and a complex number w. It returns evaluations of p in N points: w0,w1,..,wN−1.
It’s convenient to write down these N evaluations in the form of a degree N−1 polynomial which is their generating function. In that case, we can say that the algorithm returns a new polynomial, so we can also say the algorithm is a polynomial transformation. It takes a polynomial and returns a polynomial. In the rest of the article, when we say transformation we mean this transform.
Formally, the transformation is:
TN,w(p)=N−1∑j=0p(wj)xjIn general, it is hard to compute this transformation efficiently. But not if w is chosen carefully. As we’ll see in the following section, the transformation can be computed efficiently if w is an Nth root of unity, that is, if wN=1.
We’ll use a variant of Cooley-Tukey algorithm that computes this transformation when N is a power of 2.
Cooley-Tukey algorithm for multipoint evaluation
Input: N=2q, w such that wN=1 and a polynomial p.
Output: transformed polynomial p′=TN,w(p) .
Again, we solve recursively:
- the base case is N=1
- p′=p(w0)=p(1)=p0
- otherwise (N>1):
- let’s look at one element of the result p′j:
- p′j=p(wj)=N−1∑k=0wjkpk
- split on parity of powers in p, and let m=N2:
- p′j=m−1∑k=0wj2kp2k+m−1∑k=0wj(2k+1)p2k+1
- p′j=m−1∑k=0(w2)jkp2k+wjm−1∑k=0(w2)jkp2k+1
- introduce substitutions e and o, for coefficients next to even and odd powers in p, respectively:
- e=m−1∑j=0p2jxj
- o=m−1∑j=0p2j+1xj
- express p′j using e i o :
- p′j=m−1∑k=0(w2)jkek+wjm−1∑k=0(w2)jkok
- p′j=e((w2)j)+wjo((w2)j)
- here we have to recognize transformation expressions of polynomials e and o with squared w
- transform recursively e and o with w2 (note that (w2)m=1):
- e′=Tm,w2(e)
- o′=Tm,w2(o)
- now it’s clear that for 0≤j<m:
- p′j=e′j+wjo′j
- what about j≥m? It might seem impossible because we haven’t evaluated e and o in points (w2)j, but since (w2)m=1 we can use the evaluations in (w2)j−m=(w2)j.
- to summarize:
- p′j={e′j+wjo′jfor 0≤j<me′j−m+wjo′j−mfor m≤j<N
- let’s look at one element of the result p′j:
Implementing it recursively will give us an algorithm of O(NlogN) time complexity.
Cool, we can evaluate efficiently. But what about the interpolation? We can look at the interpolation as an inverse of our transformation. It turns out it is the same as the transformation with inverted w and divided by N:
T−1N,w(p′)=1NTN,w−1(p′)=1NN−1∑j=0p′(w−j)xjThe proof is left as an exercise to the reader :).
What it means is that we can use the same multipoint evaluation algorithm to do the interpolation too!
The only thing missing is the value of w…
Multiplication using Fast Fourier Transform (FFT)
Let’s use w=ei2π/N.
In that case, our transform is also known as the Discrete Fourier Transform (DFT).
If we now go back to the original multiplication problem, we can summarize the complete algorithm as follows:
c=ab=DFT−1(DFT(a)⋅DFT(b))where:
- DFT(p)=TN,ei2π/N(p),
- DFT−1(p)=1NTN,e−i2π/N(p),
- N= smallest power of two greater than 2n and
- ⋅ is element-wise multiplication operator on polynomials, i.e.: (r⋅s)i=risi.
We can use the Cooley-Tukey algorithm to compute the DFT efficiently. Such and all other fast implementations of DFT are called Fast Fourier Transformations (FFT).
The complexity of this polynomial multiplication algorithm is O(nlogn).
Multiplication using Number Theoretic Transform (NTT)
A practical disadvantage of DFT can be the fact it uses complex numbers. An alternative is to use polynomials defined over finite fields GF(p), that is, all the operations on numbers are done modulo p, where p is prime. This variant of Fourier transform is also known as Number-theoretic transform (NTT).
Similarly to DFT, the product of polynomials a and b of degree n can we written as:
c=ab=NTT−1(NTT(a)⋅NTT(b))where:
- NTT(p)=TN,g(p−1)/N(p),
- NTT−1(p)=1NTN,g−(p−1)/N(p),
- g= a primitive root modulo p (What is a primitive root and how to find one?),
- N= smallest divisor of p−1 that is greater than 2n and
- ⋅ is element-wise multiplication operator on polynomials, i.e.: (r⋅s)i=risi.
Note that (g(p−1)/N))N=1. Hence we can use the same transformation algorithm as before, the main difference being that all the computations will be done using integers and modulo p.
There is a catch though. N must divide p−1, and with the old requirement of N being a power of two, finding an appropriate N can be an impossible task.
To summarize, to apply this algorithm, p must be a prime of form k2q+1 and N must be a power of two and at most 2q.
However, it is sometimes possible to modify the Cooley-Tukey to support other values of p, by not always splitting the problem into two equal parts. Details are out of the scope of this article, but feel free to experiment.
Conclusion
We started with the polynomial multiplication problem but we also learned how to do FFT efficiently. FFT, on the other hand, is widely used for signal processing. Some big-integer libraries still use the Karatsuba algorithm, while others have opted for FFT or even fancier algorithms.
FFT/NTT-based multiplications are definitely a better choice than Karatsuba when talking about speed, but FFT can have precision problems while NTT might not be applicable. Another random benefit of the Karatsuba algorithm is that it can be used even when the division is not defined (modulo 24, for example).
Few practical tips
-
NTT-based multiplication can be used even if we need the full, “non-modulo”, result. It suffices to use a p larger than any value in the result. However, in practice we may be limited by integer data type sizes. In that case, we can repeat the multiplication using different primes. They can be small, but as long as their product is greater than any value in the result, the full result can be reconstructed using the Chinese Remainder Theorem.
-
When implementing Cooley-Tukey to do FFT be aware of the form of w and compute its powers using De Moivre’s formula.
-
To combat FFTs precision issues with big numbers we can split the polynomial’s coefficients into few smaller values. For example each coefficient ai could be expressed as ai=xiK+yi. Where K is some constant, for example K=√max(ai). If we do a similar split for the other polynomial, we can express the product as a weighted sum of 4 smaller (in terms of values) products.
-
Karatsuba algorithm will be faster than more sophisticated algorithms on smaller degree polynomials as it has a relatively low overhead factor. Similarly, the school algorithm, having the lowest overhead factor, will be the fastest for very small polynomials. So when the performance is critical, multiple approaches should be combined.
Test your understanding
This article was originally written as a competitive programming tutorial. So naturally, I’ve collected a few interesting problems available online, ranging from the ones where you just have to implement one of these algorithms, to the ones where you have to modify the algorithms and adapt them: