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