"""circlemean.py -- does a bound with both sides in closed form keep the band-end witness F_e beating the champion out to depth c r^(3/2)? QUESTION. After x -> -x the champion at rank r, depth J is P = (1 + x)^J (1 - x)^n, n = r - 1, of height H_ch, and F_e = (1 + x)^J (1 - x)^s (x + e (1 - x)^2), s = r - 3, e an integer, lies in the lattice (bandfamily.py). Bound height(F_e) above by its mean modulus on the unit circle and H_ch below by one explicit test of P. J_B(r) is the last depth at which the bound shows some F_e below H_ch: the champion misses h there, so J_ch(r) > J_B(r). Does J_B grow as r^(3/2)? BY HAND. On x = e^(it), 1 + x = e^(it/2) 2cos(t/2), 1 - x = e^(it/2) (-2i) sin(t/2), so (1 - x)^2/x = -4 sin^2(t/2). (i) |F_e| = 2^(J+s) |cos|^J |sin|^s |1 - 4e sin^2| at t/2, and with u = sin^2(t/2), mean(F_e) = (2^(J+s)/pi) B(a, b) E|1 - 4eU|, U ~ Beta(a, b), a = (s + 1)/2, b = (J + 1)/2; the mean deviation splits at u0 = 1/(4e) into regularised incomplete Betas I_u0: E|1 - cU| = 1 - c mu - 2 (1 - I_u0(a, b)) + 2 c mu (1 - I_u0(a+1, b)), mu = a/(a + b), for c = 4e > 1. (ii) For w = x^k (1 + x)^p (1 - x)^q, = (1/2pi) int P conj(W) dt; at 2k = J + n - p - q the phase cancels and the integrand is 2^(J+n+p+q) (-i)^n i^q cos^(J+p) sin^(n+q), nonzero iff n + q and J + p are even: || = (2^(J+n+p+q)/pi) B((n+q+1)/2, (J+p+1)/2), and H_ch >= ||/|w|_1, |w|_1 = 2^p at q = 0 and 2 C(p, p/2 rounded down) at q = 1. L(J) is the best over p of J's parity, q = n mod 2; at q = 0 it falls with p, so p = 0 or 1 (at n, J even the central coefficient, exactly). (iii) n even, real e with 4e = 1/E U: R = mean/L = E|U - EU|/(4 (EU)^2) ~ J/(2 sqrt(pi) s^(3/2)), so J_B ~ 2 sqrt(pi) (r - 2)^(3/2) = 3.545 (r - 2)^(3/2), 3437 at r = 101 (TRANSPLANT: Beta to Gaussian). [After the run: a = (r - 2)/2 makes it s + 1 = r - 2 throughout, so 3492 at r = 101; the frozen 3437 used s = 98.] (iv) n odd: the best p ~ J/n smooths P at its oscillation scale; L/H_ch ~ sqrt(pi/2) e^(-1/2) = 0.76 (TRANSPLANT), so J_B at even r runs about 0.76 of its odd neighbours'. PREDICTIONS (frozen before the engine). P1 J_B/(r - 2)^(3/2) at odd r rises toward 3.545, in [3.0, 3.6] at r = 101. P2 at even r J_B is 0.6 to 1.0 of the mean of its odd neighbours'. P3 J_B(100)/J_F(100) about 0.6. P4 L/H_ch at J_B >= 0.97 at odd r, 0.6 to 0.9 at even r. KILLS. K1 the slope log(J_B(r2)/J_B(r1))/log(r2/r1) at (50, 100) or (51, 101) below 1.35. K2 J_B(100)/J_F(100) < 0.4. CONTROLS, read first. C1 mean(F_e) by (i) against mpmath's quadrature of |F_e(e^(it))| at (r, J, e) = (5, 30, 2), (6, 33, 2), (8, 60, 2), relative error < 1e-12. C2 L by (ii) against /|w|_1 in integers at ranks 5..8, J = 20..40, every p <= 12 of J's parity; the parity check (a known zero, not a planted break): at J + p odd = 0 exactly. C3 mean >= exact height(F_e) and L <= exact H_ch at J_B, every rank read. C4 J_B <= J_F where J_F is known. THE ARMS. Ranks 5..14 every depth r - 1 .. J_c - 1; r = 20, 21, 30, 31, 40, 41, 50, 51, 100, 101, J scanned down from 6 r^(3/2) to the last beat (a beat above the start is excluded only by C4, where J_F is known). R is minimised over integer e (E|1 - 4eU| is convex in e) and, at odd n, over p by a scan to 4J/n + 12; mpmath at 30 digits. INSTRUMENT, before the counted run. C1 first read 2.4e-6: the quadrature did not split at the kinks of |1 - 4e sin^2|; split there it reads 4e-30. C3 first read False at r = 7, 41, 101, each where L/H_ch = 1: the central coefficient is the height and L, at 30 digits, sat a rounding above the exact integer; it now compares at a relative 1e-25. FINDINGS (copied from the printed output). (1) Controls: C1 4.14e-30, C2 1.52e-29, C3 and C4 hold at every rank read. C2's parity check never ran: J + n - p - q has the parity of J + p, so at J + p odd no k centres the phase and the loop skipped the cell first (0 executions, counted in a scratch copy); it is dropped. (2) K1 misses: the slopes over (50, 100) and (51, 101) are 1.5719 and 1.4769. K2 misses: J_B(100)/J_F(100) = 0.5518. The circle mean keeps the power r^(3/2). (3) Odd r (the central coefficient, L = H_ch at every J_B read): J_B = 18, 40, 64, 88, 128 at ranks 5, 7, 9, 11, 13 and 288, 508, 820, 1192, 3270 at r = 21, 31, 41, 51, 101; J_B/(r - 2)^(3/2) = 3.4775, 3.2529, 3.3668, 3.4752, 3.3197 at r = 21..101. P1's band holds, its rise does not: the ratio wanders with the integer e (4, 4, 5, 6, 8 at J_B). (4) Even r (the smoothed test, p = 3 .. 23): L/H_ch = 0.9583 .. 0.8207 at ranks 6..14 and 0.8269, 0.7992, 0.7882, 0.7895, 0.7771 at r = 20, 30, 40, 50, 100 (P4 holds there and at ranks 8..14; rank 6 reads 0.9583, above its band); J_B = 207, 365, 619, 815, 2423 there, J_B(r)/J_B(r + 1) = 0.6837 .. 0.7549 (P2's band, read against r + 1 alone: the arms hold no r - 1, so P2 as frozen, the mean of both neighbours, is not read), J_B/J_F = 0.5275 .. 0.6161 (P3 near). (5) The other slack: height(F_e)/mean = 0.6825 .. 0.8466 at every J_B. At ranks 5..14 J_B/J_F = 0.549 .. 0.762. THE ODD THEOREM (by hand, re-derived; the arm checks it). At odd r, x = r - 2, and even J, R < 1 is E|1 - cU| < 4 mu, c = 4e, mu = E U = a/N, N = a + b = (J + r - 1)/2. Take e nearest 1/(4 mu) = (J + 1 + x)/(4x) (at e = 0, mu > 1/2 and R < 1 at once; F_0 is nonzero), so |1 - c mu| <= 2 mu and c mu <= 1 + 2 mu; by Cauchy-Schwarz (E|1 - cU|)^2 <= (1 - c mu)^2 + c^2 Var U, and Var U = mu (1 - mu)/(N + 1) < mu/N gives c^2 Var U < (1 + 2 mu)^2/(mu N) <= 12 mu^2 once (1 + 2 mu)^2 N^2 <= 12 a^3. (1 + 2 mu) N = N + 2a is linear in J, so the condition is exactly J <= sqrt6 x^(3/2) - 3r + 5: every even J there fails the champion, and J_ch(r) > sqrt6 (r - 2)^(3/2) - 3r + 4 at every odd r. Hence J_ch(r) > 2 (r - 2)^(3/2) once g(x) = (sqrt6 - 2) x^(3/2) - 3x - 2 > 0: g increases past x = 19.8 and g(47) > 0, so at odd r >= 49. QUESTION of the arm, frozen with the kill before its code: does the exact second moment (rationals) pass at every even J <= J_hand? Q1 J_CS (best integer e, exact) /x^1.5 in [2.4, 2.85] at r = 201, rising (real c gives 2 sqrt2 = 2.83, TRANSPLANT: Var U/mu^2 ~ 1/a). Q2 J_CS/J_B in [0.7, 0.9]. Q3 J_hand/x^1.5 = 2.236 at r = 201. KILL: an even J <= J_hand where the exact second moment at the hand e fails, or R >= 1 at J_CS, over odd r = 3..201. C5 the exact second moment (frozen against quadrature; see the RUN RECORD). THE PLANT, frozen as Var U dropped, passes at every J by construction, |1 - c mu| <= 2 mu; it was replaced after the run by Var U quartered, whose scan's J must read R >= 1 at r = 5..201 (r = 3, a = 1/2, admitted after the run). It lands at 1.3 to 1.7 J_B, so it shows R catching a 30 to 70% overclaim, not one at J_CS's 0.75 to 0.89. (6) C5 1.18e-30. The kill misses: 0 failing J, 0 ranks with R >= 1 at J_CS; the quartered plant read R = 0.827 at r = 3 (its one miss) and 1.251 .. 1.593 at the printed r = 5..201, caught at 99 of 100, its J at the scan top at 0. J_hand/x^1.5 = 0.3849, 1.4801, 2.0175, 2.1461, 2.2364 at r = 5, 13, 51, 101, 201 (Q3 holds). J_CS/x^1.5 = 2.5101 .. 3.0792 over the ranks printed from 5, 2.8042 at r = 201 (Q1's band holds, its rise does not: the integer e again); J_CS/J_B = 0.7500 .. 0.8889 (Q2 holds). J_CS + 1 > 2 x^1.5 at every odd r = 5..201 (r = 3: J_CS = 0), so J_ch(r) > 2 (r - 2)^(3/2) at every odd r >= 5, exactly checked at r = 5..47 and proved past. g(47) = 1.8327. THE EXACT DEPTHS (by hand; exact_arm checks it). J = (4e - 1) x - 1 makes 1/(4 mu) = e exactly, so c mu = 1 and E(1 - cU)^2 = Var U/mu^2 < 1/(mu N) = 2/x. Odd r (J even): R < 1 once 2/x <= 16 mu^2, i.e. e <= sqrt(x/2), so J_ch(r) > (4 floor(sqrt(x/2)) - 1) x - 1 > 2 sqrt2 x^(3/2) - 5x - 1, with no computation. Even r (J odd): the test p = 4e - 1, 2k = 4e(x - 1), gives R <= 2^(-2-p) C(p, floor(p/2)) E|1 - cU|/rho, rho = B(a + 3/2, b + p/2)/ B(a, b); C(p, floor(p/2)) <= 2^p sqrt(2/(pi p)); log Gamma convex with log y - 1/y < psi(y) < log y gives y^t e^(-t/y) <= Gamma(y + t)/Gamma(y) <= (y + t)^t, so rho >= (a/N')^(3/2) (b/N')^(p/2) e^(-3/(2a) - p/(2b)), N' = N + 2e + 1, where 3/(2a) + p/(2b) = 4/x and (p/2) mu/(1 - mu) = 1/2 exactly. R < 1 once A(x, e): e sqrt(2/x) <= K sqrt(1 - 1/(4e)) (1 + delta)^(-3/2) exp(-(2e + 4.5)/x), K = sqrt(pi/2) exp(-1/2), delta = (2e + 1)/(2ex). At x >= 4e^2 the left side is at most 1/sqrt2 and the right least at x = 4e^2, where every factor rises in e: from the least e0 with RHS(e0) >= 1/sqrt2, e = floor(sqrt(x)/2) passes at every x >= 4 e0^2, and J_ch(r) > (4 floor(sqrt(x)/2) - 1) x - 1 > 2 x^(3/2) - 5x - 1. As x grows A passes every e < (K/sqrt2) sqrt(x) (1 - o(1)): liminf J_ch/r^(3/2) >= 2 sqrt2 at odd r and 2 sqrt(pi) exp(-1/2) = 2.150 at even r. Below x = 4 e0^2 the even claim is R_CS at 30 digits, its worst at most 0.9082: a computation, not a proof (the odd one is exact rationals). PREDICTIONS (frozen with the kill). S1 e0 in 8..11. S2 the exact-depth odd J/x^1.5 at r = 201 in [2.4, 2.83]. S3 at r = 100 A's largest e is 4 (J = 1469), R_CS's (the bound with the exact second moment, C and rho) 5 or 6. S4 R_CS < 1 at e = floor(sqrt(x)/2) at every even x < 4 e0^2. KILL: over odd r = 3..201 an exact depth with e <= sqrt(x/2) failing the exact second moment; over even r up to 4 e0^2 + 2 (frozen at 400, run to 486 once e0 read 11) an e with A true and R_CS >= 1, or an x < 4 e0^2 with R_CS >= 1 at floor(sqrt(x)/2). C6 rho's loggamma form against mpmath's beta; C7 R_CS >= R (30 digits) at e = floor(sqrt(x)/2); C8 the failure itself in integers, height(F_e) < H_ch at the exact depth, r = 5..41. THE PLANT: A with exp(-1/2) dropped; its largest e must read R_CS >= 1 at most even ranks. (7) C6 2.27e-28, C7 holds at r = 6..486, C8 37 of 37. The kill misses: 0 odd failures; 0 (x, e) with A true and R_CS >= 1 over even r = 4..486; R_CS at floor(sqrt(x)/2) at most 0.9082 over x < 484 (printed rounded up). e0 = 11 (RHS 0.70325 at 10, 0.70915 at 11; S1 holds at its edge). The plant caught at 216 of 238 ranks. Odd J/x^1.5 at the exact depth: 1.5396, 2.0832, 2.1399, 2.7126, 2.4807 at r = 5, 13, 51, 101, 201 (S2 holds). At r = 100 A's e is 4 (J 1469), R_CS's 5 (J 1861): S3 holds; R_CS's J/x^1.5 is 1.9183 .. 2.1363 at the printed r = 100, 200, 300, 400, 486. RUN RECORD. 34.1 s (the odd arm 4.5 s), peak commit 25.4 MB under a 512 MB guard; exact_arm 0.9 s, 21.2 MB, then behind a flag; since main runs exact_arm first, 35.2 s, 25.6 MB in all. C5 first compared with tanh-sinh quadrature and read 2e-1, then 7e-11 split at the peak: the quadrature, not the rationals, wandered at Beta(49.5, 1000.5) (INSTRUMENT). """ import math import os import sys import time from fractions import Fraction os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") import mpmath as mp # noqa: E402 import bandend # noqa: E402 import bandfamily # noqa: E402 mp.mp.dps = 30 JHI = bandfamily.JHI JC = bandfamily.JC JF = dict(JHI) JF.update({20: 336, 30: 692, 40: 1071, 50: 1525, 100: 4391}) # bandfamily.py LN2 = mp.log(2) def logB(x, y): return mp.loggamma(x) + mp.loggamma(y) - mp.loggamma(x + y) def mad(a, b, c): """E|1 - cU|, U ~ Beta(a, b).""" mu = a / (a + b) if c <= 1: return 1 - c * mu u0 = 1 / mp.mpf(c) i1 = mp.betainc(a, b, 0, u0, regularized=True) i2 = mp.betainc(a + 1, b, 0, u0, regularized=True) return 1 - c * mu - 2 * (1 - i1) + 2 * c * mu * (1 - i2) def log_mean(r, J, e): """log of F_e's circle mean.""" s = r - 3 a, b = mp.mpf(s + 1) / 2, mp.mpf(J + 1) / 2 return ((J + s) * LN2 - mp.log(mp.pi) + logB(a, b) + mp.log(mad(a, b, 4 * e))) def log_test(r, J, p): """log ||/|w|_1 at w = x^k (1 + x)^p (1 - x)^q, q = n mod 2.""" n = r - 1 q = n % 2 w1 = mp.mpf(2) ** p if q == 0 else 2 * mp.mpf(math.comb(p, p // 2)) return ((J + n + p + q) * LN2 - mp.log(mp.pi) + logB(mp.mpf(n + q + 1) / 2, mp.mpf(J + p + 1) / 2) - mp.log(w1)) def log_L(r, J): """(log L, p): the best test over p of J's parity.""" n = r - 1 p0 = J % 2 if n % 2 == 0: return log_test(r, J, p0), p0 best = None for p in range(p0, 4 * J // n + 13, 2): v = log_test(r, J, p) if best is None or v > best[0]: best = (v, p) return best def best_e(r, J): """(log min mean, e) over integer e >= 0, by convexity.""" s = r - 3 e0 = max(0, int(round((J + s + 2) / (4 * (s + 1))))) def f(e): return log_mean(r, J, e) if e >= 0 else mp.inf e, v = e0, f(e0) for step in (1, -1): while True: w = f(e + step) if w < v: e, v = e + step, w else: break return v, e def ratio(r, J): lm, e = best_e(r, J) ll, p = log_L(r, J) return float(mp.exp(lm - ll)), e, p def exact_height_F(r, J, e): V, U = bandfamily.exact_cols(bandend.Q(J, r - 3)) return max(abs(v + e * u) for v, u in zip(V, U)) def exact_Hch(r, J): return max(abs(c) for c in bandend.Q(J, r - 1)) def c1(): worst = 0.0 for r, J, e in ((5, 30, 2), (6, 33, 2), (8, 60, 2)): s = r - 3 def g(t): x = mp.exp(1j * t) return abs((1 + x) ** J * (1 - x) ** s * (x + e * (1 - x) ** 2)) t0 = 2 * mp.asin(1 / (2 * mp.sqrt(e))) # the kinks of |1 - 4e sin^2| brute = mp.quad(g, [0, t0, mp.pi, 2 * mp.pi - t0, 2 * mp.pi]) brute /= 2 * mp.pi worst = max(worst, float(abs(mp.exp(log_mean(r, J, e)) / brute - 1))) ok = worst < 1e-12 print("C1 Beta form of the circle mean against quadrature: worst" " relative error %.2e: %s" % (worst, ok), flush=True) return ok def c2(): worst = 0.0 for r in range(5, 9): n = r - 1 q = n % 2 for J in range(20, 41): P = bandend.Q(J, n) for p in range(0, 13): W = bandend.pmul([math.comb(p, i) for i in range(p + 1)], [1, -1] if q else [1]) if (J + n - p - q) % 2: continue k = (J + n - p - q) // 2 dot = sum(P[k + i] * W[i] for i in range(len(W))) w1 = sum(abs(x) for x in W) v = mp.mpf(abs(dot)) / w1 rel = abs(mp.exp(log_test(r, J, p)) / v - 1) worst = max(worst, float(rel)) ok = worst < 1e-20 print("C2 the test (ii) against /|w|_1 in integers, ranks 5..8," " J = 20..40, p <= 12: worst relative error %.2e: %s" % (worst, ok), flush=True) return ok def scan_small(r): beats = [J for J in range(r - 1, JC[r]) if ratio(r, J)[0] < 1] return max(beats) if beats else None def scan_large(r): J = int(6 * r ** 1.5) while J >= r - 1: R, e, p = ratio(r, J) if R < 1: return J, R, e, p J -= 1 return None def read(r, jb): R, e, p = ratio(r, jb) hf = exact_height_F(r, jb, e) hc = exact_Hch(r, jb) mean = mp.exp(best_e(r, jb)[0]) L = mp.exp(log_L(r, jb)[0]) # L equals H_ch when the central coefficient is the height: C2 holds # L to 1e-29, so the comparison carries a relative 1e-25 c3 = mean >= hf and L <= hc * (1 + mp.mpf(10) ** -25) return R, e, p, float(L / hc), float(hf / mean), c3 def cs_holds(x, J, e, plant=False): """Exact: E(1 - 4eU)^2 < 16 mu^2 at U ~ Beta(x/2, (J+1)/2); the plant quarters c^2 Var U (dropping it passes at every J by construction).""" a, b = Fraction(x, 2), Fraction(J + 1, 2) N = a + b mu = a / N c = 4 * e m2 = a * (a + 1) / (N * (N + 1)) if plant: return (1 - c * mu) ** 2 + c * c * (m2 - mu * mu) / 4 < 16 * mu ** 2 return 1 - 2 * c * mu + c * c * m2 < 16 * mu ** 2 def hand_e(x, J): """The integer nearest 1/(4 mu) = N/(4a) = (J + 1 + x)/(4x).""" return (2 * (J + 1 + x) + 4 * x) // (8 * x) def best_cs_e(x, J): """The integer e minimising the exact second moment (a quadratic).""" a, b = Fraction(x, 2), Fraction(J + 1, 2) N = a + b e0 = (a / N) / (4 * a * (a + 1) / (N * (N + 1))) return [math.floor(e0), math.floor(e0) + 1] def c5(): """The exact second moment against mpmath's Beta integrals at three (a, b, c).""" worst = 0.0 for x, J, e in ((5, 30, 2), (21, 300, 5), (99, 2000, 6)): a, b = Fraction(x, 2), Fraction(J + 1, 2) N = a + b c = 4 * e exact = 1 - 2 * c * a / N + c * c * a * (a + 1) / (N * (N + 1)) am, bm = mp.mpf(x) / 2, mp.mpf(J + 1) / 2 # mpmath's hypergeometric incomplete Beta, not a quadrature: at # Beta(49.5, 1000.5) tanh-sinh wanders by 1e-10 split at the peak q = sum(w * mp.betainc(am + k, bm, 0, 1) for k, w in ((0, 1), (1, -2 * c), (2, c * c))) q /= mp.beta(am, bm) worst = max(worst, float(abs(q / (mp.mpf(exact.numerator) / exact.denominator) - 1))) ok = worst < 1e-12 print("C5 the exact second moment against mpmath's Beta integrals: worst" " relative error %.2e: %s" % (worst, ok), flush=True) return ok def odd_arm(res): """THE ODD THEOREM's arm: every odd r = 3..201.""" t0 = time.time() fails_hand = 0 bad_R = 0 plant_caught = 0 miss = [] at_top = 0 nr = 0 rows = [] for r in range(3, 202, 2): x = r - 2 bound = math.sqrt(6) * x ** 1.5 - 3 * r + 5 jh = int(math.floor(bound)) jh -= jh % 2 for J in range(0, jh + 1, 2): if not cs_holds(x, J, hand_e(x, J)): fails_hand += 1 top = int(6 * x ** 1.5) + 12 top -= top % 2 jcs = next((J for J in range(top, -1, -2) if any(cs_holds(x, J, e) for e in best_cs_e(x, J))), None) jpl = next((J for J in range(top, -1, -2) if cs_holds(x, J, hand_e(x, J), plant=True)), None) Rcs = ratio(r, jcs)[0] Rpl = ratio(r, jpl)[0] bad_R += Rcs >= 1 plant_caught += Rpl >= 1 miss += [r] if Rpl < 1 else [] at_top += jpl == top nr += 1 rows.append((r, x, jh, jcs, Rcs, jpl, Rpl)) for r, x, jh, jcs, Rcs, jpl, Rpl in rows: if r <= 15 or r % 10 == 1: jb = res.get(r) print(" r = %3d J_hand %5d (/x^1.5 %.4f) J_CS %5d (/x^1.5" " %.4f, R %.4f)%s plant J %5d R %.3f" % (r, jh, jh / x ** 1.5, jcs, jcs / x ** 1.5, Rcs, " J_CS/J_B %.4f" % (jcs / jb) if jb else "", jpl, Rpl), flush=True) two = [r for r, x, jh, jcs, *_ in rows if jh + 1 <= 2 * x ** 1.5] two_cs = [r for r, x, jh, jcs, *_ in rows if jcs + 1 <= 2 * x ** 1.5] g = (math.sqrt(6) - 2) * 47 ** 1.5 - 3 * 47 - 2 print("odd r = 3..201 (%d ranks): even J <= J_hand where the exact CS" " with the hand e fails: %d; ranks with R >= 1 at J_CS: %d;" " the plant caught (R >= 1 at its J) at %d of %d, missed at r = %s," " its J at the scan top at %d" % (nr, fails_hand, bad_R, plant_caught, nr, miss, at_top), flush=True) print("J_hand + 1 <= 2 x^1.5 at r = %s; J_CS + 1 <= 2 x^1.5 at r = %s;" " g(47) = %.4f [%.1f s]" % (two, two_cs, g, time.time() - t0)) # r = 3 (a = 1/2) is the plant's one admitted miss return fails_hand == 0 and bad_R == 0 and miss in ([], [3]) K = math.sqrt(math.pi / 2) * math.exp(-0.5) def A(x, e, plant=False): """THE EVEN THEOREM's condition at J = (4e - 1) x - 1; the plant drops the factor e^(-1/2) the smoothing pays.""" delta = (2 * e + 1) / (2 * e * x) rhs = (K * math.sqrt(1 - 1 / (4 * e)) * (1 + delta) ** -1.5 * math.exp(-(2 * e + 4.5) / x)) if plant: rhs *= math.exp(0.5) return e * math.sqrt(2 / x) <= rhs def rhs_corner(e): """A's right side at x = 4e^2, where it is least over x >= 4e^2.""" return (K * math.sqrt(1 - 1 / (4 * e)) * (1 + (2 * e + 1) / (8 * e ** 3)) ** -1.5 * math.exp(-(2 * e + 4.5) / (4 * e * e))) def log_RCS(x, e): """log of the even rank's bound at J = (4e - 1) x - 1, p = 4e - 1: the exact second moment's root over the smoothed test, rho by loggamma.""" J, p = (4 * e - 1) * x - 1, 4 * e - 1 a, b = Fraction(x, 2), Fraction(J + 1, 2) N = a + b assert a / N == Fraction(1, 4 * e) c = 4 * e E2 = 1 - 2 * c * a / N + c * c * a * (a + 1) / (N * (N + 1)) am, bm = mp.mpf(x) / 2, mp.mpf(J + 1) / 2 lrho = logB(am + mp.mpf(3) / 2, bm + mp.mpf(p) / 2) - logB(am, bm) return ((-2 - p) * LN2 + mp.log(math.comb(p, p // 2)) + mp.log(mp.mpf(E2.numerator) / E2.denominator) / 2 - lrho) def exact_arm(): """The exact depths J = (4e - 1) x - 1 at both parities.""" t0 = time.time() # control: rho's loggamma form against mpmath's beta at three points worst = 0.0 for x, e in ((4, 1), (30, 3), (98, 5)): J, p = (4 * e - 1) * x - 1, 4 * e - 1 am, bm = mp.mpf(x) / 2, mp.mpf(J + 1) / 2 v = mp.exp(logB(am + 1.5, bm + p / 2) - logB(am, bm)) w = mp.beta(am + 1.5, bm + mp.mpf(p) / 2) / mp.beta(am, bm) worst = max(worst, float(abs(v / w - 1))) print("C6 rho by loggamma against mpmath's beta: worst %.2e: %s" % (worst, worst < 1e-20), flush=True) ok = worst < 1e-20 odd_fail = 0 for r in range(3, 202, 2): x = r - 2 e = math.isqrt(x // 2) if e < 1: continue J = (4 * e - 1) * x - 1 odd_fail += not cs_holds(x, J, e) if r in (5, 13, 51, 101, 201): print(" odd r = %3d e = %d, J = %d, J/x^1.5 = %.4f" % (r, e, J, J / x ** 1.5), flush=True) print("odd r = 3..201: exact depths with e <= sqrt(x/2) failing the" " exact second moment: %d" % odd_fail, flush=True) e0 = next(e for e in range(1, 10 ** 4) if rhs_corner(e) >= 2 ** -0.5) mono = all(rhs_corner(e + 1) >= rhs_corner(e) for e in range(1, 10 ** 4)) print("e0 = %d (RHS(e0 - 1) = %.5f, RHS(e0) = %.5f, 1/sqrt2 = %.5f);" " RHS rises over e = 1..10^4: %s" % (e0, rhs_corner(e0 - 1), rhs_corner(e0), 2 ** -0.5, mono), flush=True) a_bad = s4_bad = caught = nread = 0 ctl, worst_f = True, 0.0 for r in range(4, 4 * e0 * e0 + 3, 2): x = r - 2 top = 2 * math.isqrt(x) + 3 ep = max((e for e in range(1, top) if A(x, e, plant=True)), default=0) es = range(1, max(top, ep + 1)) lr = {e: log_RCS(x, e) for e in es} a_bad += sum(1 for e in es if A(x, e) and lr[e] >= 0) ea = max((e for e in es if A(x, e)), default=0) ee = max((e for e in es if lr[e] < 0), default=0) ef = math.isqrt(x) // 2 if ef >= 1: if x < 4 * e0 * e0: s4_bad += lr[ef] >= 0 worst_f = max(worst_f, float(mp.exp(lr[ef]))) J, p = (4 * ef - 1) * x - 1, 4 * ef - 1 ex = log_mean(r, J, ef) - log_test(r, J, p) ctl &= bool(lr[ef] >= ex) if ep: nread += 1 caught += lr[ep] >= 0 if r in (4, 6, 10, 14, 20, 50, 100, 200, 300, 400, 486): J = lambda e: (4 * e - 1) * x - 1 print(" even r = %3d e_A %d (J %d, /x^1.5 %.4f) e_RCS %d (J %d," " /x^1.5 %.4f) floor(sqrt x/2) = %d, R_CS there %.4f;" " plant e %d" % (r, ea, J(ea), J(ea) / x ** 1.5, ee, J(ee), J(ee) / x ** 1.5, ef, float(mp.exp(lr[ef])) if ef else 0, ep), flush=True) print("even r = 4..%d: e with A true and R_CS >= 1: %d; x < 4 e0^2 with" " R_CS >= 1 at floor(sqrt x/2): %d; C7 R_CS >= the exact R at" " floor(sqrt x/2), every rank: %s; the plant caught at %d of %d" " [%.1f s]" % (4 * e0 * e0 + 2, a_bad, s4_bad, ctl, caught, nread, time.time() - t0), flush=True) print("largest R_CS at floor(sqrt x/2) over x < 4 e0^2: %.4f (rounded" " up)" % (math.ceil(worst_f * 1e4) / 1e4)) # C8: the failure itself in integers at the exact depths, ranks 5..41 c8 = [] for r in range(5, 42): x = r - 2 e = math.isqrt(x // 2) if r % 2 else math.isqrt(x) // 2 if e < 1: continue J = (4 * e - 1) * x - 1 c8.append(exact_height_F(r, J, e) < exact_Hch(r, J)) print("C8 height(F_e) < H_ch in integers at the exact depth, r = 5..41:" " %d of %d" % (sum(c8), len(c8)), flush=True) # the plant must read R_CS >= 1 at most even ranks return (ok and ctl and all(c8) and 2 * caught > nread and not (odd_fail or a_bad or s4_bad)) def main(): t0 = time.time() if not exact_arm(): return True ok = c1() and c2() and c5() res = {} c3 = c4 = True print("ranks 5..14, every depth r - 1 .. J_c - 1", flush=True) for r in range(5, 15): jb = scan_small(r) R, e, p, lh, fm, c = read(r, jb) c3 &= c c4 &= jb <= JF[r] res[r] = jb print(" r = %2d J_B = %d (J_F %d, J_B/J_F %.3f), e = %d, p = %d," " R %.4f; L/H_ch %.4f, height(F)/mean %.4f; C3 %s [%.1f s]" % (r, jb, JF[r], jb / JF[r], e, p, R, lh, fm, c, time.time() - t0), flush=True) print("large r, J scanned down from 6 r^(3/2)", flush=True) for r in (20, 21, 30, 31, 40, 41, 50, 51, 100, 101): jb, R, e, p = scan_large(r) R, e, p, lh, fm, c = read(r, jb) c3 &= c if r in JF: c4 &= jb <= JF[r] res[r] = jb print(" r = %3d J_B = %d, J_B/(r-2)^1.5 = %.4f, J_B/r^1.5 = %.4f," " J_B/J_F %s, e = %d, p = %d, R %.4f; L/H_ch %.4f," " height(F)/mean %.4f; C3 %s [%.1f s]" % (r, jb, jb / (r - 2) ** 1.5, jb / r ** 1.5, "%.4f" % (jb / JF[r]) if r in JF else "-", e, p, R, lh, fm, c, time.time() - t0), flush=True) print("C3 mean >= height(F_e) and L <= H_ch at every J_B: %s" % c3) print("C4 J_B <= J_F wherever J_F is known: %s" % c4) for r in (20, 30, 40, 50, 100): print("P2 r = %d: J_B(r)/J_B(r + 1) = %.4f" % (r, res[r] / res[r + 1])) k1 = [] for r1, r2 in ((50, 100), (51, 101)): sl = math.log(res[r2] / res[r1]) / math.log(r2 / r1) print("K1 slope over (%d, %d): %.4f" % (r1, r2, sl)) if sl < 1.35: k1.append((r1, r2)) k2 = res[100] / JF[100] < 0.4 print("K1 %s; K2 J_B(100)/J_F(100) = %.4f (K2 %s)" % ("fires at %s" % k1 if k1 else "misses", res[100] / JF[100], "fires" if k2 else "misses")) odd = odd_arm(res) print("%.1f s" % (time.time() - t0)) return not (ok and c3 and c4 and odd) or bool(k1) or k2 if __name__ == "__main__": sys.exit(1 if main() else 0)