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 implemented.

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 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:

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})$$$. Additionally, if we substite $$$A(x)$$$ with $$$P(x)$$$, then the formula tells us how to calculate $$$P(x) \% (x^{n/2} - \sqrt{c})$$$ and $$$P(x) \% (x^{n/2} + \sqrt{c})$$$ using $$$P(x) \% (x^n - c)$$$. With this we have entire recipie for implementing a $$$O(n \log n)$$$ divide and conquer algorithm.

"""
Calculates P(x) * Q(x) % (x^n - c) in O(n log n) time

Input:
  n: integer, needs to be power of 2
  c: Non-zero complex floating point number
  P: A list of length n representing a polynomial P(x)
  Q: A list of length n representing a polynomial Q(x)
Output:
  A list of length n representing the polynomial P(x) * Q(x) % (x^n - c)
"""
def fast_polymult(P, Q, n, c):
    assert len(P) == n and len(Q) == n
    
    # Base case
    if n == 1:
        return [P[0] * Q[0]]

    assert n % 2 == 0
    import cmath
    sqrtc = cmath.sqrt(c)

    # Calulate P_minus := P mod (x^(n/2) - sqrt(c))
    #          Q_minus := Q mod (x^(n/2) - sqrt(c))

    P_minus = [p1 + sqrtc * p2 for p1,p2 in zip(P[:n//2], P[n//2:])]
    Q_minus = [q1 + sqrtc * q2 for q1,q2 in zip(Q[:n//2], Q[n//2:])]

    # Calulate P_plus := P mod (x^(n/2) + sqrt(c))
    #          Q_plus := Q mod (x^(n/2) + sqrt(c))

    P_plus = [p1 - sqrtc * p2 for p1,p2 in zip(P[:n//2], P[n//2:])]
    Q_plus = [q1 - sqrtc * q2 for q1,q2 in zip(Q[:n//2], Q[n//2:])]

    # Recursive calculate product
    PQ_minus = polymult_mod_recursive(P_minus, Q_minus, n//2, sqrtc)
    PQ_plus  = polymult_mod_recursive(P_plus,  Q_plus,  n//2, -sqrtc)

    # Calculate result PQ mod (x^n - c)
    PQ = [(m + p)/2         for m,p in zip(PQ_minus, PQ_plus)] +
         [(m + p)/(2*sqrtc) for m,p in zip(PQ_minus, PQ_plus)]
    
    return PQ

One final thing that I want to note is that this algorithm can 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 divide and conquer method from earlier.

"""
Calculates P(x) * Q(x)

Input:
  P: A list representing a polynomial P(x)
  Q: A list representing a polynomial Q(x)
Output:
  A list representing the polynomial P(x) * Q(x) % (x^n - c)
"""
def polymult(P, Q)
    # Calculate length of the list representing P*Q
    n1 = len(P)
    n2 = len(Q)
    res_len = n1 + n2 - 1
    
    # Pick n sufficiently big
    n = 1
    while n < n1 + n2 - 1:
        n *= 2

    # Pad with extra 0s to reach length n
    P = P + [0] * (n - n1)
    Q = Q + [0] * (n - n2)
    
    # Pick non-zero c arbitrarily =)
    c = 123.24

    # Calculate P*Q mod x^n - c
    PQ = polymult_mod_recursive(P, Q, n, c)

    # Remove extra 0 padding and return
    return PQ[:res_len]

(Advanced) Speeding up the algorithm

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

(Advanced) Connection between the algorithm and FFT

This algorithm is actually FFT in disguise. But it is not the same algorithm as Cooley–Tukey (standard FFT algorithm).

(Advanced) Using integers instead of complex numbers

The algorithm requires two things. Firstly it needs to be possible to divide by 2, and secondly sqrt(c) needs to exist. This means that the algorithm can be extended to working modulo a prime (or modulo a odd composite number) instead of using complex numbers.