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]

TODO

TODO

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!