Zippel's sparse interpolation in Python
Recovers a multivariate polynomial variable by variable from evaluations.
- File
zippel_sparse_interpolation.py, 109 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/zippel-sparse-interpolation/. The authors are listed in prov/GB4A8V.jsonl.
"""Zippel's probabilistic sparse interpolation over Z/pZ.
Input: a black box f: (Z/pZ)^n -> Z/pZ for a polynomial of degree at most d
in each variable, the number of variables n, the bound d, a prime p and a
random generator. Output: the polynomial as {exponent tuple: coefficient}.
The variables are added one at a time. After stage i the skeleton is the set
of monomials of f(x_1, ..., x_i, a_{i+1}, ..., a_n) for a random anchor a.
Stage i + 1 assumes that adding x_{i+1} keeps the same monomials in
x_1, ..., x_i, finds their coefficients at d + 1 values of x_{i+1} with
transposed Vandermonde systems, and interpolates each coefficient in x_{i+1}.
Cost model: operations in Z/pZ and calls of the black box. With t terms,
the method makes at most (d + 1)(1 + (n - 1) t) probes and
O(n d t (t + d)) field operations (both proved). The answer is correct with
probability at least 1 - n (n - 1) d t / (2 (p - 1)) over the random anchor
(proved in the proof file); otherwise terms can be missing, and no error is
reported. Colliding Vandermonde nodes are detected and only cost a retry.
"""
from __future__ import annotations
import random
def interp_dense(xs: list, ys: list, p: int) -> list:
"""Coefficients (low degree first) of the polynomial through (xs, ys) mod p."""
n = len(xs)
coef = list(ys)
for j in range(1, n): # Newton divided differences
for i in range(n - 1, j - 1, -1):
coef[i] = (coef[i] - coef[i - 1]) * pow(xs[i] - xs[i - j], p - 2, p) % p
poly = [0] # Newton form to monomial form, by Horner's rule
for i in range(n - 1, -1, -1):
# poly <- poly * (x - xs[i]) + coef[i]
nxt = [0] * (len(poly) + 1)
for k, c in enumerate(poly):
nxt[k + 1] = (nxt[k + 1] + c) % p
nxt[k] = (nxt[k] - c * xs[i]) % p
nxt[0] = (nxt[0] + coef[i]) % p
poly = nxt
return poly[:n]
def solve_vandermonde(nodes: list, w: list, p: int) -> list:
"""c with sum_m c_m nodes[m]^k = w[k] for k < T; the nodes are distinct."""
T = len(nodes)
lam = [1] # Lambda(z) = prod (z - v)
for v in nodes:
nxt = [0] * (len(lam) + 1)
for k, c in enumerate(lam):
nxt[k + 1] = (nxt[k + 1] + c) % p
nxt[k] = (nxt[k] - c * v) % p
lam = nxt
out = []
for v in nodes:
q = [0] * T # Lambda(z) / (z - v) by synthetic division
acc = 0
for k in range(T, 0, -1):
acc = (lam[k] + acc * v) % p
q[k - 1] = acc
num = sum(qk * wk for qk, wk in zip(q, w)) % p
den = 0
for c in reversed(q):
den = (den * v + c) % p
out.append(num * pow(den, p - 2, p) % p)
return out
def zippel(f, n: int, d: int, p: int, rng: random.Random) -> dict:
"""Sparse interpolation of the black box f with degree at most d in each variable."""
anchor = [rng.randrange(1, p) for _ in range(n)] # @line L1
xs = rng.sample(range(1, p), d + 1)
ys = [f([x] + anchor[1:]) for x in xs] # @line L2
skel = {(e,) + (0,) * (n - 1): c
for e, c in enumerate(interp_dense(xs, ys, p)) if c}
for i in range(1, n): # @line L3
S = list(skel)
T = len(S)
if T == 0:
return {}
while True: # @line L4
r = [rng.randrange(1, p) for _ in range(i)]
nodes = []
for m in S:
v = 1
for l in range(i):
v = v * pow(r[l], m[l], p) % p
nodes.append(v)
if len(set(nodes)) == T:
break
bs = rng.sample(range(1, p), d + 1)
vals = []
for b in bs: # @line L5
w = []
for k in range(T):
pt = [pow(r[l], k, p) for l in range(i)] + [b] + anchor[i + 1:]
w.append(f(pt))
vals.append(solve_vandermonde(nodes, w, p)) # @line L6
new = {}
for j, m in enumerate(S): # @line L7
cs = interp_dense(bs, [vals[jj][j] for jj in range(d + 1)], p)
for e, c in enumerate(cs):
if c:
new[m[:i] + (e,) + m[i + 1:]] = c
skel = new
return skel # @line L8
test_zippel_sparse_interpolation.py 131 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/zippel-sparse-interpolation/. The authors are listed in prov/GB4A8V.jsonl.
"""Tests for Zippel's sparse interpolation over Z/pZ.
Black boxes are built from known random sparse polynomials. The output must
equal the hidden polynomial, agree with it at many random points, and match
dense tensor-product interpolation on small cases. The number of probes is
checked against the bound (d + 1)(1 + (n - 1) t).
"""
import random
from itertools import product
import pytest
from zippel_sparse_interpolation import interp_dense, solve_vandermonde, zippel
P = 2_147_483_647 # 2^31 - 1, prime
def evaluate(poly, pt, p):
s = 0
for e, c in poly.items():
v = c
for x, k in zip(pt, e):
v = v * pow(x, k, p) % p
s = (s + v) % p
return s
class Box:
def __init__(self, poly, p):
self.poly, self.p, self.calls = poly, p, 0
def __call__(self, pt):
self.calls += 1
return evaluate(self.poly, pt, self.p)
def random_sparse(rng, n, d, t, p):
mons = set()
while len(mons) < t:
mons.add(tuple(rng.randint(0, d) for _ in range(n)))
return {m: rng.randrange(1, p) for m in mons}
def dense_reference(f, n, d, p):
"""Tensor-product Lagrange interpolation on the grid {1..d+1}^n."""
grid = list(range(1, d + 2))
vals = {pt: f(list(pt)) for pt in product(grid, repeat=n)}
for axis in range(n): # interpolate along one axis at a time
new = {}
for key in {k[:axis] + k[axis + 1:] for k in vals}:
ys = [vals[key[:axis] + (g,) + key[axis:]] for g in grid]
for e, c in enumerate(interp_dense(grid, ys, p)):
new[key[:axis] + (e,) + key[axis:]] = c
vals = new
return {e: c for e, c in vals.items() if c}
def test_helpers():
rng = random.Random(1)
for _ in range(50):
deg = rng.randint(0, 8)
poly = [rng.randrange(P) for _ in range(deg + 1)]
xs = rng.sample(range(1, P), deg + 1)
ys = [sum(c * pow(x, k, P) for k, c in enumerate(poly)) % P for x in xs]
assert interp_dense(xs, ys, P) == poly
T = rng.randint(1, 8)
nodes = rng.sample(range(P), T)
c = [rng.randrange(P) for _ in range(T)]
w = [sum(ci * pow(v, k, P) for ci, v in zip(c, nodes)) % P for k in range(T)]
assert solve_vandermonde(nodes, w, P) == c
def test_example():
# f = 3 x^2 y + 5 x y^3 - 2, d = 3
f = {(2, 1): 3, (1, 3): 5, (0, 0): P - 2}
assert zippel(Box(f, P), 2, 3, P, random.Random(7)) == f
def test_edges():
rng = random.Random(3)
assert zippel(Box({}, P), 3, 4, P, rng) == {} # zero polynomial
assert zippel(Box({(0, 0, 0): 9}, P), 3, 4, P, rng) == {(0, 0, 0): 9} # constant
one = {(5,): 2, (0,): 1}
assert zippel(Box(one, P), 1, 5, P, rng) == one # one variable
dense = {e: 1 + sum(e) for e in product(range(3), repeat=2)}
assert zippel(Box(dense, P), 2, 2, P, rng) == dense # all (d+1)^n terms
last = {(0, 0, 0, 4): 1} # only the last variable
assert zippel(Box(last, P), 4, 4, P, rng) == last
# a degree bound larger than the true degree is fine
f = {(1, 1): 4, (0, 2): 6}
assert zippel(Box(f, P), 2, 7, P, rng) == f
@pytest.mark.parametrize("seed", range(6))
def test_random_sparse(seed):
rng = random.Random(100 + seed)
for _ in range(10):
n, d = rng.randint(1, 5), rng.randint(1, 8)
t = rng.randint(1, min(20, (d + 1) ** n))
f = random_sparse(rng, n, d, t, P)
box = Box(f, P)
g = zippel(box, n, d, P, random.Random(seed * 1000 + n))
assert g == f
assert box.calls <= (d + 1) * (1 + (n - 1) * t)
for _ in range(20):
pt = [rng.randrange(P) for _ in range(n)]
assert evaluate(g, pt, P) == evaluate(f, pt, P)
@pytest.mark.parametrize("seed", range(3))
def test_against_dense(seed):
rng = random.Random(900 + seed)
for _ in range(10):
n, d = rng.randint(1, 3), rng.randint(1, 3)
f = random_sparse(rng, n, d, rng.randint(1, (d + 1) ** n), P)
ref = dense_reference(Box(f, P), n, d, P)
assert ref == f
assert zippel(Box(f, P), n, d, P, random.Random(seed)) == ref
def test_small_prime():
# with p = 10007 the failure bound is still small for these sizes
p = 10007
rng = random.Random(42)
for _ in range(30):
n, d = rng.randint(2, 3), rng.randint(1, 4)
f = random_sparse(rng, n, d, rng.randint(1, min(6, (d + 1) ** n)), p)
assert zippel(Box(f, p), n, d, p, random.Random(rng.randrange(10**6))) == f
What the next level needs
- the tests did not run in this build