Polynomials & Interpolation

Why polynomials, at all

A polynomial is one of the simplest mathematical objects that can encode a large amount of structured information. It is just a sum of powers of a variable multiplied by coefficients:

That looks simple, but polynomials are incredibly useful because we can move between two different views of the same object:

This coefficient/value duality is everywhere in a proving system. A prover may construct a polynomial in coefficient form, evaluate it over a structured domain, commit to those evaluations, and later reason about the original polynomial again.

What exactly is a polynomial

Over a field F, a polynomial in one variable is an expression

The numbers a₀, a₁, … are the coefficients. The variable x is an indeterminate. The important distinction is that a polynomial is the formal object itself; plugging a particular value into it produces an evaluation.

For example:

The degree is the largest exponent whose coefficient is non-zero. Therefore the polynomial above has degree 2.

One polynomial, several representations

Consider:

Its coefficient representation is simply:

But we could instead store values. For example:

For a degree-2 polynomial, three distinct evaluations are enough to recover it. This is not an accident. A polynomial of degree at most d has d+1 degrees of freedom, represented by its d+1 coefficients.

Polynomial arithmetic

Polynomials over a field inherit familiar arithmetic operations.

Addition

Add corresponding coefficients:

Multiplication

Multiply every term by every term, then combine equal powers:

This is convolution of coefficient vectors. That observation becomes important later because FFTs can accelerate polynomial multiplication.

Evaluation

Given x = r, evaluate by substituting:

Naively this computes powers separately. Horner's method rewrites the same polynomial as nested multiplication:

That reduces evaluation to roughly one multiplication and one addition per coefficient.

Roots and the factor theorem

A value r is a root or zero of f if:

If r is a root, then (x−r) divides the polynomial.

For example:

This simple relationship between roots and factors is one of the reasons polynomial identities are so powerful in proof systems.

Polynomial division and the remainder

Polynomial division behaves like integer division. Given polynomials f(x) and non-zero g(x), we can write:

The polynomial q is the quotient and r is the remainder.

If we divide by x−a, the remainder is exactly f(a):

Therefore:

In ZK systems this style of reasoning appears constantly. A claim that a polynomial vanishes on a domain can be converted into a divisibility statement by a polynomial that vanishes on that same domain.

What is interpolation?

Now reverse the problem.

Instead of starting with coefficients and asking for values, suppose someone gives us values:

where the xᵢ are distinct. We want to construct the unique polynomial of degree at most d satisfying:

That process is polynomial interpolation.

Lagrange interpolation from first principles

The cleanest derivation starts by constructing special polynomials. For every point xᵢ, we want a polynomial that is 1 at xᵢ and 0 at every other sample point.

Define the Lagrange basis polynomial:

Why does this work?

Evaluate it at x=xᵢ. Every numerator becomes xᵢ−xⱼ, exactly matching its denominator:

Now evaluate it at some other sample point x=xₖ, where k≠i. One numerator becomes xₖ−xₖ=0:

So each basis polynomial behaves like a selector: it selects exactly one input point.

Because we have a basis polynomial for each data point, we can combine them with the desired y-values:

At x=xₖ, every basis term disappears except Lₖ. Therefore:

Lagrange example: three points

Take the points:

There are three points, so we seek a polynomial of degree at most 2.

The first basis polynomial is:

The second:

The third:

Combine them with their y-values:

Check:

Why is the interpolating polynomial unique?

Suppose two degree-≤d polynomials, f and g, agree on d+1 distinct points.

Consider their difference:

Then h has degree at most d, but it has at least d+1 distinct roots because f and g agree at every sample point.

A non-zero polynomial of degree d over a field cannot have more than d roots. Therefore h must be the zero polynomial:

That is the mathematical reason interpolation is unique.

Interpolation over a finite field

Everything above works over a field, not just over the real numbers. That means we can interpolate directly inside a finite field.

For example, in F₁₇, division means multiplication by a modular inverse. Suppose we need:

