Tutorial: A simple O(n log n) polynomial multiplication algorithm

Revision en11, by pajenegod, 2023-07-07 13:51:58

Hi Codeforces!

I have recently come up with a really neat and simple recursive algorithm for multiplying polynomials in $$$O(n \log n)$$$ time. It is so neat and simple that I think it might possibly revolutionize the way that fast polynomial multiplication is taught and coded. You don't need to know anything about FFT to understand and implement this algorithm.

I've split this blog up into two parts. The first part is intended for anyone to be able to read and understand. The second part is advanced and goes into a ton of interesting ideas and concepts related to this algorithm.

Prerequisite: Polynomial quotient and remainder, see Wiki article and Stackexchange example.

Task:

Given two polynomials $$$P$$$ and $$$Q$$$ and an integer $$$n$$$, where degree $$$P < n$$$ and degree $$$Q < n$$$. Your task is to calculate the polynomial $$$P(x) \, Q(x) \% (x^n - c)$$$ in $$$O(n \log n)$$$ time. You may assume that $$$n$$$ is a power of two.

Solution:

We can create a divide and conquer algorithm for $$$P(x) \, Q(x) \% (x^n - c)$$$ based on the difference of squares formula. Assuming $$$n$$$ is even, then $$$(x^n - c) = (x^{n/2} - \sqrt{c}) (x^{n/2} + \sqrt{c})$$$. The idea behind the algorithm is to calculate $$$P(x) \, Q(x) \% (x^{n/2} - \sqrt{c})$$$ and $$$P(x) \, Q(x) \% (x^{n/2} + \sqrt{c})$$$ using 2 recursive calls, and then use that result to calculate $$$P(x) \, Q(x) \% (x^n - c)$$$.

So how do we actually calculate $$$P(x) \, Q(x) \% (x^n - c)$$$ using $$$P(x) \, Q(x) \% (x^{n/2} - \sqrt{c})$$$ and $$$P(x) \, Q(x) \% (x^{n/2} + \sqrt{c})$$$?

Well, we can use the following formula:

$$$ \begin{aligned} A(x) \% (x^n - c) = &\frac{1}{2} (1 + \frac{x^{n/2}}{\sqrt{c}}) (A(x) \% (x^{n/2} - \sqrt{c})) \, + \\ &\frac{1}{2} (1 - \frac{x^{n/2}}{\sqrt{c}}) (A(x) \% (x^{n/2} + \sqrt{c})). \end{aligned} $$$
Proof of the formula

This formula is very useful. If we substitute $$$A(x)$$$ by $$$P(x) Q(x)$$$, then the formula tells us how to calculate $$$P(x) \, Q(x) \% (x^n - c)$$$ using $$$P(x) \, Q(x) \% (x^{n/2} - \sqrt{c})$$$ and $$$P(x) \, Q(x) \% (x^{n/2} + \sqrt{c})$$$ in linear time. With this we have the recipie for implementing a $$$O(n \log n)$$$ divide and conquer algorithm:

Input:

  • Integer $$$n$$$ (power of 2),
  • Non-zero complex number $$$c$$$,
  • Two polynomials $$$P(x) \% (x^n - c)$$$ and $$$Q(x) \% (x^n - c)$$$.

Output:

  • The polynomial $$$P(x) \, Q(x) \% (x^n - c)$$$.

Algorithm:

Step 1. (Base case) If n = 1, then return $$$P(0) \cdot Q(0)$$$. Otherwise:

Step 2. Starting from $$$P(x) \% (x^n - c)$$$ and $$$Q(x) \% (x^n - c)$$$, in $$$O(n)$$$ time calculate

$$$ \begin{align} P(x) \% (x^{n/2} - \sqrt{c}), \\ Q(x) \% (x^{n/2} - \sqrt{c}), \\ P(x) \% (x^{n/2} + \sqrt{c}), \\ Q(x) \% (x^{n/2} + \sqrt{c}). \end{align} $$$

Step 3. Make two recursive calls to calculate $$$P(x) \, Q(x) \% (x^{n/2} - \sqrt{c})$$$ and $$$P(x) \, Q(x) \% (x^{n/2} + \sqrt{c})$$$.

Step 4. Using the formula, calculate $$$P(x) \, Q(x) \% (x^n - c)$$$ in $$$O(n)$$$ time. Return the result.

Here is a Python implementation following this recipie:

Python solution to the task

One final thing that I want to mention before going into the advanced section is that this algorithm can also be used to do fast unmodded polynomial multiplication, i.e. given polynomials $$$P(x)$$$ and $$$Q(x)$$$ calculate $$$P(x) \, Q(x)$$$. The trick is simply to pick $$$n$$$ large enough such that $$$P(x) \, Q(x) = P(x) \, Q(x) \% (x^n - c)$$$, and then use the exact same algorithm as before. $$$c$$$ can be arbitrarily picked (any non-zero complex number works).

