Note that \begin{equation} A(x) = \frac{1}{2} (1 + \frac{x^{n/2}}{\sqrt{c}}) A(x) + \frac{1}{2} (1 — \frac{x^{n/2}}{\sqrt{c}}) A(x). \end{equation}

Let $$$Q^-(x)$$$ denote the quotient of $$$A(x)$$$ divided by $$$(x^n/2 - \sqrt{c})$$$ and let $$$Q^+(x)$$$ denote the quotient of $$$A(x)$$$ divided by $$$(x^n/2 + \sqrt{c})$$$. Then

and

With this we have shown that

Here $$$A(x)$$$ is expressed as remainder + quotient times $$$(x^n - c)$$$. So we have proven the formula.

"""
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_mod(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:])]

    # Recursively calculate PQ_minus := P * Q % (x^n/2 - sqrt(c)) 
    #                       PQ_plus  := P * Q % (x^n/2 + sqrt(c))
    
    PQ_minus = fast_polymult_mod(P_minus, Q_minus, n//2, sqrtc)
    PQ_plus  = fast_polymult_mod(P_plus,  Q_plus,  n//2, -sqrtc)

    # Calculate PQ mod (x^n - c) using PQ_minus and PQ_plus
    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

"""
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 fast_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 = fast_polymult_mod(P, Q, n, c)

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

We don't actually need the assumption that $$$n$$$ is a power of 2. If $$$n$$$ ever becomes odd during the recrsion, then we have two choices: Either fall back to a $$$O(n^2)$$$ algorithm or fall back to the unmodded $$$O(n \log{n})$$$ Polynomial multiplication algorithm.

Let us discuss the run time of falling back to the $$$O(n^2)$$$ algorithm when $$$n$$$ becomes odd. Assume that $$$n = a \cdot 2^b$$$, where $$$a$$$ is an odd integer and $$$b$$$ is an integer. Think of the recursive algorithm as having layers, one layer for each possible value of $$$n$$$. The first $$$b$$$ layers will all take $$$O(n)$$$ time each. In the $$$(b+1)$$$-th layer the value of $$$n$$$ is $$$a$$$. Using the $$$O(n^2)$$$ polynomial multiplication algorithm leads to this layer taking $$$O(n/a \cdot a^2) = O(n \cdot a)$$$ time. The final time complexity comes out to be $$$O((a + b) \, n)$$$.

"""
Calculates P(x) * Q(x) % (x^n - c) in O((a + b) * n) time, where n = a*2^b.

Input:
  n: Integer
  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_mod2(P, Q, n, c):
    assert len(P) == n and len(Q) == n
    
    # Base case (n is odd)
    if n & 1:
        # Calculate the answer in O(n^2) time
        res1 = [0] * n
        res2 = [0] * n
        for i in range(n):
            for j in range(n - i):
                res1[i + j] += P[i] * Q[j]
            for j in range(n - i, n):
                res2[i + j - n] += P[i] * Q[j]
        return [r1 + c * r2 for r1,r2 in zip(res1, res2)]

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

    # Recursively calculate PQ_minus := P * Q % (x^n/2 - sqrt(c)) 
    #                       PQ_plus  := P * Q % (x^n/2 + sqrt(c))
    
    PQ_minus = fast_polymult_mod(P_minus, Q_minus, n//2, sqrtc)
    PQ_plus  = fast_polymult_mod(P_plus,  Q_plus,  n//2, -sqrtc)

    # Calculate PQ mod (x^n - c) using PQ_minus and PQ_plus
    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

The reason why this is super useful is that it allows us to speed up the fast unmodded polynomial multiplication algorithm. As long as we are fine with $$$a$$$ being less than say $$$10$$$, then we might be able to choose a significantly smaller $$$n$$$ compared to what would be possible if we are only allowed to choose powers of two. This trick has the potential of making the fast unmodded polynomial multiplication algorithm run twice as fast.

"""
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 fast_polymult2(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
    b = 0
    alim = 10
    while alim * 2**b < res_len:
        b += 1
    a = (res_len - 1) // 2**b + 1
    n = a * 2**b

    # 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 = fast_polymult_mod2(P, Q, n, c)

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

Suppose that $$$P(x)$$$ and $$$Q(x)$$$ are two real polynomial, and that we want to calculate $$$P(x) \, Q(x)$$$. As discussed earlier, we can calculate the unmodded polynomial product by picking $$$n$$$ large enough such that $$$(P(x) \, Q(x)) \% (x^n - c) = P(x) \, Q(x)$$$ (here $$$c$$$ is any non-zero complex number), and then running the divide and conquer algorithm. But it turns out there is something smarter that we can do.

If we use $$$c = \text{i}$$$ (the imaginary unit) as the inital value of $$$c$$$, then this will allow us to pick an even smaller value for $$$n$$$. The reason for this is that if we get "overflow" from $$$n$$$ being too small, then that overflow will be placed into the imaginary part of the result $$$(P(x) \, Q(x)) \% (x^n - \text{i})$$$. This means that by using $$$c = \text{i}$$$ we are allowed to to pick $$$n$$$ as half the size compared to if we weren't using $$$c=\text{i}$$$. So this trick speeds the fast unmodded polynomial multiplication algorithm up by exactly a factor of 2.

TODO

In the tail of the recursion (i.e. when $$$n$$$ reaches 1), you are calculating $$$P(x) \, Q(x) \% (x - c)$$$, for some non-zero complex number $$$c$$$. This is infact the same thing as evaluating the polynomial $$$P(x) \, Q(x)$$$ at $$$x=c$$$. Furthermore, if you initially started with $$$c=1$$$, then the $$$c$$$ in the tail will be some $$$n$$$-th root of unity. If you analyse it more carefully, then you will see that each tail corresponds to a different $$$n$$$-th root of unity. So what the algorithm is actually doing is evaluating $$$P(x) \, Q(x)$$$ in all possible $$$n$$$-th roots of unity.

A $$$n$$$-th order FFT of a polynomial is defined as the polynomial evaluated in all $$$n$$$-th roots of unity. This means that the algorithm is infact an FFT algorithm. However, if you want to use it to calculate FFT, then make sure you order the $$$n$$$-th roots of unity according to the standard order used for FFT algorithms. The standard order is $$$\exp{(\frac{2 \pi \text{i}}{n} 0)}, \exp{(\frac{2 \pi \text{i}}{n} 1)}, ..., \exp{(\frac{2 \pi \text{i}}{n} (n-1))}$$$. To get the ordering right, you will probably need to do a "bit reversal" at the end.

TODO

This algorithm requires three properties. Firstly it needs to be possible to divide by $$$2$$$, and secondly $$$\sqrt{c}$$$ needs to exist, and thirdly it needs to be possible to divide by $$$\sqrt{c}$$$. This means that it is possible to extend the algorithm to work modulo a prime (or modulo an odd composite number) instead of using complex numbers.

One of the things I dislike about NTT is that for NTT to be defined, there needs to exists a $$$n$$$-th root of unity. Usually problems involving NTT are designed so that this is never an issue. But if you want to use NTT where it hasn't been designed to magically work, then this is a really big issue. The NTT can become undefined!

Note that this algorithm does not exactly share the same drawback of being undefined. The reason for this is that if $$$sqrt(c)$$$ doesn't exist, then the algorithm can simply choose to either switch over to a $$$O(n^2)$$$ polynomial multiplication algorithm, or fall back to an unmodded $$$O(n log n)$$$ polynomial multiplication algorithm. The implications from this is that this algorithm can do fast polynomial multiplication even if it is given a relatively bad NTT prime. I just find this property to be really cool!