Since 4·13 = 52 ≡ 1 (mod 17), the inverse of 4 is 13:

Therefore the denominator in the Lagrange formula is not a floating-point division. It is exact field arithmetic.

Newton interpolation

Lagrange interpolation is conceptually clean, but there is another useful representation: Newton's form.

The coefficients cᵢ are built from divided differences. For example:

For three points:

Newton form is particularly useful when points are added incrementally because a new point can extend the representation without rebuilding every previous basis polynomial from scratch.

Why the choice of points matters

Interpolation requires distinct x-coordinates. In Lagrange interpolation, the denominators contain:

If two points have the same x-coordinate, that difference is zero and has no multiplicative inverse.

In a finite field, this is an exact algebraic requirement. There is no concept of “almost equal” points that can rescue the denominator.

This leads directly to an important ZK engineering decision: choose evaluation domains with convenient algebraic structure, often multiplicative subgroups generated by roots of unity.

Vanishing polynomials

Suppose we have a set of distinct points:

There is a natural polynomial that vanishes at every point of the domain:

For every xᵢ ∈ D, one factor becomes zero, so:

For a multiplicative subgroup generated by a primitive n-th root of unity, the domain is:

and the vanishing polynomial becomes especially simple:

This is one of the most important bridges between the previous topic and this one: roots of unity create the domain, and the domain has a compact polynomial that vanishes on every point.

Why polynomials are everywhere in ZK

A proving system needs to turn a computation into something algebraically checkable. Polynomials provide exactly that abstraction.

The conceptual chain is:

Interpolation and FFT: the important connection

Naive interpolation is mathematically straightforward but can be expensive. If we have n points, constructing all Lagrange basis polynomials directly can require roughly quadratic work.

When the interpolation points form a multiplicative subgroup:

we can exploit the structure of the domain and use the FFT / inverse FFT.

Instead of treating interpolation as a generic problem, the FFT recursively exploits the fact that powers of a root of unity form smaller root-of-unity domains.

This is why the previous page on groups and roots of unity naturally leads into this topic: the group gives the domain structure, and interpolation tells us how values on that domain correspond to a polynomial.

From the math to Rust

A minimal implementation can start with a polynomial represented by its coefficient vector. The coefficient at index i represents the coefficient of xⁱ.

// f(x) = 2 + 3x + x^2
struct Polynomial {
    coeffs: Vec<u64>,
}

impl Polynomial {
    // Horner evaluation.
    fn evaluate(&self, x: u64, modulus: u64) -> u64 {
        let mut acc = 0;

        for &coeff in self.coeffs.iter().rev() {
            acc = (acc * x + coeff) % modulus;
        }

        acc
    }
}

fn main() {
    let f = Polynomial {
        coeffs: vec![2, 3, 1],
    };

    assert_eq!(f.evaluate(2, 17), 12);
}

// Interpolate points (x_i, y_i) in F_p.
// Every x_i must be distinct.

fn mod_pow(mut base: u64, mut exp: u64, p: u64) -> u64 {
    let mut result = 1;

    while exp > 0 {
        if exp & 1 == 1 {
            result = result * base % p;
        }

        base = base * base % p;
        exp >>= 1;
    }

    result
}

fn inv(a: u64, p: u64) -> u64 {
    // Fermat: a^(p-2) = a^(-1) mod p.
    mod_pow(a, p - 2, p)
}

// Returns coefficients in ascending degree order.
fn lagrange(xs: &[u64], ys: &[u64], p: u64) -> Vec<u64> {
    assert_eq!(xs.len(), ys.len());

    let n = xs.len();
    let mut result = vec![0; n];

    for i in 0..n {
        let mut basis = vec![1u64];

        // Build ∏_{j != i} (x - x_j).
        for j in 0..n {
            if i == j {
                continue;
            }

            let xj = xs[j];
            let mut next = vec![0; basis.len() + 1];

            for k in 0..basis.len() {
                next[k] = (next[k] + p - basis[k] * xj % p) % p;
                next[k + 1] = (next[k + 1] + basis[k]) % p;
            }

            basis = next;
        }

        // Divide by ∏_{j != i} (x_i - x_j).
        let mut denom = 1;
        for j in 0..n {
            if i != j {
                denom = denom * ((xs[i] + p - xs[j]) % p) % p;
            }
        }

        let scale = ys[i] * inv(denom, p) % p;

        for k in 0..basis.len() {
            result[k] = (result[k] + scale * basis[k]) % p;
        }
    }

    result
}