Python implementation for general Fast polynomial multiplication

(Advanced) Speeding up the algorithm

This section will be about tricks that can be used to speed up the algorithm. This will in total speed it up by a factor of between 2 and 4.

$n$ doesn't actually need to be a power of 2
Imaginary-cyclic convolution
Trick to go from $\% (x^n - c)$ to $\% (x^n - 1)$

(Advanced) -is-this-fft-?

This algorithm is actually FFT in disguise. But it is also different compared to any other FFT algorithm that I've seen before (for example the Cooley–Tukey FFT algorithm).

Using this algorithm to calculate FFT
This algorithm is not the same algorithm as Cooley–Tukey

(Advanced) Connection between this algorithm and NTT

Just like how there is FFT and NTT, there are two variants of this algorithm too. One using complex floating point numbers, and the other using modulo a prime (or more generally modulo an odd composite number).

Using modulo integers instead of complex numbers
What if $sqrt(c)$ doesn't exist?
Tags fft, convolution, polynomials, recursion

History

 
 
 
 
Revisions
 
 
  Rev. Lang. By When Δ Comment
en53 English pajenegod 2023-07-21 20:29:32 631
en52 English pajenegod 2023-07-21 19:29:42 0 (published)
en51 English pajenegod 2023-07-21 19:29:23 120 (saved to drafts)
en50 English pajenegod 2023-07-10 22:53:22 1131
en49 English pajenegod 2023-07-10 20:27:39 204
en48 English pajenegod 2023-07-10 20:18:13 0 (published)
en47 English pajenegod 2023-07-10 20:14:48 435
en46 English pajenegod 2023-07-10 19:55:32 288
en45 English pajenegod 2023-07-10 19:51:08 13
en44 English pajenegod 2023-07-10 19:48:33 7
en43 English pajenegod 2023-07-10 19:41:29 1385
en42 English pajenegod 2023-07-10 03:42:37 2274
en41 English pajenegod 2023-07-09 22:40:43 571
en40 English pajenegod 2023-07-09 22:00:53 39
en39 English pajenegod 2023-07-09 21:57:38 1100
en38 English pajenegod 2023-07-09 20:46:05 384
en37 English pajenegod 2023-07-09 20:03:42 39
en36 English pajenegod 2023-07-09 19:36:59 3745
en35 English pajenegod 2023-07-09 17:00:26 121
en34 English pajenegod 2023-07-09 16:53:59 2297
en33 English pajenegod 2023-07-08 22:04:53 23
en32 English pajenegod 2023-07-08 15:24:23 3
en31 English pajenegod 2023-07-08 15:22:29 24
en30 English pajenegod 2023-07-08 15:19:05 192
en29 English pajenegod 2023-07-08 15:03:43 438
en28 English pajenegod 2023-07-08 14:50:40 16
en27 English pajenegod 2023-07-08 14:48:51 278
en26 English pajenegod 2023-07-08 14:41:17 5
en25 English pajenegod 2023-07-08 14:40:38 7540
en24 English pajenegod 2023-07-08 13:38:20 1
en23 English pajenegod 2023-07-08 13:27:45 24
en22 English pajenegod 2023-07-08 13:26:39 844
en21 English pajenegod 2023-07-08 13:08:19 1427
en20 English pajenegod 2023-07-08 12:36:17 492
en19 English pajenegod 2023-07-08 04:13:26 5
en18 English pajenegod 2023-07-08 04:10:18 30
en17 English pajenegod 2023-07-08 03:59:49 999
en16 English pajenegod 2023-07-08 03:42:36 977
en15 English pajenegod 2023-07-08 03:27:23 4195
en14 English pajenegod 2023-07-07 20:07:54 1204
en13 English pajenegod 2023-07-07 18:13:08 2
en12 English pajenegod 2023-07-07 18:08:18 39
en11 English pajenegod 2023-07-07 13:51:58 89
en10 English pajenegod 2023-07-07 13:41:19 789
en9 English pajenegod 2023-07-07 06:29:01 11
en8 English pajenegod 2023-07-07 02:22:13 4
en7 English pajenegod 2023-07-07 02:02:33 122 ffao,-is-this-fft-,meooow,ToxicPie9,algmyr,aryanc403,kostia244,nor,Spheniscine,magnus.hegdahl,jeroenodb
en6 English pajenegod 2023-07-07 01:52:20 119
en5 English pajenegod 2023-07-07 01:46:09 448
en4 English pajenegod 2023-07-07 01:35:15 2157
en3 English pajenegod 2023-07-07 00:53:13 1957
en2 English pajenegod 2023-07-07 00:15:55 2037
en1 English pajenegod 2023-07-06 23:46:15 3121 Initial revision (saved to drafts)