certify.py
16.9 kB · python · 542 lines
1import sys2import time3from fractions import Fraction4from math import factorial, isqrt56# ROUNDING78PREC = 192910TOL = 1e-521112def _scale(x):13 return PREC - (abs(x.numerator).bit_length() - x.denominator.bit_length())1415def down(x):16 if x == 0:17 return Fraction(0)18 e = _scale(x)19 if e >= 0:20 return Fraction((x.numerator << e) // x.denominator, 1 << e)21 return Fraction((x.numerator // (x.denominator << -e)) << -e)2223def up(x):24 return -down(-x)2526# INTERVALS2728def num(x):29 f = Fraction(x)30 return (down(f), up(f))3132ZERO = (Fraction(0), Fraction(0))3334ONE = (Fraction(1), Fraction(1))3536def lo(a):37 return a[0]3839def hi(a):40 return a[1]4142def iadd(a, b):43 return (down(a[0] + b[0]), up(a[1] + b[1]))4445def isub(a, b):46 return (down(a[0] - b[1]), up(a[1] - b[0]))4748def ineg(a):49 return (-a[1], -a[0])5051def imul(a, b):52 c = (a[0] * b[0], a[0] * b[1], a[1] * b[0], a[1] * b[1])53 return (down(min(c)), up(max(c)))5455def idiv(a, b):56 assert b[0] > 0 or b[1] < 0, "interval division straddles zero"57 c = (a[0] / b[0], a[0] / b[1], a[1] / b[0], a[1] / b[1])58 return (down(min(c)), up(max(c)))5960def _rpow(x, n, rnd):61 r = Fraction(1)62 b = x63 while n:64 if n & 1:65 r = rnd(r * b)66 n >>= 167 if n:68 b = rnd(b * b)69 return r7071def ipow(a, n):72 assert a[0] >= 0, "interval power needs a nonnegative base"73 return (_rpow(a[0], n, down), _rpow(a[1], n, up))7475def amax(a):76 return max(abs(a[0]), abs(a[1]))7778def relwidth(a):79 m = amax(a)80 assert m > 0, "relative width of an interval containing only zero"81 return float((a[1] - a[0]) / m)8283def below(a, b):84 return a[1] < b[0]8586# SQUARE ROOTS8788def isqrt_interval(m):89 r = isqrt(m << (2 * PREC))90 return (Fraction(r, 1 << PREC), Fraction(r + 1, 1 << PREC))9192SQRT2 = isqrt_interval(2)9394SQRT3 = isqrt_interval(3)9596# PI BY MACHIN9798def _arctan_reciprocal(n, terms):99 s = Fraction(0)100 p = Fraction(1, n)101 sq = Fraction(1, n * n)102 for k in range(terms):103 t = p / (2 * k + 1)104 s = s + t if k % 2 == 0 else s - t105 p *= sq106 e = p / (2 * terms + 1)107 return (down(s - e), up(s + e))108109PI = isub(imul(num(16), _arctan_reciprocal(5, 80)),110 imul(num(4), _arctan_reciprocal(239, 32)))111112PI2 = imul(PI, PI)113114# LOGARITHM OF TWO115116def _log_two():117 z = Fraction(1, 3)118 terms = 130119 s = Fraction(0)120 p = z121 sq = z * z122 for k in range(terms):123 s += p / (2 * k + 1)124 p *= sq125 e = p / ((2 * terms + 1) * (1 - sq))126 return (down(2 * (s)), up(2 * (s + e)))127128LOG2 = _log_two()129130# SERIES131132def _cos_order(a):133 x = float(a)134 n = 1135 t = x * x / 2.0136 while t > TOL and n < 200:137 n += 1138 t = t * x * x / ((2 * n - 1) * (2 * n))139 return n + 2140141def iversin(t):142 a = amax(t)143 n = _cos_order(a)144 assert a * a <= (2 * n + 3) * (2 * n + 4), "versine tail is not decreasing"145 tt = imul(t, t)146 p = ONE147 s = ZERO148 for k in range(1, n + 1):149 p = imul(p, tt)150 term = idiv(p, num(factorial(2 * k)))151 s = iadd(s, term) if k % 2 else isub(s, term)152 p = imul(p, tt)153 e = hi(idiv(p, num(factorial(2 * n + 2))))154 assert e >= 0, "versine tail bound went negative"155 return iadd(s, (-e, e))156157def icos_grid(j, m):158 j %= m159 if 2 * j > m:160 j = m - j161 if j == 0:162 return ONE163 return isub(ONE, iversin(imul(PI, num(Fraction(2 * j, m)))))164165def _geom_order(a):166 x = float(a)167 if x <= 0:168 return 1169 n = 1170 t = x171 while t > TOL and n < 400:172 n += 1173 t *= x174 return n + 2175176def ilog1p(p):177 a = amax(p)178 assert a < Fraction(1, 2), "log1p outside the reduced window"179 if a == 0:180 return ZERO181 n = _geom_order(a)182 q = ONE183 s = ZERO184 for k in range(1, n + 1):185 q = imul(q, p)186 term = idiv(q, num(k))187 s = iadd(s, term) if k % 2 else isub(s, term)188 e = up(_rpow(a, n + 1, up) / ((n + 1) * (1 - a)))189 return iadd(s, (-e, e))190191def ilog(x):192 assert x[0] > 0, "logarithm of a nonpositive interval"193 m = 0194 while x[0] < Fraction(2, 3):195 x = imul(x, num(2))196 m -= 1197 while x[1] > Fraction(4, 3):198 x = imul(x, num(Fraction(1, 2)))199 m += 1200 return iadd(ilog1p(isub(x, ONE)), imul(num(m), LOG2))201202def _exp_order(a):203 x = float(a)204 n = 1205 t = x206 while t > TOL and n < 400:207 n += 1208 t = t * x / n209 return n + 4210211def iexp(x):212 assert x[0] >= 0 and x[1] < 30, "exponential outside the certified window"213 n = _exp_order(amax(x))214 assert x[1] < n, "exponential tail is not geometric"215 p = ONE216 s = ONE217 for k in range(1, n + 1):218 p = idiv(imul(p, x), num(k))219 s = iadd(s, p)220 step = idiv(imul(p, x), num(n + 1))221 e = hi(idiv(step, isub(ONE, idiv(x, num(n + 2)))))222 assert e >= 0, "exponential tail bound went negative"223 return iadd(s, (Fraction(0), e))224225# BASE THREE SPECTRUM226227TAIL_LEVEL = 20228229def _versin_grid_three(kind, i):230 return iversin(imul(PI, num(Fraction(kind, 3 ** i))))231232ALPHA = {}233234BETA = {}235236for _i in range(2, TAIL_LEVEL + 1):237 ALPHA[_i] = imul(num(2), _versin_grid_three(1, _i))238 BETA[_i] = imul(num(2), _versin_grid_three(2, _i))239240def ell(i):241 a = ilog1p(ineg(idiv(ALPHA[i], num(2))))242 b = ilog1p(ineg(idiv(BETA[i], num(2))))243 return isub(a, b)244245VERS_TOP = iversin(imul(PI, num(Fraction(2, 9))))246247LOG_GUARD = idiv(ONE, isub(ONE, VERS_TOP))248249ELL_COEFF = imul(imul(num(2), PI2), LOG_GUARD)250251def _tail_ell(n):252 assert n >= 1, "the ell tail bound needs level at least one"253 return (Fraction(0), hi(imul(ELL_COEFF, num(Fraction(1, 8 * 9 ** n)))))254255def _log_r():256 s = ZERO257 for i in range(2, TAIL_LEVEL + 1):258 s = iadd(s, ell(i))259 return iadd(s, _tail_ell(TAIL_LEVEL))260261LOGR = _log_r()262263S_GUARD = idiv(ONE, isub(ONE, idiv(imul(num(2), VERS_TOP), num(40))))264265S_COEFF = imul(PI2, iadd(ONE, imul(num(4), S_GUARD)))266267def _tail_s(D, n):268 assert n >= 1, "the s tail bound needs level at least one"269 assert D >= 38, "the s tail guard assumes D at least 38"270 piece = idiv(S_COEFF, num(D - 2))271 return (Fraction(0), hi(imul(piece, num(Fraction(1, 8 * 9 ** n)))))272273def s_of(D):274 s = ZERO275 for i in range(2, TAIL_LEVEL + 1):276 p = idiv(ALPHA[i], num(D - 2))277 q = idiv(BETA[i], num(D + 2))278 s = iadd(s, isub(ilog1p(p), ilog1p(ineg(q))))279 return iadd(s, _tail_s(D, TAIL_LEVEL))280281def a_of(D):282 return ilog1p(num(Fraction(4, D - 2)))283284def k_star(D):285 return idiv(iadd(imul(num(D - 1), LOGR), s_of(D)), a_of(D))286287def k_one(D):288 top = hi(iadd(imul(num(D - 1), LOGR), s_of(D)))289 bot = lo(a_of(D))290 q = top / bot291 assert abs(q - round(q)) > Fraction(1, 10 ** 12), "K1 sits on an integer boundary at D=%d" % D292 return int(q) + 2293294# THE ENVELOPE295296C_PI = num(Fraction(7528157, 10 ** 7))297298C_ZERO = num(Fraction(7052518, 10 ** 7))299300C_TWO = num(Fraction(2266816, 10 ** 7))301302C_SUB = num(Fraction(8900159, 10 ** 7))303304def envelope(D, K):305 sigma = iexp(imul(num(2 * (K + 1)), ipow(C_SUB, D - 1)))306 head = imul(num(4 * (K - 1)), iadd(ipow(C_PI, D - 1), ipow(C_ZERO, D - 1)))307 body = iadd(head, imul(num(2), ipow(C_TWO, D - 1)))308 return imul(body, sigma)309310def delta_of(D, K):311 t = imul(PI, num(Fraction(D - 2, 3 ** (K + 1))))312 return (Fraction(0), hi(idiv(imul(t, t), num(2))))313314# CHECK THE FREQUENCY SEPARATION AT BASE FIVE315316def g_five_small(i):317 a = iversin(imul(PI, num(Fraction(2, 5 ** i))))318 b = iversin(imul(PI, num(Fraction(4, 5 ** i))))319 return isub(num(4), imul(num(2), iadd(a, b)))320321def g_five_grid(n, m):322 return imul(num(2), iadd(icos_grid(n, m), icos_grid(2 * n, m)))323324def c_factor(j):325 a = imul(num(5), PI2)326 x = idiv(a, num(24 * 25 ** j))327 assert x[1] < 1, "the C_j product bound diverges at j=%d" % j328 return idiv(ONE, isub(ONE, x))329330def check_sep5():331 g2 = g_five_small(2)332 g3 = g_five_small(3)333 assert lo(g2) >= Fraction(36897796, 10 ** 7), "ghat_2 fell below its stated bound"334 c2 = c_factor(2)335 c3 = c_factor(3)336 assert hi(c2) <= Fraction(10033008, 10 ** 7), "C_2 above its stated bound"337 assert hi(c3) <= Fraction(10001317, 10 ** 7), "C_3 above its stated bound"338 worst = ZERO339 count = 0340 for n in range(1, 25):341 if n % 5 == 0 or n % 25 in (1, 24):342 continue343 count += 1344 v = g_five_grid(n, 25)345 m = (amax(v), amax(v))346 if hi(m) > hi(worst):347 worst = m348 assert count == 18, "the level-two enumeration lost a residue"349 assert hi(worst) <= Fraction(28242670, 10 ** 7), "the level-two maximum moved"350 ratio = idiv(worst, g2)351 assert hi(ratio) <= Fraction(7654303, 10 ** 7), "the j=2 ratio moved"352 branch2 = imul(ratio, c2)353 assert hi(branch2) <= Fraction(7679580, 10 ** 7), "the j=2 branch moved"354 branches = [branch2]355 for j in (3, 4):356 top = iadd(ONE, idiv(imul(num(12), PI), num(5 ** j)))357 span = imul(PI, num(Fraction(2, 5 ** j)))358 bot = isub(num(4), imul(num(5), imul(span, span)))359 assert bot[0] > 0, "the g_5 quadratic bound went nonpositive at j=%d" % j360 b = imul(idiv(top, bot), c_factor(j))361 target = Fraction(3264722, 10 ** 7) if j == 3 else Fraction(2651481, 10 ** 7)362 assert hi(b) <= target, "the j=%d branch moved" % j363 branches.append(b)364 best = max(hi(b) for b in branches)365 assert best <= Fraction(7679580, 10 ** 7), "the separation maximum moved"366 assert best < Fraction(768, 1000), "the separation constant r = 0.768 is no longer valid"367 return "lem:sep5, three branches and 18 level-two residues", g2, g3368369# CHECK THE EXIT AND SUBTREE CONSTANTS AT BASE THREE370371def ghat3(i):372 return imul(num(2), icos_grid(1, 3 ** i))373374def arm3(j):375 return imul(num(2), icos_grid(3 ** (j - 1) - 1, 2 * 3 ** j))376377def check_exit3():378 exit_open = idiv(icos_grid(7, 54), icos_grid(2, 54))379 assert hi(exit_open) <= Fraction(7052518, 10 ** 7), "the 0-track exit constant moved"380 exit_half = idiv(icos_grid(8, 54), icos_grid(2, 54))381 assert hi(exit_half) <= Fraction(6137010, 10 ** 7), "the pi-track exit constant moved"382 stray = idiv(icos_grid(2, 9), icos_grid(1, 9))383 assert hi(stray) <= Fraction(2266816, 10 ** 7), "the n = 2, 7 mod 9 constant moved"384 prefix = ONE385 entries = {}386 for j in range(3, 7):387 i = j - 1388 prefix = imul(prefix, idiv(icos_grid(1, 2 * 3 ** i), icos_grid(1, 3 ** i)))389 entries[j] = imul(prefix, idiv(arm3(j), ghat3(j)))390 assert hi(entries[3]) <= Fraction(7528157, 10 ** 7), "the j=3 entry constant moved"391 assert hi(entries[4]) <= Fraction(66966, 10 ** 5), "the j=4 entry constant moved"392 tail_entry = imul(iexp(LOGR), idiv(arm3(5), ghat3(5)))393 assert hi(tail_entry) <= Fraction(6419, 10 ** 4), "the j>=5 entry constant moved"394 sub = idiv(SQRT3, ghat3(3))395 assert hi(sub) <= Fraction(8900159, 10 ** 7), "the subtree constant moved"396 return "lem:exit3 four constants and lem:subtree3"397398# CHECK THE CROSSING TAILS399400def check_cross3():401 worst_lt = Fraction(0)402 worst_st = Fraction(0)403 for K in range(1, 8):404 lt = imul(_tail_ell(K + 1), num(9 ** K))405 assert hi(lt) < Fraction(19, 10), "the LT tail exceeds 1.9 at K=%d" % K406 worst_lt = max(worst_lt, hi(lt))407 for D in (38, 400):408 st = imul(imul(_tail_s(D, K + 1), num(9 ** K)), num(D - 2))409 assert hi(st) < 7, "the ST tail exceeds 7 at K=%d D=%d" % (K, D)410 worst_st = max(worst_st, hi(st))411 return "lem:cross3 tails, LT <= %.4f and ST <= %.4f" % (float(worst_lt), float(worst_st))412413# CHECK THE BASE FIVE CRITERION414415def check_crit5(g2, g3):416 r = num(Fraction(768, 1000))417 rows = []418 for D in range(18, 33, 2):419 K = 2420 kappa = imul(idiv(num(D + 4), iadd(num(D), g2)),421 ipow(idiv(num(D + 4), iadd(num(D), g3)), K - 1))422 left = imul(imul(num(2 * 5 ** K - 1), ipow(r, D - 1)), kappa)423 right = icos_grid(D - 2, 2 * 5 ** (K + 1))424 assert below(left, right), "the base-5 criterion fails at D=%d" % D425 assert 5 ** (K + 1) >= 4 * (D - 2), "the base-5 window condition fails at D=%d" % D426 rows.append((D, float(hi(left)), float(lo(right))))427 D = 16428 kappa = imul(idiv(num(D + 4), iadd(num(D), g2)), idiv(num(D + 4), iadd(num(D), g3)))429 left = imul(imul(num(2 * 5 ** 2 - 1), ipow(r, D - 1)), kappa)430 right = icos_grid(D - 2, 2 * 5 ** 3)431 assert lo(left) > hi(right), "the D=16 criterion no longer fails, so the case split moved"432 log5 = ilog(num(5))433 def log5_of(x):434 return idiv(ilog(num(x)), log5)435 b34 = imul(iadd(ONE, idiv(isub(num(4), g2), iadd(num(34), g2))),436 iexp(idiv(imul(isub(num(4), g3), isub(log5_of(4 * 34), ONE)),437 iadd(num(34), g3))))438 assert hi(b34) <= Fraction(10089188, 10 ** 7), "B(34) moved"439 thresh = idiv(imul(SQRT2, num(Fraction(1, 2))), imul(num(8), b34))440 assert lo(thresh) >= Fraction(876069, 10 ** 7), "the D>=34 threshold moved"441 h34 = imul(num(34), ipow(r, 33))442 assert hi(h34) <= Fraction(56028, 10 ** 7), "h(34) moved"443 assert below(h34, thresh), "the D>=34 tail bound fails at D=34"444 slope = isub(idiv(iadd(num(34), g3), imul(num(34), ilog(num(5)))), isub(log5_of(136), ONE))445 assert hi(slope) < 0, "the B monotonicity slope is no longer negative at D=34"446 assert lo(idiv(ONE, ineg(ilog(r)))) >= Fraction(37883, 10 ** 4), "1/|ln r| moved"447 return "prop:crit5, eight rows at even D = 18..32 and the D >= 34 chain", rows448449# CHECK THE MONOTONE DOMINATION BEYOND FOUR HUNDRED450451def check_domination(k400):452 growth = Fraction(402 * 401, 400 * 399)453 assert growth <= Fraction(10104, 10 ** 4), "the K1 growth ratio moved"454 step = imul(imul(num(growth), ipow(C_PI, 2)), num(Fraction(1001, 1000)))455 assert hi(step) <= Fraction(573, 1000), "the envelope per-step ratio moved"456 margin_step = imul(num(Fraction(402, 404)), isub(ONE, num(Fraction(1, 3 ** 12))))457 assert lo(margin_step) >= Fraction(9950, 10 ** 4), "the margin per-step ratio moved"458 assert k400 >= 12, "the 3^-K1 bound assumed a deeper certificate at D=400"459 return "thm:base3 tail, growth 1.0104, envelope step 0.573, margin step 0.9950"460461# THE STAR SCAN462463def star_scan():464 worst_slack = None465 worst_slack_at = None466 worst_width = 0.0467 worst_width_at = None468 rows = {}469 k400 = None470 for D in range(38, 401, 2):471 K = k_one(D)472 assert K - 1 >= hi(k_star(D)), "K1 fails K1 >= K* + 1 at D=%d" % D473 E = envelope(D, K)474 margin = isub(num(Fraction(4, D + 2)), delta_of(D, K))475 assert below(E, margin), "(star) fails at D=%d" % D476 slack = idiv(margin, E)477 pair = isub(ONE, ipow(num(Fraction(D - 2, D + 2)), 2))478 assert below(imul(num(2), E), pair), "the transient window bound fails at D=%d" % D479 w = max(relwidth(E), relwidth(margin))480 assert w < 1e-40, "interval width is no longer negligible at D=%d" % D481 if D == 38:482 assert lo(slack) > Fraction(111, 100), "the D=38 slack fell below 1.11"483 assert w < float(lo(slack)) * 1e-40, "the D=38 width is not far below its slack"484 if worst_slack is None or lo(slack) < worst_slack:485 worst_slack = lo(slack)486 worst_slack_at = D487 if w > worst_width:488 worst_width = w489 worst_width_at = D490 rows[D] = (K, float(lo(slack)))491 if D == 400:492 k400 = K493 return worst_slack, worst_slack_at, worst_width, worst_width_at, rows, k400494495# RUN496497def run():498 results = []499 t = time.time()500 dom, g2, g3 = check_sep5()501 results.append(("certify_sep5", dom, time.time() - t))502 t = time.time()503 dom = check_exit3()504 results.append(("certify_exit3", dom, time.time() - t))505 t = time.time()506 dom = check_cross3()507 results.append(("certify_cross3", dom, time.time() - t))508 t = time.time()509 dom, rows5 = check_crit5(g2, g3)510 results.append(("certify_crit5", dom, time.time() - t))511 t = time.time()512 slack, slack_at, width, width_at, rows, k400 = star_scan()513 results.append(("certify_star", "182 inequalities, even D = 38..400", time.time() - t))514 t = time.time()515 dom = check_domination(k400)516 results.append(("certify_domination", dom, time.time() - t))517 summary = ("lem:sep5, lem:exit3, lem:subtree3, lem:cross3, prop:crit5, and the 182 "518 "inequalities of (star) at even D = 38..400; worst slack %.4f at D = %d, "519 "worst relative interval width %.1e" % (float(slack), slack_at, width))520 return summary, results, rows, rows5521522# MAIN523524def main():525 t0 = time.time()526 summary, results, rows, rows5 = run()527 for name, dom, secs in results:528 print("%s: PASS (%s, %.1f s)" % (name, dom, secs))529 print("total: %.1f s" % (time.time() - t0))530 print(summary)531 print("base 5 criterion, even D = 18..32")532 print(" D left right")533 for D, left, right in rows5:534 print(" %-4d %-9.4f %.4f" % (D, left, right))535 print("base 3 certificate depth and (star) slack, even D = 38..400")536 print(" D K1 slack")537 for D in (38, 40, 42, 50, 100, 200, 300, 400):538 K, s = rows[D]539 print(" %-4d %-6d %.4g" % (D, K, s))540541if __name__ == "__main__":542 main()