"""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())