recip.py
Checks: COLLIDE: The reciprocal pentanomials theorem.
Cited on: COLLIDE.
To run it: save recip.py; needs flint, mpmath, numpy.
"""recip.py -- the reciprocal 0/1 pentanomials read by Mahler measure:
an effective bound M(f) <= 4^(1/a0) M(F) puts every ray but finitely
many out of reach of two non-cyclotomic factors, given no reciprocal
measure in (1, g] at some g > sqrt M(F), and Capelli decides the rest.
QUESTION. f = 1 + x^q + x^p + x^(2p-q) + x^(2p), 0 < q < p, the
reciprocal 0/1 pentanomials. With a = p, b = p - q, k = gcd(a, b) and
(a, b) = k(a0, b0), x^(-p) f = F(x^a, x^b) for F = 1 + u + 1/u + v + 1/v
(Boyd's polynomial), and f = f0(x^k), f0 = 1 + x^(a0-b0) + x^a0 +
x^(a0+b0) + x^(2a0). Is the non-cyclotomic part of f ever reducible?
H(g): no reciprocal non-cyclotomic irreducible polynomial in Z[x] has
Mahler measure in (1, g]. Lehmer's conjecture in its strong form
(Lehmer's number least) gives H(1.176280); for g below Lehmer's number,
H(g) is weaker than that form, and with Smyth stronger than the usual
form (some gap exists). C. J. Smyth (known): a non-reciprocal
non-cyclotomic irreducible has measure at least theta0 = 1.324718.
THE BOUND (by hand, before the engine). On the circle F(e^(i t1),
e^(i t2)) = 1 + 2cos t1 + 2cos t2, so m(f0) = (1/2pi) int log|1 +
2cos(a0 t) + 2cos(b0 t)| dt. Put phi = a0 t: at the a0 preimages of phi
the second angle runs over psi0 + 2pi r/a0, each r once (gcd(a0, b0) =
1), psi0 = b0 phi/a0. With c = 1 + 2cos phi = -(w + 1/w), |c + 2cos psi|
= |e^(i psi) - w||e^(i psi) - 1/w|, and prod_r (X - zeta^r z) = X^a0 -
z^a0 gives, exactly,
m(f0) = (1/2pi) int (1/a0) log|w^a0 + w^(-a0) - 2cos(b0 phi)| dphi.
Where c > 2 (|phi| < pi/3) w is real, |w| = e^s, s = arccosh(c/2), and
the integrand is at most s + (2/a0) log(1 + e^(-a0 s)) <= s + (2/a0) log
2; where -1 <= c <= 2, |w| = 1 and it is at most (1/a0) log 4. Jensen in
v gives m(F) = (1/2pi) int_(|phi| < pi/3) s dphi, so
m(f0) <= m(F) + (log 4)/a0, M(f) = M(f0) <= 4^(1/a0) M(F).
THE REDUCTION, given H(g).
(1) Two non-cyclotomic irreducible factors of f, with multiplicity,
give M(f) > g^2: a non-reciprocal one w brings w~, distinct, and
M(w) >= theta0 > g. So every a0 with 4^(1/a0) M(F) < g^2 is clean
at every k: A(g) is the least such a0, which exists exactly when
g > sqrt M(F) = 1.133902.
(2) For a0 < A(g), f0 is factored. If it has two non-cyclotomic
factors, with multiplicity, its ray fails at k = 1; if none, its
ray is cyclotomic.
If one, w, with a root alpha and K = Q(alpha): the factors of
f0(x^k) off the cyclotomic ones are those of w(x^k), and by
Capelli w(x^k) is reducible iff alpha is in K^l for a prime l | k,
or 4 | k and alpha is in -4K^4. alpha = beta^l, beta in K, forces
Q(beta) = K and M(beta) = M(w)^(1/l) with beta no root of unity,
so M(w)^(1/l) > g (H(g), or theta0 > g when beta's polynomial is
not reciprocal): l < L(a0) = log(4^(1/a0) M(F))/log g. So the
ray is clean iff w(x^4) is irreducible (it holds alpha in K^2 and
alpha in -4K^4) and w(x^l) is, for every odd prime l <= L(a0).
PREDICTIONS, frozen before the engine.
P1 the bound holds at every primitive pair swept.
P2 the true error is far below (log 4)/a0: a0 |m(f0) - m(F)| tends to
0 (carried over from stationary phase on the smooth curve 1 +
2cos t1 + 2cos t2 = 0, which suggests a0^(-3/2); not derived).
P3 the primitives with M(f0) >= 1.176280^2 all have a0 <= 8 (a0 < 19
is forced by the bound).
P4 every w(x^4), w(x^l) of step (2) is irreducible.
P5 g can go toward sqrt(M(F)) ~ 1.134 at the price of a longer step
(2): degree 8 a0 at a0 < A(g).
KILLS (frozen). K1 "bound violations: N" over every primitive pair, N
> 0 beyond 1e-9 kills the bound. K2 "two non-cyclotomic factors: N" in
step (2); N > 0 is a reciprocal NO on that ray. K3 m(F) by quadrature
and m(f0) at the largest a0 differing by more than (log 4)/a0 is an
instrument fault (a loose cross-check of the quadrature: the bound is
one-sided).
CONTROLS. C1 the Chebyshev-form measure on Lehmer's L(x)^2 and on
L(x) S8(x), S8 the degree-8 Salem polynomial; L^2's roots are double,
found to about sqrt(eps), so its bar is 1e-6, L S8's 1e-9. C2 that
measure against numpy's roots of f0 itself, every primitive a0 <= 14. C3
the step-(2) reader on w(x^4) times L(x) at one planted primitive: two
factors. C4 the bound check with log 4 replaced by 0 reports violations.
METHOD. x^(-a0) f0 = h(s), h = 1 + 2T_a0(s) + 2T_b0(s) in Chebyshev
form, s = (x + 1/x)/2; a root s of h carries the root pair
x = s +- sqrt(s^2 - 1), product 1, so m(f0) = sum over the roots of h
of log max |x|. numpy's chebroots (colleague matrix, a0 x a0). m(F) by
mpmath quadrature at 30 digits. Step (2) by flint's exact factoring.
FINDINGS. 9 of 9 checks. C1 L^2 "M = 1.383636622 (known 1.383636563,
|error in log M| 4.3e-08", L S8 "1.506390098, |error in log M|
7.3e-15"; C2 "largest difference 2.99e-14". "m(F) =
0.25133043371325223137, M(F) = 1.2857348642919862749". SWEEP to a0 =
200: "K1 bound violations: 0", the kill misses; C4 fires 7314 times.
Per decade of a0 (2..9, 10..99, 100..200) max a0|e| 1.005, 0.528,
0.209 and max a0^1.5|e| 2.011, 2.100, 2.107: P2 held, the rate
a0^(-3/2), its constant 2.0 to 2.1 unexplained; the bound slack by
"max a0 e/log 4 0.342" at worst. P3 held: one primitive reaches
Lehmer^2, "(a0, b0) = (3, 2): M(f0) = 1.506136". K3 within the bound.
CAPELLI at g = 1.15: "A(g) = 50 (2 log g - bound at A(g):
4.676e-04)", "f0 cyclotomic 3, one non-cyclotomic factor 750, two or
more 0; w(x^l) read at l in [3, 4, 5], degree <= 392", "K2 two
non-cyclotomic factors: 0": P4 held. At --gap 1.176280818, just below
Lehmer's number, "A(g) = 19", 98 with one factor, l in [3, 4], K2 0.
P5: g goes down to sqrt M(F) only as A(g) grows without bound. So given
H(1.15), no reciprocal 0/1 pentanomial has a reducible
non-cyclotomic part; what is computed is m(F)'s value, to a margin of
4.7e-04 in log, and flint's exact factoring; the sweep reads only P2
and P3.
RUN RECORD. The bare run: 9 of 9, 51.0 s, peak 78.0 MB; --rehearse
(the sweep to 40, Capelli to a0 <= 12) 0.1 s; --top 40 --gap
1.176280818 0.2 s. --gap 1.13 stops at the gap control ("g = 1.13,
sqrt M(F) = 1.133902"). C1 L^2 failed at a 1e-9 bar
(4.3e-08): its roots are double. --top 40 --gap 1.14 (A(g) = 130) was
stopped past 10 min, past the a0^4 scaling of the 1.15 run's 11.0 s
(about 8.4 min): its cost is unmeasured.
RESOURCE. Eigenvalues of an a0 x a0 matrix per pair: about 0.3 A^2
pairs. Step (2) factors w(x^4) of degree up to 8 a0 at a0 < A(g).
"""
import os
import sys
import time
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
from math import gcd, log # noqa: E402
import mpmath # noqa: E402
import numpy as np # noqa: E402
from flint import fmpz_poly # noqa: E402
mpmath.mp.dps = 30
CHECKS = []
REHEARSE = "--rehearse" in sys.argv
def arg(name, default, cast):
return cast(sys.argv[sys.argv.index(name) + 1]) if name in sys.argv \
else default
A = 40 if REHEARSE else arg("--top", 200, int)
GAP = arg("--gap", 1.15, float)
LEHMER = 1.176280818259918
L10 = [1, 1, 0, -1, -1, -1, -1, -1, 0, 1, 1]
S8 = [1, 0, 0, -1, -1, -1, 0, 0, 1]
TOL = 1e-9
def check(name, ok, detail=""):
CHECKS.append(bool(ok))
tag = "ok " if ok else "FAIL"
print(f" [{tag}] {name}" + (f" ({detail})" if detail else ""))
def control(name, ok, detail=""):
check(name, ok, detail)
if not ok:
print(" control failed: run stopped")
raise SystemExit(1)
def mF():
"""m(F) = (1/pi) int_0^(pi/3) arccosh((1 + 2cos phi)/2) dphi."""
return mpmath.quad(lambda t: mpmath.acosh((1 + 2 * mpmath.cos(t)) / 2),
[0, mpmath.pi / 3]) / mpmath.pi
def m_cheb(c):
"""log M of the monic reciprocal polynomial whose x^(-d) times is
sum c_j T_j(s), s = (x + 1/x)/2, from the roots s."""
s = np.polynomial.chebyshev.chebroots(c).astype(complex)
z = s + np.sqrt(s * s - 1)
r = np.abs(z)
return float(np.sum(np.log(np.maximum(r, 1 / r))))
def cheb_of_recip(p):
"""Chebyshev coefficients of x^(-d) p(x), p reciprocal of degree 2d."""
d = (len(p) - 1) // 2
return [p[d]] + [2 * p[d + j] for j in range(1, d + 1)]
def h_of(a0, b0):
c = [0] * (a0 + 1)
c[0] += 1
c[a0] += 2
c[b0] += 2
return c
def f0_coeffs(a0, b0):
c = [0] * (2 * a0 + 1)
for e in (0, a0 - b0, a0, a0 + b0, 2 * a0):
c[e] += 1
return c
def m_direct(c):
r = np.abs(np.roots(c[::-1]))
return float(np.sum(np.log(np.maximum(r, 1))))
def polymul(p, q):
return list(np.convolve(p, q))
def primitives(top):
for a0 in range(2, top + 1):
for b0 in range(1, a0):
if gcd(a0, b0) == 1:
yield a0, b0
def section_controls():
print("CONTROLS")
for name, p, want, tol in (("L^2", polymul(L10, L10), 2 * log(LEHMER),
1e-6),
("L S8", polymul(L10, S8),
log(LEHMER) + log(1.280638156267758), 1e-9)):
got = m_cheb(cheb_of_recip(p))
control(f"C1 {name}: M = {np.exp(got):.9f}", abs(got - want) < tol,
f"known {np.exp(want):.9f}, |error in log M| "
f"{abs(got - want):.1e} "
f"against {tol:.0e}")
worst = 0.0
for a0, b0 in primitives(14):
c = f0_coeffs(a0, b0)
worst = max(worst, abs(m_cheb(h_of(a0, b0)) - m_direct(c)))
control("C2 Chebyshev form against numpy's roots of f0, a0 <= 14",
worst < 1e-8, f"largest difference {worst:.2e}")
def section_sweep(mf):
print()
print(f"SWEEP -- every primitive 0 < b0 < a0 <= {A}")
t0 = time.time()
lf = float(mf)
rows = {}
viol = viol0 = 0
hot = []
for a0, b0 in primitives(A):
m = m_cheb(h_of(a0, b0))
e = m - lf
viol += e > log(4) / a0 + TOL
viol0 += e > TOL
if m >= 2 * log(LEHMER):
hot.append((a0, b0, np.exp(m)))
dec = len(str(a0))
r = rows.setdefault(dec, [0, 0.0, 0.0, 0.0, 0.0, -9.0])
r[0] += 1
r[1] = max(r[1], abs(np.exp(m) - np.exp(lf)))
r[2] = max(r[2], a0 * abs(e))
r[3] = max(r[3], a0 ** 1.5 * abs(e))
r[4] = max(r[4], a0 * e / log(4))
r[5] = max(r[5], m)
if a0 == A:
last = m
for dec in sorted(rows):
n, d1, d2, d3, d4, d5 = rows[dec]
lo, hi = 10 ** (dec - 1), min(10 ** dec - 1, A)
print(f" a0 {max(lo, 2):4d}..{hi:4d}: {n:6d} pairs, max |M(f0) - "
f"M(F)| {d1:.3e}, max a0|e| {d2:.3f}, max a0^1.5|e| {d3:.3f}, "
f"max a0 e/log 4 {d4:.3f}, max M(f0) {np.exp(d5):.6f}")
print(f" {time.time() - t0:.1f} s")
check(f"K1 bound violations: {viol}", viol == 0,
f"m(f0) > m(F) + log 4/a0 + {TOL}")
check(f"C4 with log 4 replaced by 0 the check fires: {viol0}",
viol0 > 0)
print(f" M(f0) >= Lehmer^2 = {LEHMER ** 2:.6f}: {len(hot)} primitives,"
f" largest a0 {max(a for a, _, _ in hot)}")
for a0, b0, M in hot:
print(f" (a0, b0) = ({a0}, {b0}): M(f0) = {M:.6f}")
check(f"K3 m(F) and m(f0) at a0 = {A} within the bound",
abs(last - lf) <= log(4) / A,
f"m(F) {lf:.12f}, the a0 = {A} row's last {last:.12f}")
def spread(p, k):
c = [0] * ((len(p) - 1) * k + 1)
for i, v in enumerate(p):
c[i * k] = v
return c
def odd_primes_upto(x):
return [p for p in range(3, int(x) + 1)
if all(p % d for d in range(2, int(p ** 0.5) + 1))]
def noncyc(f):
_, fl = f.factor()
return [g for g, e in fl for _ in range(e)
if g.degree() > 0 and not g.is_cyclotomic()]
def section_capelli(mf, plant=None, limit=None):
print()
lf = float(mf)
lg = log(GAP)
top = next(a for a in range(2, 10 ** 6) if log(4) / a + lf < 2 * lg)
margin = 2 * lg - (log(4) / top + lf)
print(f"CAPELLI -- g = {GAP}: A(g) = {top} (2 log g - bound at A(g): "
f"{margin:.3e}); every primitive a0 < "
f"{'A(g)' if not limit else min(top, limit + 1)}")
t0 = time.time()
count = [0, 0, 0]
bad, deg, ells = [], 0, set()
for a0, b0 in primitives(min(top - 1, limit or top)):
nc = noncyc(fmpz_poly(f0_coeffs(a0, b0)))
if len(nc) != 1:
count[len(nc) and 1] += 1
if len(nc) >= 2:
bad.append((a0, b0, 1))
continue
count[2] += 1
w = [int(c) for c in nc[0].coeffs()]
lmax = (log(4) / a0 + lf) / lg
for ell in [4] + odd_primes_upto(lmax):
ells.add(ell)
f = fmpz_poly(spread(w, ell))
if plant == (a0, b0, ell):
f = f * fmpz_poly(L10)
deg = max(deg, f.degree())
if len(noncyc(f)) >= 2:
bad.append((a0, b0, ell))
print(f" f0 cyclotomic {count[0]}, one non-cyclotomic factor "
f"{count[2]}, two or more {count[1]}; w(x^l) read at l in "
f"{sorted(ells)}, degree <= {deg}; "
f"{time.time() - t0:.1f} s")
return bad
def main():
t0 = time.time()
section_controls()
mf = mF()
print(f" m(F) = {mpmath.nstr(mf, 20)}, M(F) = "
f"{mpmath.nstr(mpmath.exp(mf), 20)}")
section_sweep(mf)
print()
control("the gap is above sqrt M(F) and below Lehmer's number",
2 * log(GAP) > float(mf) and GAP < LEHMER,
f"g = {GAP}, sqrt M(F) = {np.exp(float(mf) / 2):.6f}")
bad = section_capelli(mf, plant=(5, 2, 4), limit=6)
control("C3 a planted L(x) at (5, 2, l = 4) is read", bad == [(5, 2, 4)],
f"{bad}")
bad = section_capelli(mf, limit=12 if REHEARSE else None)
check(f"K2 two non-cyclotomic factors: {len(bad)}", not bad, f"{bad}")
print()
print(f"{sum(CHECKS)} of {len(CHECKS)} checks, {time.time() - t0:.1f} s")
return 0 if all(CHECKS) else 1
if __name__ == "__main__":
sys.exit(main())