Storing a polynomial

P = 3x⁴ − 2x + 5 stored twice: every coefficient in a strip, zeros included, and only its three terms in a chain.

⋯3 4struct term {5    double coeff;6    int exp;7};8 9struct pnode {10    double coeff;11    int exp;12    struct pnode *next;13};⋯49}50 51struct pnode *from_terms(const struct term *t, int k)52{53    struct pnode *head = NULL;54    for (int i = k - 1; i >= 0; i--) {55        struct pnode *n = malloc(sizeof *n);56        if (n == NULL)57            exit(1);58        n->coeff = t[i].coeff;59        n->exp = t[i].exp;60        n->next = head;61        head = n;62    }63    return head;⋯75int main(void)76{77    const double cp[5] = { 5, -2, 0, 0, 3 };78    const struct term tp[3] = { { 3, 4 }, { -2, 1 }, { 5, 0 } };79    struct pnode *lp = from_terms(tp, 3);80    printf("P: dense 5 cells, 3 terms; P(2) = %g %g %g\n",81           eval_dense(cp, 4, 2), eval_terms(tp, 3, 2), eval_list(lp, 2));82    double *cq = calloc(1001, sizeof *cq);83    if (cq == NULL)84        return 1;85    cq[0] = 1;86    cq[1000] = 1;87    const struct term tq[2] = { { 1, 1000 }, { 1, 0 } };88    struct pnode *lq = from_terms(tq, 2);89    printf("Q: dense 1001 cells, 2 terms; Q(1.001) = %.6g %.6g %.6g\n",90           eval_dense(cq, 1000, 1.001), eval_terms(tq, 2, 1.001), eval_list(lq, 1.001));91    free(cq);92    free_poly(lp);93    free_poly(lq);94    return 0;95}

Three ways to store it

A polynomial is a sum of terms cₑxᵉ. P = 3x⁴ − 2x + 5 has three terms, and its degree, the largest exponent, is 4.

Dense. An array of degree + 1 coefficients, the coefficient of xᵉ in cell e, zeros included:

const double cp[5] = { 5, -2, 0, 0, 3 };

Terms in an array. Only the non-zero terms, as (coefficient, exponent) pairs, largest exponent first:

struct term {
    double coeff;
    int exp;
};

const struct term tp[3] = { { 3, 4 }, { -2, 1 }, { 5, 0 } };

A struct term is 16 bytes: the double, the int, and 4 bytes of padding.

Terms in a list. The same terms in linked nodes:

struct pnode {
    double coeff;
    int exp;
    struct pnode *next;
};

A node is 24 bytes, in a 32-byte block. from_terms builds the list from an array of terms, working from the last term back and pushing each on the front, so the list keeps the array’s order.

The array of terms has a fixed size; the list grows and shrinks one term at a time, which adding polynomials needs.

The code and its cost

Evaluating the dense form

The dense form is evaluated by Horner’s rule: start with the top coefficient, then repeatedly multiply by x and add the next one down.

double eval_dense(const double *c, int deg, double x)
{
    double v = c[deg];
    for (int i = deg - 1; i >= 0; i--)
        v = v * x + c[i];
    return v;
}

That is deg multiplications. For P(2) the running value goes 3, 6, 12, 22, 49: 3; 3 × 2 + 0 = 6; 6 × 2 + 0 = 12; 12 × 2 − 2 = 22; 22 × 2 + 5 = 49.

Evaluating the terms

The term forms add up cₑ × xᵉ term by term, with a fast power for each xᵉ:

double eval_terms(const struct term *t, int k, double x)
{
    double v = 0.0;
    for (int i = 0; i < k; i++)
        v += t[i].coeff * power(x, t[i].exp);
    return v;
}

eval_list is the same loop, following next instead of an index. All three print the same value:

P: dense 5 cells, 3 terms; P(2) = 49 49 49

A polynomial with few terms

Q = x¹⁰⁰⁰ + 1 has degree 1,000 but only two terms. Dense, it needs 1,001 coefficients, 999 of them zero. The program asks calloc for them, 8,008 bytes, and sets two:

double *cq = calloc(1001, sizeof *cq);
cq[0] = 1;
cq[1000] = 1;

As terms it needs two struct terms, 32 bytes, or two list nodes, 64 bytes of blocks. A polynomial of degree 1,000 does not need 1,000 numbers if most of its coefficients are zero. All three forms give

Q: dense 1001 cells, 2 terms; Q(1.001) = 3.71692 3.71692 3.71692

since 1.001¹⁰⁰⁰ is about 2.71692, close to e. The three ways do their arithmetic in different orders and could differ in the last bits of a double; %.6g prints six significant digits, and those agree.

Which to use

Dense storage wins when most coefficients are non-zero: no exponents to store, and Horner’s rule takes degree steps. Term storage wins when few are: memory and work follow the number of terms, not the degree. Use an array of terms when the terms are fixed, and a list when terms come and go.

Running it yourself

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

Worked example

Horner by hand. P(2): 3, then 3 × 2 + 0 = 6, 6 × 2 + 0 = 12, 12 × 2 + (−2) = 22, 22 × 2 + 5 = 49.

Memory for Q. Dense: 1,001 × 8 = 8,008 bytes. Terms: 2 × 16 = 32 bytes. List: 2 blocks of 32 bytes, 64 bytes, plus the head pointer.

Work for Q. Horner on Q makes 1,000 multiplications, almost all by a coefficient of 0. The term form computes x¹⁰⁰⁰ with a fast power: 1000 is 1111101000 in binary, ten digits with six 1s, so the loop squares 10 times and multiplies 6 times.

Two checks. How many cells does 7x³ need in dense form? (4.) Which store suits adding a term often? (The list.)

Where it’s used

Computer algebra. Systems that manipulate formulas keep polynomials with many variables as lists of terms, since most possible terms are absent.

Filter design. A digital filter’s transfer function is a ratio of two polynomials in z⁻¹ (the z-transform); their coefficients are usually all non-zero, so they are stored densely.

Error-correcting codes. Codes on discs and in QR codes compute with polynomials whose coefficients are bits or bytes.

Reference card

Next, adding polynomials merges two lists of terms.

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.