recip.py
Checks: COLLIDE: The reciprocal pentanomials theorem; theorem; rule; rule; theorem; pattern.
Cited on: COLLIDE.
To run it: save recip.py; needs flint, mpmath, numpy, sympy.
"""recip.py -- the reciprocal 0/1 pentanomials: every ray with a0 <= 90
is clean at every dilation with no hypothesis, by a measure floor and a
power-residue certificate; past it an effective bound M(f) <= 4^(1/a0)
M(F) puts every ray out of reach of two non-cyclotomic factors, given no
reciprocal measure in (1, g] at some g > sqrt(4^(1/91) M(F)).
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). Call f clean when its non-cyclotomic part is 1 or
irreducible. Is every f clean?
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). 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 ray with 4^(1/a0) M(F) < g^2
is clean at every k: A(g) is the least such a0. A(g) <= 91 exactly when
g > G = sqrt(4^(1/91) M(F)) = 1.142572408; so given H(g) at a g above
G, with THE DILATIONS below, every reciprocal f is clean.
THE DILATIONS, a0 <= 90, no hypothesis (by hand, before the engine).
(1) f0 is factored exactly. Two non-cyclotomic factors, with
multiplicity: the ray fails at k = 1. None: f0(x^k) is cyclotomic
at every k. One, h: f0 reciprocal makes h~ = +-h a factor too, so
h is reciprocal, of even degree d >= 2 (degree 1 is x +- 1), and
h(0) = +-1.
(2) Let alpha be a root of h, K = Q(alpha). A root of h(x^k) is no
root of unity, so the non-cyclotomic part of f0(x^k) is h(x^k),
irreducible over Q iff x^k - alpha is irreducible over K (degree
count, [K(alpha^(1/k)) : Q] = dk). Capelli (known; S. Lang,
Algebra, VI.9.1): x^k - alpha is reducible over K iff alpha is in
K^l for a prime l | k, or 4 | k and alpha is in -4K^4. So the ray is
clean at every k iff alpha is in no K^l, l prime, and not in -4K^4.
(3) THE FLOOR. If alpha = beta^l, beta in K, then beta is an algebraic
integer (a root of x^l - alpha), no root of unity, and Q(beta) = K
(alpha in Q(beta)), so beta has degree d and its d conjugates are
the sigma(beta), whose l-th powers are alpha's: M(alpha) =
M(beta)^l. P. Voutier, An effective lower bound for the height of
algebraic numbers, Acta Arith. 74 (1996) 81-95, Corollary 2, as
printed: a non-zero algebraic number of degree d >= 2, not a root
of unity, has log M > 2/(log(3d))^3 (derived there from his theorem
with Boyd's computations and Smyth's theorem). And M(alpha) = M(h)
<= M(f0) <= ||f0||_2 = sqrt 5: Jensen's inequality for the concave
log gives m(f0) = (1/2) int log|f0|^2 <= (1/2) log int |f0|^2, which
is log ||f0||_2 by Parseval. So
l < L(d) = (log 5)(log 3d)^3 / 4.
By hand: d = 180, (log 540)^3 = 249.04, L = 100.2, so l <= 97;
d = 4, L = 6.17; d = 2, L = 2.31, so L > 2 at every d and l = 2 is
always read. No hypothesis, and no computed constant.
(4) THE CERTIFICATE, at each prime l < L. Let p be a prime, p = 1 mod
l (p odd at l = 2), and r a simple root of h mod p. By Hensel r
lifts to a root of h in Z_p, an embedding K -> Q_p with alpha ->
that root. alpha is a unit, so r != 0. If alpha = beta^l, beta is
integral, so its image lies in Z_p and reduces to b with b^l = r in
F_p: r^((p-1)/l) = 1. So one (p, r) with r^((p-1)/l) != 1 proves
alpha not in K^l. At -4K^4, p odd: alpha = -4 gamma^4 puts gamma's
image in Z_p^x (|gamma^4|_p = |alpha/4|_p = 1), so -r/4 is a fourth
power in F_p; one (p, r) with c^((p-1)/gcd(4, p-1)) != 1, c = -r/4
mod p, proves alpha not in -4K^4. This replaces factoring h(x^l),
degree up to 97 * 180.
(5) A pair (h, l), or (h, -4), with no certificate among the first 50
primes p that qualify (p = 1 mod l, h with a simple root mod p;
every simple root tried) is factored outright by flint: h(x^l), or
h(x^4) for -4; two non-cyclotomic factors are a NO on the ray, and
one is counted undecided, conservatively.
So "verified a0 <= 90, every k" rests on flint's exact factoring of
the primitives and its roots mod p, Capelli and Voutier's Corollary 2;
m(F) enters only
THE REDUCTION, as G.
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 primitive a0 <= 90 has at most one non-cyclotomic factor (an
earlier run, its floor from THE BOUND, l < 33.3, found none with two).
P5 every (h, l) and (h, -4) is certified; no fallback. Landau's floor
adds the primes 37 to 97 over the bound's floor, each a fresh chance
of about 1 - 1/l per root.
P6 the largest L is 100.2 (d = 180), l up to 97; the pairs about 2.5
times the bound floor's 25733.
P7 the dilations take under 6 minutes (the bound floor's 131 s times
about 2.5).
KILLS (frozen). K1 "bound violations: N" over every primitive pair, N
> 0 beyond 1e-9 kills the bound. K2 "NO triples (a0, b0, l): N" over
every primitive a0 <= 90, l = 1 the factor step; N > 0 is a reciprocal
counterexample (a find). "undecided pairs (h, l): N" > 0 leaves the rule
unproved at those rays (instrument). 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 dilation sweep's own path at a0 <= 5 with one primitive, (5, 2),
read in place of f0 as: f0 L(x), a NO at l = 1 (the factor step); h3,
the minimal polynomial of lambda^3, lambda a root of L(x), a NO at l =
3 (floor, certificate, fallback); h2, of lambda^2, a NO at l = 2; h4,
of -4 lambda^4, a NO at -4; each the only triple. Then h3 read with the
floor cut to 2.5 misses its NO, by construction: it reads only ells (a
plant at the floor's top would need M(h) above sqrt 5). C4 the
bound check with log 4 replaced by 0 reports violations. C5 L(x) is
certified at l = 3 (the certificate fires). C6 a cross-check: certified
pairs whose h(x^l) splits, counted; it counts the planted (h3, 3) under
a certifier that passes everything, then 0 over every pair at a0 <= 10.
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. The dilations: flint's exact factoring
of f0 and of h(x^l) at a fallback, flint's roots of h mod p for the
certificates, sympy's resultant for the planted powers of lambda.
FINDINGS. SWEEP to a0 = 200: 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". "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.
GAP "G = sqrt(4^(1/91) M(F)) = 1.14257240829; at g = 1.142573: A(g) =
91 (2 log g - bound at A(g): 1.036e-06)", m(F) to 30 digits.
DILATIONS. C3 each plant read as its one NO ("[(5, 2, 1)]", "[(5, 2,
3)]", "[(5, 2, 2)]", "[(5, 2, -4)]"), the cut floor missing h3's;
C5 "p = 67, 1 tried"; C6 "1 planted split counted", then "0 of 230
pairs". To a0 = 90: "f0 cyclotomic 3, one non-cyclotomic factor 2476,
two or more 0; h of degree <= 180, largest L 100.21, l up to 97";
"pairs (h, l) and (h, -4): 55020, certified 55020 (53098 at the first
p); most p tried for a certificate 8"; "K2 NO triples (a0, b0, l): 0",
"K2 undecided pairs (h, l): 0". The kill misses: every reciprocal 0/1
pentanomial with a0 <= 90 is clean at every k, with no hypothesis
(flint, Capelli, Voutier's Corollary 2). P4, P5 and P7 held, and P6's
L and l; its pairs, 2.1 times 25733, not 2.5. No float enters the
verdict but L, read with a 1e-9 widening, far above double precision;
the printed margin to the next prime above L is the harmless side (an
extra prime only adds work).
RUN RECORD. The bare run 17 of 17, 333.6 s (the sweep 40.1 s, the
dilations 292.6 s), peak 92.4 MB; --rehearse 3.1 s. C1 L^2 failed at a
1e-9 bar (4.3e-08): its roots are double. The first bare run, under a
clock of three times the script's former 51.0 s, was killed at 154 s.
RESOURCE. Eigenvalues of an a0 x a0 matrix per pair: about 0.3 A^2
pairs. The dilations factor f0 of degree <= 180 and find roots mod p of
h of degree <= 180, 55020 pairs.
"""
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, nmod_poly # noqa: E402
import sympy # 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)
TOP = 30 if REHEARSE else 90
GAP = arg("--gap", 1.142573, 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
TRIES = 50
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 [int(v) for v in 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 section_gap(mf):
"""A(g), and the least g with A(g) <= 91."""
print()
lf = float(mf)
lg = log(GAP)
a = next(a for a in range(2, 10 ** 6) if log(4) / a + lf < 2 * lg)
g91 = mpmath.sqrt(mpmath.exp(mf + mpmath.log(4) / 91))
print(f"GAP -- G = sqrt(4^(1/91) M(F)) = {mpmath.nstr(g91, 12)}; at g = "
f"{GAP}: A(g) = {a} (2 log g - bound at A(g): "
f"{2 * lg - (log(4) / a + lf):.3e})")
control("the gap is above G and below Lehmer's number",
a <= 91 and GAP < LEHMER, f"g = {GAP}, A(g) = {a}")
def is_prime(n):
return n > 1 and all(n % d for d in range(2, int(n ** 0.5) + 1))
def ells(L):
"""Every prime l < L, and -4."""
return [n for n in range(2, int(L) + 1) if is_prime(n) and n < L] + [-4]
def floor_L(d):
"""L(d): alpha = beta^l forces l < log M(h) (log 3d)^3 / 2 (Voutier,
Cor. 2), with M(h) <= ||f0||_2 = sqrt 5."""
return log(5) * log(3 * d) ** 3 / 4
def spread(p, k):
c = [0] * ((len(p) - 1) * k + 1)
for i, v in enumerate(p):
c[i * k] = v
return c
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 splits(h, ell):
"""h(x^l), or h(x^4) at -4, has two non-cyclotomic factors."""
return len(noncyc(fmpz_poly(spread(h, 4 if ell == -4 else ell)))) >= 2
def power_of_roots(h, e):
"""The monic polynomial whose roots are c * (the roots of h)^e, as
fmpz coefficients, by the resultant in y of h(y) and x - c y^e."""
c, e = e if isinstance(e, tuple) else (1, e)
d = len(h) - 1
y, x = sympy.symbols("y x")
hy = sympy.Poly(list(reversed(h)), y)
res = sympy.resultant(hy.as_expr(), x - c * y ** e, y)
p = sympy.Poly(res, x)
co = [int(v) for v in reversed(p.all_coeffs())]
lead = co[-1]
assert all(v % lead == 0 for v in co) and len(co) == d + 1
return [v // lead for v in co]
def certify(h, ell):
"""(certified?, qualifying p tried, p used). ell a prime, or -4."""
step = 2 if ell in (2, -4) else 2 * ell # odd p = 1 mod ell
tried, p = 0, 1
while tried < TRIES:
p += step
if not is_prime(p):
continue
roots = nmod_poly([c % p for c in h], p).roots()
simple = [int(r) for r, m in roots if m == 1 and int(r) != 0]
if not simple:
continue
tried += 1
for r in simple:
if ell == -4:
c = (-r * pow(4, -1, p)) % p
if pow(c, (p - 1) // gcd(4, p - 1), p) != 1:
return True, tried, p
elif pow(r, (p - 1) // ell, p) != 1:
return True, tried, p
return False, tried, None
def read_h(h, stats, L=None):
"""Every prime l < L (by default L(deg h)) and -4: the list of
(l, certified, NO), a NO when the fallback factors a split."""
if L is None:
L = floor_L(len(h) - 1)
q = next(n for n in range(int(L) + 1, 10 ** 6) if is_prime(n))
stats["L"] = max(stats.get("L", 0), L)
stats["margin"] = min(stats.get("margin", 9.0), (q - L) / L)
out = []
for ell in ells(L * (1 + 1e-9)):
ok, tried, _ = certify(h, ell)
no = not ok and splits(h, ell)
stats["pairs"] = stats.get("pairs", 0) + 1
stats["cert"] = stats.get("cert", 0) + ok
stats["first"] = stats.get("first", 0) + (ok and tried == 1)
stats["tried"] = max(stats.get("tried", 0), tried if ok else 0)
stats["lmax"] = max(stats.get("lmax", 0), ell)
out.append((ell, ok, no))
return out
def dilate(top, plant=None):
"""Every primitive a0 <= top: (counts, stats, NO triples, undecided).
plant = (a0, b0, poly) reads poly in place of f0 at that primitive."""
count, stats, nos, undecided = [0, 0, 0], {"d": 0}, [], []
for a0, b0 in primitives(top):
f = fmpz_poly(f0_coeffs(a0, b0))
if plant and plant[:2] == (a0, b0):
f = fmpz_poly(plant[2])
nc = noncyc(f)
count[min(len(nc), 2)] += 1
if len(nc) >= 2:
nos.append((a0, b0, 1))
if len(nc) != 1:
continue
h = [int(c) for c in nc[0].coeffs()]
stats["d"] = max(stats["d"], len(h) - 1)
for ell, ok, no in read_h(h, stats):
if no:
nos.append((a0, b0, ell))
elif not ok:
undecided.append((a0, b0, ell))
return count, stats, nos, undecided
def crosscheck(pairs, cert):
"""How many pairs cert passes whose h(x^l) splits."""
return sum(bool(cert(h, ell)[0]) and splits(h, ell) for h, ell in pairs)
def section_dilate_controls():
print()
print("DILATION CONTROLS")
h3 = power_of_roots(L10, 3)
for name, poly, ell in (("f0 L(x)", polymul(f0_coeffs(5, 2), L10), 1),
("h3", h3, 3),
("h2", power_of_roots(L10, 2), 2),
("h4", power_of_roots(L10, (-4, 4)), -4)):
_, _, nos, und = dilate(5, plant=(5, 2, poly))
control(f"C3 {name} at (5, 2): the one NO at l = {ell}",
nos == [(5, 2, ell)] and not und, f"{nos}, undecided {und}")
out = read_h(h3, {}, L=2.5)
control("C3 h3 with the floor cut to 2.5: its NO missed",
not any(no for _, _, no in out), f"{out}")
ok, tried, p = certify(L10, 3)
control("C5 L(x) certified at l = 3", ok, f"p = {p}, {tried} tried")
planted = crosscheck([(h3, 3)], lambda h, ell: (True,))
control(f"C6 a certifier passing everything: {planted} planted split "
f"counted", planted == 1)
pairs = []
for a0, b0 in primitives(10):
nc = noncyc(fmpz_poly(f0_coeffs(a0, b0)))
if len(nc) == 1:
h = [int(c) for c in nc[0].coeffs()]
pairs += [(h, ell) for ell in ells(floor_L(len(h) - 1))]
bad = crosscheck(pairs, certify)
control(f"C6 certified and split at a0 <= 10: {bad} of {len(pairs)} "
f"pairs", bad == 0 and pairs)
def section_dilate():
print()
print(f"DILATIONS -- every primitive 0 < b0 < a0 <= {TOP}, every k")
t0 = time.time()
count, s, nos, und = dilate(TOP)
print(f" f0 cyclotomic {count[0]}, one non-cyclotomic factor "
f"{count[1]}, two or more {count[2]}; h of degree <= {s['d']}, "
f"largest L {s['L']:.2f}, l up to {s['lmax']}; least (next prime "
f"- L)/L {s['margin']:.2e}")
print(f" pairs (h, l) and (h, -4): {s['pairs']}, certified "
f"{s['cert']} ({s['first']} at the first p); most p tried for a "
f"certificate {s['tried']}; {time.time() - t0:.1f} s")
check(f"K2 NO triples (a0, b0, l): {len(nos)}", not nos, f"{nos[:10]}")
check(f"K2 undecided pairs (h, l): {len(und)}", not und, f"{und[:10]}")
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)
section_gap(mf)
section_dilate_controls()
section_dilate()
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())