Point-value representation in Python
Stores the values of the polynomial at n distinct points; multiplication becomes pointwise.
- File
point_value_representation.py, 107 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/point-value-representation/. The authors are listed in prov/T0CK84.jsonl.
"""The point-value representation of polynomials over a field.
Fix N distinct points u_0, ..., u_(N-1). A polynomial f of degree less than N
is stored as the list of values [f(u_0), ..., f(u_(N-1))]. Addition and
multiplication are pointwise; conversion to coefficients is interpolation.
Coefficient lists are lowest degree first and normalized (no zero on top).
The elements may be Fraction, a class for Z/pZ, or anything with + - * / and
"== 0". Costs are counted in field operations.
"""
from __future__ import annotations
from typing import Any, List, Sequence
Values = List[Any]
def from_coefficients(a: Sequence[Any], pts: Sequence[Any]) -> Values:
"""The values of a at the points, by Horner's rule. Cost: O(N len(a))."""
out = []
for u in pts:
s: Any = 0
for c in reversed(a):
s = s * u + c
out.append(s)
return out
def add(y: Sequence[Any], z: Sequence[Any]) -> Values:
"""The values of f + g. Cost: N additions."""
return [y[i] + z[i] for i in range(len(y))] # @line L4
def sub(y: Sequence[Any], z: Sequence[Any]) -> Values:
"""The values of f - g. Cost: N subtractions."""
return [y[i] - z[i] for i in range(len(y))]
def mul(y: Sequence[Any], z: Sequence[Any]) -> Values:
"""The values of f g. Cost: N multiplications.
The result represents f g only when deg f + deg g < N.
"""
return [y[i] * z[i] for i in range(len(y))] # @line L5
def weights(pts: Sequence[Any]) -> List[Any]:
"""The barycentric weights w_i = 1 / prod_(j != i) (u_i - u_j). Cost: O(N^2)."""
n = len(pts)
w: List[Any] = []
for i in range(n): # @line L1
p: Any = pts[i] - pts[i] + 1 # the 1 of the field, not the int 1
for j in range(n): # @line L2
if j != i:
p = p * (pts[i] - pts[j])
if p == 0:
raise ValueError("the points must be distinct")
w.append(1 / p) # @line L3
return w
def evaluate_at(y: Sequence[Any], pts: Sequence[Any], w: Sequence[Any], t: Any) -> Any:
"""The value f(t) by the barycentric Lagrange formula. Cost: O(N)."""
for k in range(len(pts)): # @line L6
if t - pts[k] == 0:
return y[k]
ell: Any = 1
for u in pts: # @line L7
ell = ell * (t - u)
s: Any = 0
for i in range(len(pts)): # @line L8
s = s + w[i] * y[i] / (t - pts[i])
return ell * s
def to_coefficients(y: Sequence[Any], pts: Sequence[Any]) -> List[Any]:
"""The coefficient list of the unique f with deg f < N and f(u_i) = y_i.
Newton's divided differences, then nested multiplication. Cost: O(N^2).
"""
n = len(pts)
if n == 0:
return []
d = list(y) # @line L9
for k in range(1, n): # @line L10
for i in range(n - 1, k - 1, -1):
d[i] = (d[i] - d[i - 1]) / (pts[i] - pts[i - k])
c: List[Any] = [d[n - 1]] # @line L11
for i in range(n - 2, -1, -1): # @line L12
# c <- c (x - u_i) + d_i
new = [0] * (len(c) + 1)
for j, cj in enumerate(c):
new[j + 1] = new[j + 1] + cj
new[j] = new[j] - cj * pts[i]
new[0] = new[0] + d[i]
c = new
while c and c[-1] == 0: # @line L13
c.pop()
return c
def extend(y: Sequence[Any], pts: Sequence[Any], new_pts: Sequence[Any]) -> Values:
"""The values of the same polynomial at new points. Cost: O(N^2 + N K) for K new points."""
return from_coefficients(to_coefficients(y, pts), new_pts)
test_point_value_representation.py 184 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/point-value-representation/. The authors are listed in prov/T0CK84.jsonl.
import random
from fractions import Fraction
import pytest
from point_value_representation import (
add,
evaluate_at,
extend,
from_coefficients,
mul,
sub,
to_coefficients,
weights,
)
class ModP:
"""An element of Z/pZ, enough for the tests."""
def __init__(self, v: int, p: int) -> None:
self.p = p
self.v = v % p
def _c(self, o):
return o.v if isinstance(o, ModP) else o
def __add__(self, o):
return ModP(self.v + self._c(o), self.p)
__radd__ = __add__
def __sub__(self, o):
return ModP(self.v - self._c(o), self.p)
def __rsub__(self, o):
return ModP(self._c(o) - self.v, self.p)
def __mul__(self, o):
return ModP(self.v * self._c(o), self.p)
__rmul__ = __mul__
def __neg__(self):
return ModP(-self.v, self.p)
def __truediv__(self, o):
return ModP(self.v * pow(self._c(o) % self.p, -1, self.p), self.p)
def __rtruediv__(self, o):
return ModP(self._c(o) * pow(self.v, -1, self.p), self.p)
def __eq__(self, o):
return self.v == self._c(o) % self.p
def __hash__(self):
return hash(self.v)
# Independent schoolbook oracle on plain numbers.
def school_mul(a, b, p=None):
if not a or not b:
return []
c = [0] * (len(a) + len(b) - 1)
for i, x in enumerate(a):
for j, y in enumerate(b):
c[i + j] += x * y
return strip(c, p)
def strip(a, p=None):
a = [x % p for x in a] if p is not None else list(a)
while a and a[-1] == 0:
a.pop()
return a
def raw(xs):
return [x.v if isinstance(x, ModP) else x for x in xs]
def test_example():
pts = [Fraction(0), Fraction(1), Fraction(2)]
f, g = [1, 2], [3, -1] # 1 + 2x and 3 - x
yf, yg = from_coefficients(f, pts), from_coefficients(g, pts)
assert yf == [1, 3, 5] and yg == [3, 2, 1]
yh = mul(yf, yg)
assert yh == [3, 6, 5]
assert to_coefficients(yh, pts) == [3, 5, -2]
w = weights(pts)
assert w == [Fraction(1, 2), -1, Fraction(1, 2)]
assert evaluate_at(yh, pts, w, Fraction(3)) == 3 + 15 - 18
assert evaluate_at(yh, pts, w, Fraction(1)) == 6
def test_too_few_points():
pts = [Fraction(0), Fraction(1)]
y = from_coefficients([0, 1], pts) # x
sq = mul(y, y) # x^2 needs 3 points
assert to_coefficients(sq, pts) == [0, 1] # wrong: this is x, not x^2
pts3 = pts + [Fraction(2)]
assert to_coefficients(mul(extend(y, pts, pts3), extend(y, pts, pts3)), pts3) == [0, 0, 1]
def test_edges():
pts = [Fraction(5)]
assert to_coefficients([Fraction(7)], pts) == [7]
assert to_coefficients([Fraction(0)], pts) == []
assert to_coefficients([], []) == []
pts = [Fraction(x) for x in (0, 1, 2, 3)]
assert to_coefficients([0, 0, 0, 0], pts) == []
assert to_coefficients([4, 4, 4, 4], pts) == [4]
assert from_coefficients([], pts) == [0, 0, 0, 0]
with pytest.raises(ValueError):
weights([Fraction(1), Fraction(1)])
# far apart degrees: a constant times a degree 5 polynomial
big = [Fraction(x) for x in (1, -2, 0, 3, 0, 1)]
pts6 = [Fraction(x) for x in range(-3, 3)]
prod = mul(from_coefficients([Fraction(-4)], pts6), from_coefficients(big, pts6))
assert to_coefficients(prod, pts6) == [-4 * x for x in big]
FIELDS = [
("F2", 2),
("F3", 3),
("F7", 7),
("F998244353", 998244353),
("Q", None),
]
@pytest.mark.parametrize("name,p", FIELDS)
def test_random_against_schoolbook(name, p):
rng = random.Random(sum(map(ord, name)) * 13)
for _ in range(500):
n_pts = rng.randint(1, 9 if p is None or p > 9 else p)
if p is None:
vals = rng.sample(range(-20, 20), n_pts)
pts = [Fraction(v, rng.choice([1, 2, 3])) for v in vals]
if len(set(pts)) < n_pts:
continue
make = lambda x: Fraction(x)
else:
pts = [ModP(v, p) for v in rng.sample(range(min(p, 10**6)), n_pts)]
make = lambda x: ModP(x, p)
da = rng.randint(0, n_pts - 1)
db = rng.randint(0, n_pts - 1 - da)
a_int = strip([rng.randint(-9, 9) for _ in range(da + 1)], p)
b_int = strip([rng.randint(-9, 9) for _ in range(db + 1)], p)
a = [make(x) for x in a_int]
b = [make(x) for x in b_int]
ya, yb = from_coefficients(a, pts), from_coefficients(b, pts)
want = school_mul(a_int, b_int, p)
got = to_coefficients(mul(ya, yb), pts)
assert strip(raw(got), p) == want if p is not None else got == want
s = to_coefficients(add(ya, yb), pts)
n = max(len(a_int), len(b_int))
want_s = strip([(a_int[i] if i < len(a_int) else 0) + (b_int[i] if i < len(b_int) else 0) for i in range(n)], p)
assert strip(raw(s), p) == want_s if p is not None else s == want_s
d = to_coefficients(sub(ya, ya), pts)
assert d == []
# barycentric evaluation at a fresh point agrees with Horner
w = weights(pts)
t = make(rng.randint(-50, 50))
yab = mul(ya, yb)
assert evaluate_at(yab, pts, w, t) == from_coefficients(got, [t])[0]
# round trip
assert to_coefficients(ya, pts) == a or strip(raw(to_coefficients(ya, pts)), p) == a_int
def test_integer_coefficients_through_rationals():
rng = random.Random(77)
for _ in range(300):
a = strip([rng.randint(-50, 50) for _ in range(rng.randint(0, 6))])
b = strip([rng.randint(-50, 50) for _ in range(rng.randint(0, 6))])
n = max(len(a) + len(b) - 1, 1)
pts = [Fraction(i) for i in range(n)]
got = to_coefficients(mul(from_coefficients(a, pts), from_coefficients(b, pts)), pts)
assert got == school_mul(a, b)
What the next level needs
- the tests did not run in this build