The second implementation follows the mathematics almost line by line:

1. Build the numerator polynomial ∏(x−xⱼ) for each interpolation point.
2. Compute its denominator ∏(xᵢ−xⱼ).
3. Invert that denominator inside the field.
4. Scale the basis polynomial by yᵢ.
5. Add all scaled basis polynomials together.

Naive interpolation vs. production techniques

The direct Lagrange algorithm is excellent for understanding the mathematics, but a production prover cannot blindly use it for millions of points.

The engineering lesson is not that Lagrange is “bad”. It is that the representation and domain determine which algorithm is appropriate.

Polynomial identities and why they are useful in proofs

Suppose two polynomials f and g are claimed to be equal. Define:

The claim f=g is equivalent to h=0. If h is non-zero and has degree at most d, it cannot vanish at more than d points over a field.

This gives the basic intuition behind polynomial identity testing: checking values at carefully chosen or random points can provide strong evidence that two polynomials are the same, provided the field and degree bounds are handled correctly.

Complete Rust implementation

Now let's turn the mathematics into a small, self-contained Rust implementation. This version works over the prime field F₁₇ and implements the core operations we have derived: modular arithmetic, polynomial construction, evaluation, addition, multiplication, and Lagrange interpolation.

// Polynomial arithmetic and interpolation over F_17.
// Coefficients are stored in ascending degree order:
// [a0, a1, a2] represents a0 + a1*x + a2*x^2.

const MODULUS: u64 = 17;

fn mod_add(a: u64, b: u64) -> u64 {
    (a + b) % MODULUS
}

fn mod_sub(a: u64, b: u64) -> u64 {
    (a + MODULUS - b) % MODULUS
}

fn mod_mul(a: u64, b: u64) -> u64 {
    (a * b) % MODULUS
}

fn mod_pow(mut base: u64, mut exp: u64) -> u64 {
    let mut result = 1;

    while exp > 0 {
        if exp & 1 == 1 {
            result = mod_mul(result, base);
        }

        base = mod_mul(base, base);
        exp >>= 1;
    }

    result
}

fn mod_inv(a: u64) -> u64 {
    assert!(a != 0, "division by zero");

    // Fermat's little theorem:
    // a^(-1) = a^(p-2) mod p for prime p.
    mod_pow(a, MODULUS - 2)
}

#[derive(Debug, Clone)]
struct Polynomial {
    coeffs: Vec<u64>,
}

impl Polynomial {
    fn new(coeffs: Vec<u64>) -> Self {
        Self {
            coeffs: coeffs
                .into_iter()
                .map(|x| x % MODULUS)
                .collect(),
        }
    }

    // Remove unnecessary zero coefficients.
    fn normalize(&mut self) {
        while self.coeffs.len() > 1
            && *self.coeffs.last().unwrap() == 0
        {
            self.coeffs.pop();
        }
    }

    // Degree of the polynomial.
    fn degree(&self) -> usize {
        self.coeffs.len() - 1
    }

    // Evaluate f(x) using Horner's method.
    fn evaluate(&self, x: u64) -> u64 {
        let mut result = 0;

        for &coeff in self.coeffs.iter().rev() {
            result = mod_add(mod_mul(result, x), coeff);
        }

        result
    }

    // Polynomial addition.
    fn add(&self, other: &Self) -> Self {
        let n = self.coeffs.len().max(other.coeffs.len());
        let mut coeffs = vec![0; n];

        for i in 0..n {
            let a = self.coeffs.get(i).copied().unwrap_or(0);
            let b = other.coeffs.get(i).copied().unwrap_or(0);

            coeffs[i] = mod_add(a, b);
        }

        Self::new(coeffs)
    }

