Naive Divide & Conquer: Polynomial Multiplication [WIP]

dc_poly_mult.py

import math

import numpy as np


def dc_poly_mult(A, B, n, a, b):
    """
    A(x) = D1(x) * x^n/2 + D0(x)
    where
    D1(x) = a_n-1 * x^(n/2-1) + a_n-2 * x^(n/2-2) + ... + a_(n/2)
    D0(x) = a_(n/2-1) * x^(n/2-1) + a_(n/2-2) * x^(n/2-2) + ... + a_0

    B(x) = E1(x) * x^n/2 + E0(x)
    where
    E1(x) = b_n-1 * x^(n/2-1) + b_n-2 * x^(n/2-2) + ... + b_(n/2)
    E0(x) = b_(n/2-1) * x^(n/2-1) + b_(n/2-2) * x^(n/2-2) + ... + b_0
    AB = (D1 * x^n/2 + D0) * (E1 * x^n/2 + E0)
       = (D1 * E1) * x^n + (D1 * E0 + D0 * E1) * x^(n/2) + D0 * E0
    """
    C = np.zeros(2 * n - 1, dtype=int)
    if n == 1:
        return A[a] * B[b]
    C[0 : n - 1] = dc_poly_mult(A, B, math.ceil(n / 2), a, b)  # product[0..n-2]
    C[n : 2 * n - 1] = dc_poly_mult(
        A, B, n // 2, a + (n // 2), b + (n // 2)
    )  # product[n..2n-2]
    D0E1 = dc_poly_mult(A, B, n // 2, a, b + (n // 2))
    D1E0 = dc_poly_mult(A, B, n // 2, a + (n // 2), b)
    C[(n // 2) : (n + n // 2 - 1)] += D0E1 + D1E0
    return C