Multiplying polynomials

(1 + 2x + 3x²)(4 + 5x): three terms times two terms. Watch the six products fill a grid and collect into four coefficients.

⋯3static int products = 0;4 5void poly_mul(const int *a, int na, const int *b, int nb, int *c)6{7    for (int k = 0; k < na + nb - 1; k++)8        c[k] = 0;9    for (int i = 0; i < na; i++) {10        for (int j = 0; j < nb; j++) {11            c[i + j] += a[i] * b[j];12            products++;13        }14    }15}16 17int main(void)18{19    const int a[3] = { 1, 2, 3 };20    const int b[2] = { 4, 5 };21    int c[4];22    poly_mul(a, 3, b, 2, c);23    printf("c:");24    for (int k = 0; k < 4; k++)25        printf(" %d", c[k]);26    printf(", %d products\n", products);27    return 0;28}

Every term times every term

To multiply two polynomials, multiply every term of one by every term of the other and collect the products by power. For a = 1 + 2x + 3x² and b = 4 + 5x, whose terms are 4 and 5x:

Six products, but only four powers: the two x terms add to 13x, and the two x² terms to 22x². The product is 4 + 13x + 22x² + 15x³. It does not have m · n terms: like terms collect.

The code and its cost

The double loop

With dense coefficient arrays, term i of a times term j of b lands on power i + j:

void poly_mul(const int *a, int na, const int *b, int nb, int *c)
{
    for (int k = 0; k < na + nb - 1; k++)
        c[k] = 0;
    for (int i = 0; i < na; i++) {
        for (int j = 0; j < nb; j++) {
            c[i + j] += a[i] * b[j];
            products++;
        }
    }
}

The product of a degree-(na − 1) and a degree-(nb − 1) polynomial has degree na + nb − 2, so na + nb − 1 coefficients. The program prints

c: 4 13 22 15, 6 products

m terms times n terms take m · n products.

A convolution

Collect the products by power instead of by row: coefficient k gathers every a[i] · b[j] with i + j = k, that is

c[k] = Σᵢ a[i] · b[k − i],

over the i for which both indices are in range. That is the convolution of the two coefficient arrays: multiplying polynomials and filtering a signal are the same computation. In the grid, each sum runs along an anti-diagonal. The Toeplitz matrices lesson writes the same sum as a matrix product.

Sparse polynomials

For polynomials stored as lists of terms, multiply b by one term of a at a time, and add the partial products with the merge of adding polynomials:

struct pnode *poly_mul(const struct pnode *a, const struct pnode *b)
{
    struct pnode *result = NULL;
    for (const struct pnode *p = a; p != NULL; p = p->next)
        result = poly_add(result, times_term(b, p->coeff, p->exp));
    return result;
}

times_term builds a new list, b’s terms each multiplied by one term: the coefficients multiply and the exponents add. The program poly-mul-sparse.c multiplies (x¹⁰⁰ + 1)(x¹⁰⁰ − 1):

product: +1x^200 -1x^0

That takes 4 products, and the two x¹⁰⁰ terms cancel inside poly_add. With dense arrays of 101 coefficients each, the double loop would make 101 × 101 = 10,201 products, nearly all of them 0.

Faster methods

For dense polynomials of n coefficients the double loop takes n² products. Two faster methods exist. Split each polynomial into a low half and a high half; the product needs the three half-size products lo × lo, hi × hi and (lo + hi) × (lo + hi), instead of four, and repeating the split down to single terms takes about n^1.585 products. This is Karatsuba’s method, a later lesson. Or evaluate both polynomials with the FFT at as many points as the product has coefficients, 2n − 1, rounded up to a power of two in practice, multiply the values, and convert back, about n log n steps: fast convolution.

Running it yourself

gcc -std=c11 -Wall -Wextra -pedantic poly-mul.c -o poly-mul
./poly-mul
gcc -std=c11 -Wall -Wextra -pedantic poly-mul-sparse.c -o poly-mul-sparse
./poly-mul-sparse

Worked example

A square. (1 + x + x² + x³)² makes 4 × 4 = 16 products. The anti-diagonals of the 4 × 4 grid have 1, 2, 3, 4, 3, 2 and 1 cells, so the coefficients are 1 2 3 4 3 2 1.

One coefficient. The coefficient of x² in the program’s product: a[0] · b[2] + a[1] · b[1] + a[2] · b[0] = 1 · 0 + 2 · 5 + 3 · 4 = 22, since b has no x² term.

Large polynomials. Two polynomials of degree 1,023 have 1,024 coefficients each: 1,024² = 1,048,576 products by the double loop. Splitting in halves ten times, down to single terms, takes 3¹⁰ = 59,049 products.

Two checks. How many coefficients does a degree-4 by degree-2 product have? (7.) How many products for 5 terms by 3 terms? (15.)

Where it’s used

Computer algebra. Multiplying symbolic polynomials, with term lists for sparse ones.

Large integers. A number’s digits are the coefficients of a polynomial evaluated at x = 10; multiplying very large integers is multiplying polynomials and carrying, with the fast methods for numbers of thousands of digits.

Filtering. Filtering a signal is the same convolution (convolution), done fast with the FFT for long filters (fast convolution).

Reference card

Back to adding polynomials.

Platform: C11, gcc on x86-64 Linux: int 4 bytes, double 8, pointers 8, little-endian. Addresses in the pictures are short stand-ins, by one rule for the whole course; sizes are the real ones.