variation.py
Cited on: READING.
To run it: save variation.py, with edge.py, excess.py, heat.py, majorant.py, real_support.py, remainder.py, route.py beside it; needs mpmath, numpy, scipy.
"""
variation.py -- h's total variation D_1 bounded at every P >= 399 at
once, each piece a closed form or a P-free constant, and the larger
sup's excess X it leaves.
QUESTION Q. Is there a bound D1U(P) >= D_1 holding at every P >= 399,
each piece a closed form or a P-free constant, whose X stays below 1/4
- 0.0636 at every P >= 399? Notation as excess.py: T_u(j) = C(n, u)
2^-n j/(2n), n = 2u - j, j <= u; U = floor(P^2/16), M = UP.
BY HAND.
(Q1) The bumps' variation in closed form. For 1 <= j <= u - 1,
T(j+1)/T(j) = 2(u - j)(j + 1)/(j(2u - j - 1)) > 1 iff j^2 + j < 2u;
for j = -i, i >= 2, |T(-i+1)|/|T(-i)| > 1 iff i^2 - i > 2u. So T is
0 at j = 0, rises to its max at j* = min{j : j^2 + j >= 2u}, falls
to T(u) = 2^-(u+1) and jumps to 0; on j < 0 it is negative, |T|
largest at i* = max{i : i^2 - i <= 2u}, falling to 0. Exactly
V_u = 2 T(j*) + 2 |T(-i*)|.
Robbins with its correction negative (n >= m + 1) and nH(m/n) <= n
ln 2 - 2(m - n/2)^2/n give, for 1 <= u <= n - 1,
C(n, u) 2^-n <= sqrt(n/(2 pi u (n - u))) e^(-j^2/(2n)),
and sup_x x e^(-x^2/(2n)) = sqrt(n) e^(-1/2), so
u V_u <= K(u) = (2 pi e)^(-1/2) [(1 - j*/u)^(-1/2)
+ (1 + i*/u)^(-1/2)],
-> kappa. j*/u <= (1 + sqrt(1 + 8u))/(2u), falling in u, and the
second bracket term is below 1, so past a = 4e4
sum_(a < u <= U) V_u <= K_a ln(U/a),
K_a = (2 pi e)^(-1/2) [(1 - (1 + sqrt(1 + 8a))/(2a))^(-1/2) + 1]
= 0.484804 by hand; below a the table V_u (excess.py).
(Q2) The far bumps, sharper than excess.py's 2 sum |T|. Bump u = U + s
meets m <= M only at j <= -sP, and there i = -j has i^2 - i > 2u
(s^2 P^2 - sP > 2(U + s) as U <= P^2/16), so |T_u| is monotone on
that window: its variation over m in [0, M] is at most
|T_(U+s)(-sP)|. far' = sum_(s >= 1) |T_(U+s)(-sP)|, by (Q1)'s bound
with u(n - u) = u(u + i): |T_u(-i)| <= i/(2 sqrt(n)) (2 pi u (u +
i))^(-1/2) e^(-i^2/(2n)). By hand about 16 e^-4/(sqrt(pi) P^2),
1.0e-6 at P = 399, falling like P^-2; s past 50 under Hoeffding
(|T| <= e^(-i^2/(2n))/2) below e^(-2000).
(Q3) The edge, exactly: 1.5 + c_0 eps/beta, eps = 2 - beta = beta^-d,
c_0 = (1 - eps)/(2 - (d + 1) eps) in (0, 1); at most 1.5 +
(2/3)^d. L_s - L = 2 c_0 - 1 = (d - 1) eps/(2 - (d + 1) eps).
(Q4) The root tail, a P-free constant. Modes with theta = 2 pi k/P <=
0.35: the disc |mu| <= sigma theta^2/P (majorant.py W2, sigma read
at theta0 = 2 pi k/P0 for k <= 22 and at theta = 0.35 past, the
self-map's left side rising in theta, falling in P) holds the root,
and |P mu + theta^2 + i theta^3| <= rhobar, so Re(P mu) <= -theta^2
(1 - r), r = rhobar/theta^2 read at (theta0, P0) (rhobar/theta^2
rises in theta, falls in P). With |c| <= (theta/P)C1, C1 = sqrt(1 +
2.25 theta0^2) + E_c/P0; |1 - x| <= theta C2, C2 = expm1(theta0 +
mubar0)/theta0; 1/(1 - |x|) <= 1 + P/(theta^2 (1 - r)) (1 - e^-a >=
a/(1 + a)); U theta^2 >= (pi^2 k^2/4)(1 - 16/P^2):
mode k <= 2 C1 C2 [theta0^2/P0 + 1/(1 - r)]
exp(-(pi^2 k^2/4)(1 - 16/P0^2)(1 - r)),
whatever P >= P0. Modes with theta > 0.35: phi >= 2 pi k/d >
0.349125, Q >= Q0 = 5 - 4 cos 0.349125, |c| <= 2/P (|(d+1)x - 2d| >=
P), |1 - x| <= 2, at most d/2 of them: (4d/P)(1 + 2d/ln Q0)
e^(-UP ln Q0/(2d)), falling in P (e^-1072 at 399). Sum: R_r =
2 sum_k e^(-pi^2 k^2/4) (1 + 1%) = 0.1715 by hand (k = 1: C1
1.00095, C2 1.0079, r 3.5e-4), the exact tail 0.1698 at d = 400.
That the roots are one per sector, so the fixed points are all of
them, is (Q6).
(Q5) D1U(P) = edge + (table to U, K_a ln(U/a) past) + far' + R_r, so
past P = 800, D1U <= 2K_a ln P + 0.9069 + far' (2K_a = 0.969606,
U = floor(P^2/16) <= P^2/16). X(P) as excess.py
(Z1) with L >= W >= (P^2 - 1)/(4P) (E >= 0): X falls as L rises
((2/beta)D1U <= 2L + 1 - L_s), rises with D1U, so both substitutions
are safe. Past the read range, t1 >= (P/4 - 2)/D1U - 1 grows like
P/ln P and X <= P eps + 1.99^-t1 2.02 D1U falls.
(Q6) The roots, one per sector (written after the first two runs, read
by Q6's check in the third). On the closed unit disc, Re(2 - x) >= 1 and
|2 - x| >= 1, so g_k(x) = e^(2 pi i k/d) (2 - x)^(-1/d) (principal
branch) maps the disc into itself with |g_k'| <= 1/d: one fixed
point x_k each, k = 0..d - 1 (Banach). A root in the disc has x^d =
1/(2 - x), so it is x_k for exactly one k; x_0 = 1, and |x_k| = 1
forces |2 - x_k| = 1, so x_k = 1: the P roots inside are x_1..x_P,
one per sector, dphi - 2 pi k = arg(1/(2 - x)) in [-pi/6, pi/6]
(the largest |Im x|/(2 - Re x) on the disc is 1/sqrt 3, at cos phi
= 1/2). A root in majorant.py's disc about i theta_k, theta_k <=
0.35, has dphi - 2 pi k = theta_k + d Im mu, within 0.35 + sigma
theta^2 (d/P) < 0.6 of 0, so in no other sector's window (2 pi -
pi/6 away): it is x_k. The owed clause of (Q4) and of majorant.py's
W2 is this.
PREDICTIONS, CONTROL AND KILL (frozen before code).
Q0 CONTROL, read first. (a) (Q1)'s identity: 2T(j*) + 2|T(-i*)|
equals excess.py's V_table at every u = 1..4e4, 1e-12 relative;
PLANT: j* + 1 in its place fails at some u. (b) far''s monotone
window: the windowed variation of bump U + s over m in [0, M]
equals |T(-sP)| - |T(-uP)| at d = 40, 100, s = 1, 2, 3, to 1e-15
absolute. (c) X recomputed from excess.py's D_1^+ = 6.713490 and its
L at d = 400 gives its 2.049e-4.
Q1 K(u) >= u V_u at every u = 4..4e4; K(4e4)/(4e4 V) - 1 in [0, 1e-4];
K_a in [0.4847, 0.4849].
Q2 far' at P = 399 in [5e-7, 2e-6], below excess.py's far 1.1e-4.
Q3 R_r in [0.1698, 0.175], every mode bound above the exact root term
at d = 400, 800, 1600 (k <= 22, 30 digits not needed: doubles);
PLANT: the exponent's 16/P^2 read as 16/P fails at d = 400.
Q4 D1U(399) - D_1^+(399) in [0, 0.01] (D_1^+ = 6.713490).
Q5 X(P) at every P = 399..10^6: largest within t1 = 14's plateau, at
most 2.1e-4; D1U(P) <= 0.9696 ln P + 0.91 past P = 800.
KILL, refrozen from ROAD's (its "below the exact D_1^+" assumed the
far piece unchanged; (Q2) moves it, so the wrong-bound test is per
piece): (K1) a piece's P-uniform bound below its exact counterpart
at d = 400 or 800 (Sum V_u, far' against exact binomials, the root
tail, the edge), or D1U below the exact D_1 there; (K2) X(399) from
D1U at or above 1/4 - 0.0636, or X(P) at or above it at any P read.
AMENDED BEFORE THE RUN: Q3's plant as frozen moves the bound up (1 -
16/P widens the exponent's slack), so it could not fail; the plant
reads the cut as U = P^2/16 + P, the factor 1 + 16/P0.
AFTER THE SECOND RUN, frozen before its own: Q6 g_k's fixed points,
iterated 40 times from 0, equal heat.roots' x_k, k = 1..d/2, at d =
400 and 1600 to 1e-12; |d arg x_k - 2 pi k| <= pi/6 at both; the
window 0.35 + sigma(0.35) 0.35^2 (d/P) at P = 399 in [0.5, 0.6), by
hand 0.544; PLANT: g_(k+1)'s fixed points in x_k's place miss at the
float arm's 1e-9 bar.
AFTER THE FIRST RUN: Q0a's 1e-12 bar missed, 3.2e-10 at u = 36242;
against 40 digits the table and the closed form err there by -1.3e-10
and +1.8e-10, and by 1e-11 to 1e-10 both from u = 10^4 on (a 40-digit
diagnostic run beside the script, not kept): lgamma's
float wall near 8e5 (ulp 1.2e-10), not the identity. The float arm's
bar is now 1e-9, and the identity is read against a 40-digit brute
variation at u = 10, 100, 1000.
FINDINGS (copied from the second run's prints, Q6's from the fourth).
- Q0 holds: the identity to 3.2e-10 against the float table (the
wall) and to 2.3e-40 against 40 digits; j* + 1 misses by 0.23; the
far bumps' windowed variation equals the end value to 2.1e-25 at d
= 40, 100, s = 1, 2, 3; excess.py's X at d = 400 recomputes,
2.0488e-4.
- Q1 holds: K(u) >= u V_u at every u = 4..4e4, 0.666854 against
0.550781 at u = 4, 0.483951 against 0.483947 at 4e4 (7.6e-6); K_a =
0.484803, kappa 0.483941.
- Q2 holds: far' is 1.0910e-6 at P = 399 (exact 1.0907e-6), 2.6550e-7
at 799, P^2 far' 0.1737 and 0.1695: a hundred times below excess.py's
far.
- Q3 holds: R_r = 0.171483, mode 1's bound 0.171376 against the exact
0.169726, 0.169639, 0.169617 at d = 400, 800, 1600 (C1 1.000953, C2
1.007956, r 3.56e-4); the exact tails 0.169831, 0.169743, 0.169721;
past k = 22 and past theta = 0.35 below e^-1000 (e^-1063). The
plant (cut read as P^2/16 + P) falls below at modes 1 to 4.
- Q4 holds: D1U(399) = 6.715030, 0.00154 above excess.py's D_1^+ and
0.1397 above the exact D_1 (6.575283); at d = 800 7.387126 against
7.247336.
- Q5 holds: X is largest at P = 408, 2.0559e-4 (t1 = 14), 2.0493e-4 at
399, 5.50e-8 at 799, 1.70e-76 at 9999, 10^-622.7 at 10^5; D1U - 0.9696
ln P lies in [0.9069, 0.9081] at the P printed, at most 0.9070 past
800.
- The KILL does not fire: D1U, far' and the root tail sit above the
exact at d = 400, 800 (Sum V_u is the table there, the same numbers;
the edge's bound 1.5 + (2/3)^d is read as 1.5, below the doubles'
resolution), and X stays below 1/4 - 0.0636 at every P = 399..10^6.
Past 10^6 by (Q5) by hand, not run. Unproved: the table and lgamma
in doubles; terms below the doubles' resolution left out of the
sums, the edge's (2/3)^d, L_s - L's P (2/3)^d, and beta read as 2 in
eps_t1 = beta^-(t1+1) (B = beta exactly) and in 2/beta, which lower
X by a factor near 1 + (t1 + 2) 2^-(d+1).
- Q6 holds: g_k's fixed points equal heat.roots' x_k
to 7.3e-16 at d = 400 and 7.5e-16 at 1600; |d arg x_k - 2 pi k| is
at most 0.5228 and 0.5234 against pi/6 = 0.5236; sigma(0.35) =
1.5822, the window 0.5443; g_(k+1)'s fixed points miss x_k by at
least 1.6e-2 and 3.9e-3, the plant caught at each d (its nearest,
3.9e-3). The plant reads that the fixed points are indexed by k
alone: g_(k+1)'s is the root x_(k+1), caught by construction.
RUN RECORD: 15/15 checks; 60.0 s under the 512 MB guard, peak commit
162.1 MB. The first run stopped at its control, Q0a's 1e-12 bar
(above). The second read 13/13 in 59.5 s, 164.0 MB; K1's check was
renamed after it to the pieces it reads, its condition unchanged, and
Q6's two checks were added after it, frozen before the third run. The
third read 15/15 with the plant's nearest taken as the larger of the
two d's (1.6e-2), so a plant missed at one d passed; the fourth takes
the smaller.
"""
import math
import os
import sys
import time
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import numpy as np # noqa: E402
import mpmath as mp # noqa: E402
from heat import roots, coef # noqa: E402
from majorant import mode_const # noqa: E402
from excess import V_table, edge, excess # noqa: E402
from route import stream # noqa: E402
CHECKS, NAMES = [], []
MISSED = ()
P0, A_TAB, TH_CAP = 399, 40000, 0.35
G = (2 * math.pi * math.e) ** -0.5
BUDGET_X = 0.25 - 0.0636
D1P_399 = 6.713490 # excess.py, d = 400
def check(name, ok, detail=""):
CHECKS.append(bool(ok))
NAMES.append(name)
print(f" [{'PASS' if ok else 'FAIL'}] {name}" +
(f" -- {detail}" if detail else ""), flush=True)
def logT(u, j):
"""ln |T_u(j)| exactly (lgamma), j != 0, j <= u."""
n = 2 * u - j
return (math.lgamma(n + 1) - math.lgamma(u + 1) - math.lgamma(n - u + 1)
- n * math.log(2) + math.log(abs(j) / (2 * n)))
def logT_rob(u, j):
"""(Q1)'s Robbins bound on ln |T_u(j)|, 1 <= u <= n - 1."""
n = 2 * u - j
return (math.log(abs(j) / (2 * n)) + 0.5 * math.log(n / (2 * math.pi * u
* (n - u))) - j * j / (2 * n))
def jstar(u):
j = max(1, math.ceil((-1 + math.sqrt(1 + 8 * u)) / 2) - 1)
while j * j + j < 2 * u:
j += 1
return j
def istar(u):
i = math.floor((1 + math.sqrt(1 + 8 * u)) / 2) + 1
while i * i - i > 2 * u:
i -= 1
return i
def V_closed(u, shift=0):
j = min(jstar(u) + shift, u)
return 2 * math.exp(logT(u, j)) + 2 * math.exp(logT(u, -istar(u)))
def K_of(u):
return G * ((1 - jstar(u) / u) ** -0.5 + (1 + istar(u) / u) ** -0.5)
def K_a(a):
return G * ((1 - (1 + math.sqrt(1 + 8 * a)) / (2 * a)) ** -0.5 + 1)
def far_prime(P, exact=False):
"""sum_(s <= 50) |T_(U+s)(-sP)|, by the Robbins bound or exactly."""
U = P * P // 16
tot = 0.0
for s in range(1, 51):
u, i = U + s, s * P
assert i * i - i > 2 * u
tot += math.exp((logT if exact else logT_rob)(u, -i))
return tot
def mode_bound(k, th0, plant=False):
sg, mb, ec, rho = mode_const(th0, P0)
r = rho / th0 ** 2
C1 = math.sqrt(1 + 2.25 * th0 * th0) + ec
C2 = math.expm1(th0 + mb) / th0
cut = (1 + 16 / P0) if plant else (1 - 16 / P0 ** 2)
return (2 * C1 * C2 * (th0 * th0 / P0 + 1 / (1 - r))
* math.exp(-(math.pi * k) ** 2 / 4 * cut * (1 - r))), r, C1, C2
def q_part(P):
"""ln of the modes past theta = 0.35: (4d/P)(1 + 2d/ln Q0) e^(-UP lnQ0/(2d))."""
d, U = P + 1, P * P // 16
lq = math.log(5 - 4 * math.cos(TH_CAP * P0 / (P0 + 1)))
return (math.log(4 * d / P) + math.log1p(2 * d / lq)
- U * P * lq / (2 * d))
def root_const(plant=False):
K4 = int(TH_CAP * P0 / (2 * math.pi))
per = [mode_bound(k, 2 * math.pi * k / P0, plant)[0]
for k in range(1, K4 + 1)]
sg, mb, ec, rho = mode_const(TH_CAP, P0)
rm = rho / TH_CAP ** 2
C1 = math.sqrt(1 + 2.25 * TH_CAP ** 2) + ec
C2 = math.expm1(TH_CAP + mb) / TH_CAP
k = K4 + 1
ex = (math.pi * k) ** 2 / 4 * (1 - 16 / P0 ** 2) * (1 - rm)
ratio = math.exp(-(math.pi ** 2 / 4) * (2 * k + 1) * (1 - 16 / P0 ** 2)
* (1 - rm))
past = 2 * C1 * C2 * (TH_CAP ** 2 / P0 + 1 / (1 - rm)) * math.exp(-ex) \
/ (1 - ratio)
return per, past, rm
def exact_modes(d):
P = d - 1
U = P * P // 16
x, wt, b, _ = roots(d)
c = coef(x, d)
ax = np.abs(x)
return wt * np.abs(c) * np.abs(1 - x) * np.exp(U * P * np.log(ax)) / (1 - ax)
def section_control(V):
print("Q0 CONTROL", flush=True)
t0 = time.time()
err = max(abs(V_closed(u) / V[u] - 1) for u in range(1, A_TAB + 1))
check("Q0a V_u = 2T(j*) + 2|T(-i*)| at every u = 1..4e4 (1e-9 rel)",
err < 1e-9, f"largest {err:.1e} ({time.time() - t0:.0f} s)")
mp.mp.dps = 40
worst = 0
for u in (10, 100, 1000):
lo = -int(30 * math.sqrt(u)) - 100
T = [mp.binomial(2 * u - j, u) * mp.mpf(2) ** (j - 2 * u) * j
/ (2 * (2 * u - j)) for j in range(lo, u + 1)]
var = abs(T[0]) + sum(abs(T[k + 1] - T[k]) for k in range(len(T) - 1)) + abs(T[-1])
cl = 2 * T[jstar(u) - lo] + 2 * abs(T[-istar(u) - lo])
worst = max(worst, abs(var / cl - 1))
check("Q0a the identity by 40-digit brute variation at u = 10, 100, 1000",
worst < 1e-30, f"largest {mp.nstr(worst, 3)}")
errp = max(abs(V_closed(u, 1) / V[u] - 1) for u in range(1, A_TAB + 1))
check("Q0a PLANT: j* + 1 in its place fails", errp >= 1e-12,
f"largest {errp:.1e}")
worst = 0.0
for d in (40, 100):
P = d - 1
U = P * P // 16
for s in (1, 2, 3):
u = U + s
i = np.arange(u * P, s * P - 1, -1)
n = 2 * u + i
lg = np.array([math.lgamma(v + 1) for v in n]) - math.lgamma(u + 1) \
- np.array([math.lgamma(v - u + 1) for v in n]) \
- n * math.log(2) + np.log(i / (2.0 * n))
T = np.exp(lg)
var = np.abs(np.diff(T)).sum()
gap = abs(var - (T[-1] - T[0]))
worst = max(worst, gap)
print(f" d = {d}, s = {s}: windowed variation {var:.6e}, "
f"|T(-sP)| - |T(-uP)| {T[-1] - T[0]:.6e}", flush=True)
check("Q0b the far bumps' variation is the end value at d = 40, 100",
worst < 1e-15, f"largest gap {worst:.1e}")
L, E, h0, D, tail, b, integ = stream(400)
L = L + tail
X, XR, t0_, t1 = excess(400, L, h0, D1P_399, b)
check("Q0c X from excess.py's D_1^+ at d = 400 is its 2.049e-4",
abs(X - 2.049e-4) < 5e-8, f"{X:.4e}, t0 {t0_:.3f}")
return L, h0, b, D[1] + 2 * tail
def section_bumps(V):
print("Q1 THE BUMPS: K(u) >= u V_u, K_a past the table", flush=True)
bad = [u for u in range(4, A_TAB + 1) if K_of(u) < u * V[u] * (1 - 1e-13)]
u = A_TAB
rel = K_of(u) / (u * V[u]) - 1
ka = K_a(A_TAB)
for v in (4, 10, 100, 1000, 10000, 40000):
print(f" u = {v}: j* {jstar(v)}, i* {istar(v)}, u V_u "
f"{v * V[v]:.6f}, K(u) {K_of(v):.6f}", flush=True)
print(f" K_a at a = {A_TAB}: {ka:.6f}; kappa {2 * G:.6f}", flush=True)
check("Q1 K(u) >= u V_u at every u = 4..4e4, K(4e4) within 1e-4, K_a "
"in [0.4847, 0.4849]", not bad and 0 <= rel <= 1e-4 and
0.4847 <= ka <= 0.4849,
f"below at {bad[:5]}, K(4e4)/(u V) - 1 {rel:.2e}, K_a {ka:.6f}")
return ka
def section_far():
print("Q2 THE FAR BUMPS by their end values", flush=True)
rows = {}
for P in (399, 799):
fb, fe = far_prime(P), far_prime(P, exact=True)
rows[P] = (fb, fe)
print(f" P = {P}: far' bound {fb:.4e}, exact {fe:.4e}, P^2 far' "
f"{P * P * fb:.4f}", flush=True)
check("Q2 far' at P = 399 in [5e-7, 2e-6]", 5e-7 <= rows[399][0] <= 2e-6,
f"{rows[399][0]:.3e}")
return rows
def section_roots():
print("Q3 THE ROOT TAIL, a P-free constant", flush=True)
per, past, rm = root_const()
qp = q_part(P0)
R = sum(per) + past + math.exp(qp)
print(f" per mode k = 1, 2, 3: {per[0]:.6f} {per[1]:.3e} {per[2]:.3e};"
f" past k = 22 at theta 0.35 (r {rm:.4f}): {past:.1e}; theta >"
f" 0.35: e^{qp:.0f}", flush=True)
b1, r1, C1, C2 = mode_bound(1, 2 * math.pi / P0)
print(f" k = 1: C1 {C1:.6f}, C2 {C2:.6f}, r {r1:.2e}; R_r {R:.6f}, "
f"2 sum e^(-pi^2 k^2/4) "
f"{2 * sum(math.exp(-(math.pi * k) ** 2 / 4) for k in range(1, 9)):.6f}",
flush=True)
low, plant_low, tails = [], [], {}
pp, _, _ = root_const(plant=True)
for d in (400, 800, 1600):
ex = exact_modes(d)
tails[d] = float(ex.sum())
K4 = len(per)
low += [(d, k + 1) for k in range(K4) if per[k] < ex[k]]
plant_low += [(d, k + 1) for k in range(K4) if pp[k] < ex[k]]
print(f" d = {d}: exact tail {tails[d]:.6f}, mode 1 {ex[0]:.6f} "
f"(bound {per[0]:.6f}), modes past 22 {ex[K4:].sum():.1e}",
flush=True)
check("Q3 R_r in [0.1698, 0.175], every mode bound above the exact "
"term at d = 400, 800, 1600", 0.1698 <= R <= 0.175 and not low
and all(tails[d] <= R for d in tails),
f"R_r {R:.6f}, below at {low}")
check("Q3 PLANT: the cut read as P^2/16 + P fails at d = 400",
any(d == 400 for d, k in plant_low), f"below at {plant_low[:4]}")
return R, tails
def fixed_points(d, shift=0):
"""x_k = g_(k+shift)'s fixed point, iterated from 0, k = 1..d/2."""
w = np.exp(2j * np.pi * (np.arange(1, d // 2 + 1) + shift) / d)
x = np.zeros(d // 2, dtype=complex)
for _ in range(40):
x = w * (2 - x) ** (-1.0 / d)
return x
def section_sectors():
print("Q6 THE ROOTS, one per sector, by g_k's contraction", flush=True)
worst, plant = 0.0, math.inf
arg = 0.0
for d in (400, 1600):
x, _, _, _ = roots(d)
k = np.arange(1, d // 2 + 1)
err = float(np.abs(fixed_points(d) - x).max())
mis = float(np.abs(fixed_points(d, 1) - x).min())
dev = float(np.abs(d * np.angle(x) - 2 * np.pi * k).max())
worst, plant, arg = max(worst, err), min(plant, mis), max(arg, dev)
print(f" d = {d}: |fixed point - heat.roots| at most {err:.1e}; "
f"g_(k+1)'s at least {mis:.1e}; |d arg x_k - 2 pi k| at most "
f"{dev:.4f}", flush=True)
sg = mode_const(TH_CAP, P0)[0]
win = TH_CAP + sg * TH_CAP ** 2 * (P0 + 1) / P0
print(f" sigma(0.35) {sg:.4f}; the window 0.35 + sigma theta^2 (d/P) "
f"{win:.4f}; pi/6 {math.pi / 6:.4f}", flush=True)
check("Q6 g_k's fixed points are heat.roots' x_k at d = 400, 1600 "
"(1e-12), each within pi/6 of its sector, the window in [0.5, 0.6)",
worst < 1e-12 and arg <= math.pi / 6 and 0.5 <= win < 0.6,
f"largest {worst:.1e}, {arg:.4f}, window {win:.4f}")
check("Q6 PLANT: g_(k+1)'s fixed points miss x_k at the 1e-9 bar",
plant >= 1e-9, f"nearest {plant:.1e}")
def D1U_all(Ps, Vc, ka, R):
U = Ps * Ps // 16
bumps = np.where(U <= A_TAB, Vc[np.minimum(U, A_TAB)],
Vc[A_TAB] + ka * np.log(np.maximum(U, 1) / A_TAB))
far = np.zeros(len(Ps))
for s in range(1, 51):
u = (U + s).astype(float)
i = (s * Ps).astype(float)
n = 2 * u + i
far += np.exp(np.log(i / (2 * n)) + 0.5 * np.log(n / (2 * np.pi * u
* (u + i))) - i * i / (2 * n))
return 1.5 + bumps + far + R, bumps, far
def X_all(Ps, D):
L = (Ps * Ps - 1) / (4.0 * Ps)
Ls = L # + 2 c_0 - 1 <= P (2/3)^d
t1 = np.floor((Ls - 1) / D)
B = 2.0 # (1 - 2^-d)/(1/2), beta = 2 in doubles
lg = -t1 * math.log(2) - math.log(B) + np.log(np.minimum(2 * L + 1 - Ls, D))
return np.exp(lg), t1, lg / math.log(10)
def section_read(Vc, ka, R, far_rows, ctrl, tails):
print("Q4, Q5 and KILL: D1U and X at every P = 399..10^6", flush=True)
L400, h0, b, D1_400 = ctrl
_, _, _, D, tail, _, _ = stream(800)
D1_800 = D[1] + 2 * tail
Ps = np.arange(399, 10 ** 6 + 1, dtype=np.int64)
D, bumps, far = D1U_all(Ps, Vc, ka, R)
X, t1, lx = X_all(Ps, D)
for P in (399, 400, 799, 800, 1599, 3199, 9999, 10 ** 5, 10 ** 6):
k = P - 399
print(f" P = {P}: D1U {D[k]:.6f} (bumps {bumps[k]:.6f}, far' "
f"{far[k]:.2e}), - 0.9696 ln P {D[k] - 0.9696 * math.log(P):.4f};"
f" t1 {t1[k]:.0f}, X {X[k]:.3e} (log10 {lx[k]:.1f})", flush=True)
for d, Dex in ((400, D1_400), (800, D1_800)):
k = d - 1 - 399
fb, fe = far_rows[d - 1]
print(f" d = {d}: D1U {D[k]:.6f} against the exact D_1 {Dex:.6f};"
f" far' {fb:.3e} against {fe:.3e}; R_r {R:.6f} against "
f"{tails[d]:.6f}", flush=True)
wrong = [d for d, Dex in ((400, D1_400), (800, D1_800))
if D[d - 400] < Dex or far_rows[d - 1][0] < far_rows[d - 1][1]
or R < tails[d]]
check("KILL K1 not fired: D1U, far' and the root tail above the exact "
"at d = 400, 800", not wrong, f"wrong at {wrong}")
g = D[0] - D1P_399
check("Q4 D1U(399) - D_1^+(399) in [0, 0.01]", 0 <= g <= 0.01, f"{g:.5f}")
km = int(np.argmax(X))
past = Ps > 800
sl = float((D[past] - 0.9696 * np.log(Ps[past])).max())
check("Q5 X largest within t1 = 14's plateau, at most 2.1e-4; D1U <= "
"0.9696 ln P + 0.91 past 800", t1[km] == 14 and X[km] <= 2.1e-4
and sl <= 0.91, f"largest X {X[km]:.4e} at P = {Ps[km]} (t1 "
f"{t1[km]:.0f}), D1U - 0.9696 ln P at most {sl:.4f}")
over = Ps[X >= BUDGET_X]
check("KILL K2 not fired: X below 1/4 - 0.0636 at every P read",
X[0] < BUDGET_X and len(over) == 0,
f"X(399) {X[0]:.4e}, at or above at {over[:5]}")
def main():
t0 = time.time()
V = V_table(A_TAB)
Vc = np.cumsum(V)
ctrl = section_control(V)
if not all(ok for n, ok in zip(NAMES, CHECKS) if n not in MISSED):
print("control failed; stopping")
sys.exit(1)
ka = section_bumps(V)
far_rows = section_far()
R, tails = section_roots()
section_sectors()
section_read(Vc, ka, R, far_rows, ctrl, tails)
print(f"{sum(CHECKS)}/{len(CHECKS)} checks passed, {time.time() - t0:.0f} s")
sys.exit(as_recorded(MISSED))
def as_recorded(missed):
"""0 iff the checks that failed are exactly the record's misses, by name."""
got = {n for n, ok in zip(NAMES, CHECKS) if not ok}
print("misses %s; the record's %s: %s" % (
"; ".join(sorted(got)) or "none", "; ".join(sorted(missed)) or "none",
"as recorded" if got == set(missed) else "NOT AS RECORDED"))
return 0 if got == set(missed) else 1
if __name__ == "__main__":
main()