verify.py
10.8 kB · python · 319 lines
1import math2import time3from array import array4from operator import mul56# TRANSFORM78def factors(q, a0, ts):9 inv = 1.0 / (q - 1)10 c = 2 * a0 - q + 111 sin = math.sin12 cos = math.cos13 sqrt = math.sqrt14 pi = math.pi15 out = []16 for t in ts:17 d = sin(pi * q * t) / sin(pi * t)18 e = d * d - 2.0 * d * cos(pi * c * t) + 1.019 out.append(sqrt(e) * inv if e > 0.0 else 0.0)20 return out2122def direct(q, a0, t):23 re = 0.024 im = 0.025 for a in range(q):26 if a != a0:27 re += math.cos(2.0 * math.pi * a * t)28 im += math.sin(2.0 * math.pi * a * t)29 return math.hypot(re, im) / (q - 1)3031def closed_form(pts=1500):32 worst = 0.033 for q, a0 in ((7, 0), (10, 5), (21, 0), (21, 10), (34, 16)):34 ts = [(j + 0.37) / pts for j in range(pts)]35 for t, v in zip(ts, factors(q, a0, ts)):36 worst = max(worst, abs(v - direct(q, a0, t)))37 assert worst < 1e-11, f"closed form off by {worst}"38 print(f"closed form against the character sum, five sets, {pts} arguments each: worst gap {worst:.2e}")3940def one_digit_sum():41 for q in range(3, 40):42 for a0 in range((q + 1) // 2):43 ts = [a / q for a in range(1, q)]44 total = 1.0 + sum(factors(q, a0, ts))45 assert abs(total - 2.0) < 1e-9, f"one-digit sum {total} at base {q} digit {a0}"46 print("one-digit grid sum equals 2 exactly, every base 3 to 39 and every distinct digit")4748def domination(pts=4000):49 slack = 0.050 for q in (3, 10, 21, 34, 126, 200):51 ts = [(j + 0.5) / pts for j in range(pts)]52 cap = [min(1.0, (abs(math.sin(math.pi * q * t) / math.sin(math.pi * t)) + 1.0) / (q - 1)) for t in ts]53 for a0 in range((q + 1) // 2):54 for v, u in zip(factors(q, a0, ts), cap):55 assert v <= u + 1e-12, f"domination fails at base {q} digit {a0}"56 slack = max(slack, u - v)57 print(f"domination holds on 6 bases, every distinct digit, {pts} arguments: widest slack {slack:.4f}")5859# WINDOW MACHINE6061def weights(q, a0, nd, m, guard=1e-12):62 W = q ** nd63 c = 2 * a0 - q + 164 inv = 1.0 / (q - 1)65 lip = 2.0 * math.pi * (q * (q - 1) // 2 - a0) * inv66 slack = lip / (2.0 * W * m) + guard67 hi = array("d", bytes(8 * W))68 lo = array("d", bytes(8 * W))69 step = 1.0 / (W * m)70 sin = math.sin71 cos = math.cos72 sqrt = math.sqrt73 pi = math.pi74 for w in range((W + 1) // 2):75 top = 0.076 bot = 2.077 base = w * m78 for r in range(m):79 t = (base + r + 0.5) * step80 d = sin(pi * q * t) / sin(pi * t)81 e = d * d - 2.0 * d * cos(pi * c * t) + 1.082 v = sqrt(e) * inv if e > 0.0 else 0.083 if v > top:84 top = v85 if v < bot:86 bot = v87 u = top + slack88 u = u if u < 1.0 else 1.089 b = bot - slack90 b = b if b > 0.0 else 0.091 hi[w] = u92 lo[w] = b93 hi[W - 1 - w] = u94 lo[W - 1 - w] = b95 return hi, lo, slack9697def iterate(G, q, nd, seed, steps):98 S = q ** (nd - 1)99 P = q ** (nd - 2)100 y = list(seed) if seed else [1.0] * S101 lam = 0.0102 for _ in range(steps):103 z = [sum(map(mul, G[v * q:v * q + q], y[(v % P) * q:(v % P) * q + q])) for v in range(S)]104 top = max(z)105 if top <= 0.0:106 return 0.0, y107 y = [x / top for x in z]108 lam = top109 return lam, y110111def lift(y, q):112 return [y[v // q] for v in range(len(y) * q)]113114def cw_up(G, y, q, nd):115 S = q ** (nd - 1)116 P = q ** (nd - 2)117 best = 0.0118 for v in range(S):119 b = (v % P) * q120 r = sum(map(mul, G[v * q:v * q + q], y[b:b + q])) / y[v]121 if r > best:122 best = r123 return best124125def cw_low(G, y, q, nd):126 S = q ** (nd - 1)127 P = q ** (nd - 2)128 best = 0.0129 for tau in (0.0, 1e-13, 1e-9, 1e-5, 1e-3):130 low = float("inf")131 seen = 0132 for v in range(S):133 if y[v] <= tau:134 continue135 b = (v % P) * q136 acc = 0.0137 for c in range(q):138 if y[b + c] > tau:139 acc += G[v * q + c] * y[b + c]140 r = acc / y[v]141 if r < low:142 low = r143 seen += 1144 if seen and low > best:145 best = low146 return best147148def climb(q, a0, nd, m, side):149 start = min(nd, 3)150 pair = weights(q, a0, start, m)151 G = pair[0] if side else pair[1]152 _, y = iterate(G, q, start, None, 300)153 for level in range(start + 1, nd + 1):154 pair = weights(q, a0, level, m)155 G = pair[0] if side else pair[1]156 _, y = iterate(G, q, level, lift(y, q), 30)157 return G, y158159def upper(q, a0, nd, m):160 G, y = climb(q, a0, nd, m, True)161 return math.log(cw_up(G, y, q, nd)) / math.log(q)162163def lower(q, a0, nd, m):164 G, y = climb(q, a0, nd, m, False)165 mu = cw_low(G, y, q, nd)166 return None if mu <= 0.0 else math.log(mu) / math.log(q)167168# THE LADDER169170def below_the_floor():171 t0 = time.time()172 lanes = 0173 worst = (1.0, None)174 used = {}175 for q in range(3, 21):176 for a0 in range((q + 1) // 2):177 got = None178 for nd in (2, 3, 4):179 e = lower(q, a0, nd, 8)180 if e is not None and e > 0.2502:181 got = (e, nd)182 break183 assert got, f"no lower bound above a quarter at base {q} digit {a0}"184 used[got[1]] = used.get(got[1], 0) + 1185 lanes += 1186 if got[0] < worst[0]:187 worst = (got[0], (q, a0, got[1]))188 q, a0, nd = worst[1]189 print(f"every one of the {lanes} distinct sets of every base 3 to 20 certifies alpha_1 > 1/4")190 print(f" windows used: {used}; tightest base {q} digit {a0} at {nd} window digits, alpha_1 > {math.floor(worst[0] * 1e7) / 1e7:.7f}")191 print(f" {time.time() - t0:.1f}s")192193def the_first_base():194 t0 = time.time()195 four = (lower(21, 0, 4, 8), upper(21, 0, 4, 8))196 assert four[0] < 0.25 < four[1], "four window digits should not decide base 21"197 five = upper(21, 0, 5, 2)198 assert five < 0.25, f"base 21 digit 0 reads {five} at five window digits"199 print(f"base 21 missing 0: four window digits give [{math.floor(four[0] * 1e7) / 1e7:.7f}, {math.ceil(four[1] * 1e7) / 1e7:.7f}], undecided")200 print(f" five window digits give alpha_1 < {math.ceil(five * 1e7) / 1e7:.7f}, clearing 1/4 by {0.25 - five:.2e}")201 print(f" {time.time() - t0:.1f}s")202203WITNESS = ((21, 9), (22, 10), (23, 7), (24, 11), (25, 11), (26, 8), (27, 12), (28, 13), (29, 9), (30, 14), (31, 14), (32, 10), (33, 15))204205def the_family_floor():206 t0 = time.time()207 print("every base 21 to 33 carries an excluded digit with alpha_1 > 1/4, one witness each")208 for q, a0 in WITNESS:209 nd = 4 if q == 33 else 3210 e = lower(q, a0, nd, 8)211 assert e is not None and e > 0.25, f"the witness digit {a0} of base {q} reads {e}"212 print(f" base {q} missing {a0}: alpha_1 > {math.floor(e * 1e7) / 1e7:.7f} at {nd} window digits")213 tall = (0.0, None)214 for a0 in range(17):215 e = upper(34, a0, 3, 8)216 nd = 3217 if e >= 0.25:218 e = upper(34, a0, 4, 2)219 nd = 4220 assert e < 0.25, f"base 34 digit {a0} reads {e}"221 print(f" base 34 missing {a0}: alpha_1 < {math.ceil(e * 1e7) / 1e7:.7f} at {nd} window digits")222 if e > tall[0]:223 tall = (e, (a0, nd))224 a0, nd = tall[1]225 print(f"all 17 distinct sets of base 34 certify alpha_1 < 1/4, the largest bound at digit {a0} and {nd} window digits, alpha_1 < {math.ceil(tall[0] * 1e7) / 1e7:.7f}")226 print(f" {time.time() - t0:.1f}s")227228def the_calibration():229 t0 = time.time()230 lo = lower(10, 5, 5, 8)231 hi = upper(10, 5, 5, 8)232 assert lo is not None and hi < 27.0 / 77.0, f"base 10 missing 5 reads [{lo}, {hi}]"233 print(f"base 10 missing 5 at five window digits: alpha_1 in [{math.floor(lo * 1e7) / 1e7:.7f}, {math.ceil(hi * 1e7) / 1e7:.7f}], under 27/77 = {27 / 77:.7f}")234 print(f" {time.time() - t0:.1f}s")235236# THE LEBESGUE CONSTANT237238GAMMA = 0.5772156649015328606239C1 = 2.0 / math.pi240C0 = 0.98241242def lebesgue(M):243 n = (M + 1) // 2244 total = 2.0 * sum(1.0 / math.sin((2 * j + 1) * math.pi / (2 * M)) for j in range(n))245 return total - (M % 2)246247def lebesgue_scan(offsets=97):248 top = 0.0249 for M in list(range(2, 200)) + [256, 400, 1000, 2048, 4096]:250 best = 0.0251 for k in range(1, offsets + 1):252 th = k / (2.0 * offsets)253 s = abs(math.sin(math.pi * th))254 acc = 0.0255 for j in range(M):256 d = (j + th) / M257 d = min(d, 1.0 - d) if d <= 1.0 else d258 acc += 1.0 / math.sin(math.pi * min(d, 1.0 - d))259 best = max(best, s * acc)260 exact = lebesgue(M)261 assert best <= exact + 1e-9, f"the half offset is not the maximum at M = {M}"262 cap = M * (C1 * math.log(M) + GAMMA * C1 + C1 * math.log(8.0 / math.pi)) + 2.0 / math.pi263 assert exact <= cap, f"Lebesgue bound fails at M = {M}"264 if M >= 100:265 top = max(top, (exact / M - C1 * math.log(M)))266 assert top < C0, f"the uniform constant {top} does not sit under {C0}"267 print(f"Lebesgue constant: the half offset maximises at every M scanned, and L_M/M - (2/pi) log M <= {top:.7f} < {C0} for M >= 100")268269# THE THRESHOLD270271def margin(q, c0=C0):272 w = q ** 0.25 * (1.0 - 1.0 / q)273 return (w - 1.0) ** 3 - C1 * math.log(q) * w - c0 * (w - 1.0)274275def root(q, c0=C0):276 lo, hi = 1.0 + 1e-15, 64.0277 for _ in range(200):278 mid = 0.5 * (lo + hi)279 if (mid - 1.0) ** 3 - C1 * math.log(q) * mid - c0 * (mid - 1.0) < 0.0:280 lo = mid281 else:282 hi = mid283 return hi284285def threshold():286 assert margin(125) < 0.0, "the certificate should fail at 125"287 thin = (1.0, None)288 for q in range(126, 3000):289 m = margin(q)290 assert m > 0.0, f"the certificate fails at {q}"291 if m < thin[0]:292 thin = (m, q)293 for q in (211, 500, 3000, 10 ** 6, 10 ** 12):294 cap = 1.0 + math.sqrt(2.0 * C1 * math.log(q) + C0)295 assert cap < q ** 0.25 * (1.0 - 1.0 / q), f"the closed-form cap fails at {q}"296 print(f"threshold: the certificate fails at 125 and holds on [126, 3000), tightest margin {thin[0]:.3e} at q = {thin[1]}")297 for q in (126, 200, 1000, 10 ** 6):298 e = math.log(root(q) * q / (q - 1.0)) / math.log(q)299 assert e < 0.25 or q < 126, "the exponent bound should clear a quarter"300 print(f" q = {q}: alpha_1 < {math.ceil(e * 1e6) / 1e6:.6f}")301 assert abs(GAMMA * C1 + C1 * math.log(8.0 / math.pi) - 0.9625228) < 1e-7302303# VERIFY304305def main():306 t0 = time.time()307 closed_form()308 one_digit_sum()309 domination()310 lebesgue_scan()311 threshold()312 below_the_floor()313 the_first_base()314 the_family_floor()315 the_calibration()316 print(f"first-base-below-a-quarter: every check green in {time.time() - t0:.1f}s")317318if __name__ == "__main__":319 main()