circlemean.py
Checks: FLATTEN: The takeover depth grows at least as r^(3/2) theorem; rule; pattern.
Cited on: FLATTEN.
To run it: save circlemean.py, with band.py, bandend.py, bandfamily.py, flatten.py beside it; needs mpmath, numpy, sympy.
"""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, <P, w> = (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: |<P, w>| = (2^(J+n+p+q)/pi) B((n+q+1)/2, (J+p+1)/2),
and H_ch >= |<P, w>|/|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 <P, w>/|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 <P, w> = 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 |<P, w>|/|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 <P, w>/|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)