    // Polynomial multiplication.
    // Coefficients are multiplied by convolution.
    fn mul(&self, other: &Self) -> Self {
        let mut coeffs =
            vec![0; self.coeffs.len() + other.coeffs.len() - 1];

        for i in 0..self.coeffs.len() {
            for j in 0..other.coeffs.len() {
                let product =
                    mod_mul(self.coeffs[i], other.coeffs[j]);

                coeffs[i + j] =
                    mod_add(coeffs[i + j], product);
            }
        }

        Self::new(coeffs)
    }
}

// Lagrange interpolation.
// Given distinct x_i and corresponding y_i, return
// the unique polynomial f of degree < n such that
// f(x_i) = y_i for every i.
fn lagrange_interpolate(xs: &[u64], ys: &[u64]) -> Polynomial {
    assert_eq!(xs.len(), ys.len());
    assert!(!xs.is_empty());

    let mut result = Polynomial::new(vec![0]);

    for i in 0..xs.len() {
        // numerator = ∏_{j != i} (x - x_j)
        let mut basis = Polynomial::new(vec![1]);

        for j in 0..xs.len() {
            if i == j {
                continue;
            }

            let factor = Polynomial::new(vec![
                mod_sub(0, xs[j]),
                1,
            ]);

            basis = basis.mul(&factor);
        }

        // denominator = ∏_{j != i} (x_i - x_j)
        let mut denominator = 1;

        for j in 0..xs.len() {
            if i == j {
                continue;
            }

            denominator = mod_mul(
                denominator,
                mod_sub(xs[i], xs[j]),
            );
        }

        // Scale L_i(x) by y_i / denominator.
        let scale = mod_mul(
            ys[i],
            mod_inv(denominator),
        );

        let scaled_basis = Polynomial::new(
            basis.coeffs
                .iter()
                .map(|&c| mod_mul(c, scale))
                .collect(),
        );

        result = result.add(&scaled_basis);
    }

    let mut result = result;
    result.normalize();
    result
}

fn main() {
    // f(x) = 1 + x + x^2 over F_17.
    let f = Polynomial::new(vec![1, 1, 1]);

    assert_eq!(f.degree(), 2);
    assert_eq!(f.evaluate(0), 1);
    assert_eq!(f.evaluate(1), 3);
    assert_eq!(f.evaluate(2), 7);

    // Recover the same polynomial from its evaluations.
    let xs = [0, 1, 2];
    let ys = [1, 3, 7];

    let recovered =
        lagrange_interpolate(&xs, &ys);

    assert_eq!(recovered.coeffs, vec![1, 1, 1]);

    // Verify interpolation at additional field points.
    for x in 0..MODULUS {
        assert_eq!(
            f.evaluate(x),
            recovered.evaluate(x)
        );
    }

    println!("f(x) = 1 + x + x^2 over F_17");
    println!("interpolation verified");
}

What the implementation is actually doing

1. Represent the polynomial: [a₀, a₁, a₂, ...] stores coefficients in increasing degree order.
2. Evaluate: Horner's method evaluates the polynomial in O(n) field operations.
3. Add: coefficients with the same degree are added.
4. Multiply: every coefficient is multiplied with every other coefficient, producing convolution.
5. Construct Lagrange bases: each basis is ∏(x−xⱼ) for all j ≠ i.
6. Normalize the basis: divide by ∏(xᵢ−xⱼ), implemented as multiplication by a field inverse.
7. Scale by yᵢ: this makes the basis contribute exactly the desired value at xᵢ.
8. Sum the bases: the result is the unique polynomial passing through all supplied points.

The whole idea in one map

If finite fields answer “where does the arithmetic happen?”, groups and roots of unity answer “what structured points can we use?”, and polynomials answer “how do we encode and manipulate the computation algebraically?” Interpolation is the bridge between the polynomial's coefficients and its values.