Vincent-Collins-Akritas algorithm in Python
Isolates real roots by repeated Descartes-sign counts and Mobius transformations of the interval.
- File
vincent_collins_akritas.py, 174 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/vincent-collins-akritas/. The authors are listed in prov/S6W6AA.jsonl.
"""The Vincent-Collins-Akritas (Descartes) method for real root isolation.
Polynomials are lists of coefficients, lowest degree first; the zero
polynomial is []. The method works with ints only: Taylor shifts by 1,
multiplications by powers of two, and sign tests.
"""
from __future__ import annotations
from fractions import Fraction
from math import gcd
from typing import Sequence
Interval = tuple[Fraction, Fraction]
def _strip(a: list) -> list:
while a and a[-1] == 0:
a.pop()
return a
def _taylor_shift(a: list[int], c: int) -> list[int]:
"""Return the coefficients of a(x + c), by the classical O(n^2) scheme."""
a = list(a)
n = len(a) - 1
for i in range(n):
for j in range(n - 1, i - 1, -1):
a[j] += c * a[j + 1]
return a
def _variations(a: Sequence[int]) -> int:
"""Number of sign changes in a, zeros ignored."""
v, last = 0, 0
for x in a:
if x:
if last and (x > 0) != (last > 0):
v += 1
last = x
return v
def _descartes_test(a: list[int]) -> int:
"""Var of (x + 1)^n a(1 / (x + 1)): bounds the roots of a in (0, 1)."""
return _variations(_taylor_shift(a[::-1], 1))
def _poly_rem(a: list[Fraction], b: list[Fraction]) -> list[Fraction]:
a = list(a)
while len(a) >= len(b):
q = a[-1] / b[-1]
s = len(a) - len(b)
for i, y in enumerate(b):
a[s + i] -= q * y
_strip(a)
return a
def _poly_gcd(a: list[Fraction], b: list[Fraction]) -> list[Fraction]:
while b:
a, b = b, _poly_rem(a, b)
return a
def _poly_divexact(a: list[Fraction], b: list[Fraction]) -> list[Fraction]:
a = list(a)
q = [Fraction(0)] * (len(a) - len(b) + 1)
while len(a) >= len(b):
t = a[-1] / b[-1]
s = len(a) - len(b)
q[s] = t
for i, y in enumerate(b):
a[s + i] -= t * y
_strip(a)
return q
def _to_primitive_int(a: Sequence[Fraction]) -> list[int]:
den = 1
for x in a:
den = den * x.denominator // gcd(den, x.denominator)
b = [int(x * den) for x in a]
c = 0
for x in b:
c = gcd(c, x)
if b[-1] < 0:
c = -c
return [x // c for x in b]
def _derivative(a: Sequence) -> list:
return [i * a[i] for i in range(1, len(a))]
def square_free_part(f: Sequence[int | Fraction]) -> list[int]:
"""Primitive integer polynomial with the same distinct roots as f.
Precondition: f is not the zero polynomial.
Cost: O(n^2) operations on rationals, by Euclid's algorithm.
"""
a = _strip([Fraction(x) for x in f])
if not a:
raise ValueError("the zero polynomial has no isolating intervals")
if len(a) == 1:
return [1]
g = _poly_gcd(a, _strip(_derivative(a)))
return _to_primitive_int(_poly_divexact(a, g))
def vca_isolate(f: Sequence[int]) -> list[Interval]:
"""Isolating intervals for the real roots of a square-free f in Z[x].
Precondition: f has int coefficients, is not zero, and is square-free
(no repeated complex root); leading zeros in the list are allowed.
Returns a sorted list of pairs (lo, hi) of Fractions. A pair with
lo < hi is an open interval holding exactly one root; a pair with
lo == hi is an exact rational root. Distinct pairs are disjoint.
Cost: O(n^3 (L + log n)) additions, shifts and sign tests, where L is
the coefficient length, from the recursion tree bound of Eigenwillig,
Sharma and Yap and O(n^2) work per node.
"""
f = _strip([int(x) for x in f])
if not f:
raise ValueError("the zero polynomial has no isolating intervals")
n = len(f) - 1
if n == 0:
return []
if square_free_part(f) != _to_primitive_int([Fraction(x) for x in f]):
raise ValueError("f must be square-free; use isolate_real_roots")
bound = 1 + max(Fraction(abs(x), abs(f[n])) for x in f[:n]) # @line L1
big = 1
while big < bound:
big *= 2
scaled = [x * big**i for i, x in enumerate(f)] # @line L2
scaled = _taylor_shift(scaled, -1)
a0 = [x << i for i, x in enumerate(scaled)]
stack: list[tuple[list[int], int, int]] = [(a0, 0, 0)] # @line L3
out: list[Interval] = []
width = 2 * big
def ends(c: int, k: int) -> Interval:
lo = Fraction(-big) + Fraction(c * width, 2**k)
return lo, lo + Fraction(width, 2**k)
while stack: # @line L4
a, c, k = stack.pop() # @line L5
v = _descartes_test(a) # @line L6
if v == 0: # @line L7
continue
if v == 1: # @line L8
out.append(ends(c, k))
continue
left = [x << (n - i) for i, x in enumerate(a)] # @line L9
right = _taylor_shift(left, 1)
if right[0] == 0: # @line L10
mid = ends(c, k)[0] + Fraction(width, 2 ** (k + 1))
out.append((mid, mid))
stack.append((right, 2 * c + 1, k + 1)) # @line L11
stack.append((left, 2 * c, k + 1))
out.sort() # @line L12
return out
def isolate_real_roots(f: Sequence[int | Fraction]) -> list[Interval]:
"""Isolating intervals for the distinct real roots of any nonzero f.
Precondition: f is nonzero, with int or Fraction coefficients.
Reduces f to its square-free part, then calls vca_isolate.
Cost: O(n^2) rational operations for the reduction, plus vca_isolate.
"""
return vca_isolate(square_free_part(f))
test_vincent_collins_akritas.py 344 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/vincent-collins-akritas/. The authors are listed in prov/S6W6AA.jsonl.
"""Tests for the Vincent-Collins-Akritas method.
Real root isolation needs an ordered field, so the inputs are over Z and
over Q (Fraction); prime fields Z/pZ have no order and are not used. The
oracle is an independent Sturm sequence count in exact rational arithmetic.
"""
from __future__ import annotations
import random
from fractions import Fraction
from math import comb
import pytest
from vincent_collins_akritas import (
_descartes_test,
_taylor_shift,
isolate_real_roots,
square_free_part,
vca_isolate,
)
# ------------------------------------------------------------ the oracle
def _trim(a):
a = list(a)
while a and a[-1] == 0:
a.pop()
return a
def _rem(a, b):
a = [Fraction(x) for x in a]
while len(a) >= len(b) and a:
q = a[-1] / b[-1]
s = len(a) - len(b)
for i, y in enumerate(b):
a[s + i] -= q * y
a = _trim(a)
return a
def _sturm(f):
seq = [_trim([Fraction(x) for x in f])]
d = _trim([i * seq[0][i] for i in range(1, len(seq[0]))])
if not d:
return seq
seq.append(d)
while True:
r = _rem(seq[-2], seq[-1])
if not r:
return seq
seq.append([-x for x in r])
def _ev(a, x):
s = Fraction(0)
for c in reversed(a):
s = s * x + c
return s
def _var(vals):
v = [x for x in vals if x != 0]
return sum(1 for i in range(len(v) - 1) if (v[i] > 0) != (v[i + 1] > 0))
def _count(seq, lo, hi):
"""Distinct roots in (lo, hi] by Sturm; lo, hi not multiple roots."""
return _var([_ev(p, lo) for p in seq]) - _var([_ev(p, hi) for p in seq])
def _total(seq):
neg = _var([p[-1] * (-1) ** (len(p) - 1) for p in seq])
pos = _var([p[-1] for p in seq])
return neg - pos
def _sqfree(f):
"""f / gcd(f, f') over Q, computed with the oracle's own remainders."""
a = _trim([Fraction(x) for x in f])
b = _trim([i * a[i] for i in range(1, len(a))])
g = a
while b:
g, b = b, _rem(g, b)
q = [Fraction(0)] * (len(a) - len(g) + 1)
r = list(a)
while len(r) >= len(g):
c = r[-1] / g[-1]
s = len(r) - len(g)
q[s] = c
for i, y in enumerate(g):
r[s + i] -= c * y
r = _trim(r)
assert not r
return q
def _check(f, out):
"""The output isolates exactly the distinct real roots of f.
Endpoints of open intervals may be roots that are reported as exact
points; the Sturm count on (lo, hi] is corrected for that.
"""
h = _sqfree(f)
seq = _sturm(h)
assert len(out) == _total(seq)
for lo, hi in out:
if lo == hi:
assert _ev(h, lo) == 0
else:
assert lo < hi
inside = _count(seq, lo, hi) - (1 if _ev(h, hi) == 0 else 0)
assert inside == 1
for (lo1, hi1), (lo2, hi2) in zip(out, out[1:]):
assert hi1 <= lo2
if lo1 == hi1 and lo2 == hi2:
assert lo1 < lo2
for lo, hi in out:
if lo < hi:
for r in (lo, hi):
if _ev(h, r) == 0:
assert (r, r) in out
def _from_roots(roots):
"""Primitive integer polynomial with the given rational roots."""
f = [Fraction(1)]
for r in roots:
g = [Fraction(0)] * (len(f) + 1)
for i, x in enumerate(f):
g[i + 1] += x
g[i] -= r * x
f = g
den = 1
for x in f:
den = den * x.denominator
return [int(x * den) for x in f]
def _contains(out, r):
return sum(1 for lo, hi in out if (lo == hi == r) or lo < r < hi)
# ------------------------------------------------------------ helpers
def test_taylor_shift_against_binomial_expansion():
rng = random.Random(11)
for _ in range(300):
n = rng.randint(0, 9)
a = [rng.randint(-20, 20) for _ in range(n + 1)]
c = rng.choice([-3, -1, 1, 2, 5])
ref = [0] * (n + 1)
for i, x in enumerate(a):
for j in range(i + 1):
ref[j] += x * comb(i, j) * c ** (i - j)
assert _taylor_shift(a, c) == ref
def test_descartes_test_is_an_upper_bound_with_parity():
rng = random.Random(12)
for _ in range(300):
roots = [Fraction(rng.randint(-9, 9), rng.randint(1, 9)) for _ in range(rng.randint(1, 5))]
roots = sorted(set(roots))
a = _from_roots(roots)
inside = sum(1 for r in roots if 0 < r < 1)
v = _descartes_test(a)
assert v >= inside and (v - inside) % 2 == 0
# ------------------------------------------------------------ worked examples
def test_worked_example_cubic():
out = vca_isolate([1, -3, 0, 1])
assert out == [
(Fraction(-4), Fraction(0)),
(Fraction(0), Fraction(1)),
(Fraction(1), Fraction(2)),
]
def test_lagrange_cubic_close_roots():
out = vca_isolate([7, -7, 0, 1])
assert out == [
(Fraction(-8), Fraction(0)),
(Fraction(1), Fraction(3, 2)),
(Fraction(3, 2), Fraction(2)),
]
def test_exact_midpoint_roots():
out = vca_isolate([0, -1, 0, 1])
assert (Fraction(0), Fraction(0)) in out
_check([0, -1, 0, 1], out)
# ------------------------------------------------------------ edge cases
def test_constants_and_zero():
assert vca_isolate([5]) == []
assert vca_isolate([-3, 0, 0]) == []
with pytest.raises(ValueError):
vca_isolate([])
with pytest.raises(ValueError):
vca_isolate([0, 0])
with pytest.raises(ValueError):
isolate_real_roots([])
def test_degree_one_monic_and_non_monic():
for f in ([3, 2], [-1, 1], [0, 1], [0, -7], [5, 3], [-6, 4, 0, 0]):
out = vca_isolate(f)
_check(f, out)
assert len(out) == 1
def test_no_real_roots():
assert vca_isolate([1, 0, 1]) == []
assert vca_isolate([1, 0, 0, 0, 1]) == []
assert vca_isolate([2, 1, 1]) == []
def test_non_square_free_rejected_and_reduced():
f = [0, 0, -1, 1] # x^2 (x - 1)
with pytest.raises(ValueError):
vca_isolate(f)
out = isolate_real_roots(f)
_check(f, out)
assert len(out) == 2
g = _from_roots([Fraction(1, 2)] * 3 + [Fraction(-2)] * 2 + [Fraction(5)])
out = isolate_real_roots(g)
_check(g, out)
assert len(out) == 3
def test_fraction_coefficients():
f = [Fraction(-1, 3), Fraction(0), Fraction(3, 4)] # roots +- 2/3
out = isolate_real_roots(f)
assert len(out) == 2
assert _contains(out, Fraction(2, 3)) == 1
assert _contains(out, Fraction(-2, 3)) == 1
def test_square_free_part():
assert square_free_part([0, 0, -1, 1]) == [0, -1, 1]
assert square_free_part([4, 4, 1]) == [2, 1]
assert square_free_part([7]) == [1]
# ------------------------------------------------------------ hard inputs
def test_wilkinson_roots_one_to_twelve():
roots = [Fraction(i) for i in range(1, 13)]
f = _from_roots(roots)
out = vca_isolate(f)
_check(f, out)
for r in roots:
assert _contains(out, r) == 1
def test_mignotte_polynomial_close_roots():
# x^7 - 2 (10 x - 1)^2 has two roots very close to 1/10.
f = [0] * 8
f[7] = 1
for i, c in enumerate([1, -20, 100]):
f[i] -= 2 * c
out = vca_isolate(f)
_check(f, out)
assert out == [
(Fraction(819, 8192), Fraction(3277, 32768)),
(Fraction(3277, 32768), Fraction(1639, 16384)),
(Fraction(2), Fraction(4)),
]
def test_chebyshev_t12_all_roots_real():
t = [[1], [0, 1]]
for _ in range(11):
a, b = t[-2], t[-1]
c = [0] + [2 * x for x in b]
for i, x in enumerate(a):
c[i] -= x
t.append(c)
f = t[12]
out = vca_isolate(f)
_check(f, out)
assert len(out) == 12
def test_very_different_degrees_in_one_input():
f = [-2] + [0] * 29 + [1] # x^30 - 2
out = vca_isolate(f)
_check(f, out)
assert len(out) == 2
def test_rational_roots_with_large_denominators():
roots = [Fraction(1, 97), Fraction(2, 97), Fraction(-50, 3), Fraction(13, 7)]
f = _from_roots(roots)
out = vca_isolate(f)
_check(f, out)
for r in roots:
assert _contains(out, r) == 1
# ------------------------------------------------------------ random oracle
def test_random_small_polynomials_against_sturm():
rng = random.Random(2024)
checked = 0
for _ in range(700):
n = rng.randint(1, 6)
f = [rng.randint(-6, 6) for _ in range(n)] + [rng.choice([-3, -2, -1, 1, 2, 3])]
out = isolate_real_roots(f)
_check(f, out)
checked += 1
assert checked >= 500
def test_random_products_of_known_roots():
rng = random.Random(7)
for _ in range(150):
roots = sorted({Fraction(rng.randint(-30, 30), rng.randint(1, 6)) for _ in range(rng.randint(1, 6))})
extra = rng.choice([[], [1, 0, 1], [-2, 0, 1], [3, 1, 1]])
f = _from_roots(roots)
if extra:
g = [0] * (len(f) + len(extra) - 1)
for i, x in enumerate(f):
for j, y in enumerate(extra):
g[i + j] += x * y
f = g
out = isolate_real_roots(f)
_check(f, out)
for r in roots:
assert _contains(out, r) == 1
What the next level needs
- the tests did not run in this build