verify.py
24.3 kB · python · 831 lines
1import math2import os3import sys4import time5from fractions import Fraction6from itertools import product7from math import comb89sys.dont_write_bytecode = True10sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))11import certify1213# POLYNOMIALS1415def polymul(a, b):16 out = [0] * (len(a) + len(b) - 1)17 for i, x in enumerate(a):18 if x:19 for j, y in enumerate(b):20 if y:21 out[i + j] += x * y22 return out2324def polypow(a, e):25 r = [1]26 for _ in range(e):27 r = polymul(r, a)28 return r2930def middle(q):31 assert q % 2 == 1 and q >= 3, "base must be odd and at least 3"32 return (q - 1) // 23334def alpha_poly(q):35 m = middle(q)36 return [0 if d == m else 1 for d in range(q)]3738def digit_poly(D, q):39 m = middle(q)40 A = alpha_poly(q)41 second = A[:]42 second[m] += D43 return polymul(polypow(A, D - 1), second)4445def digit_poly_base3_binomial(D):46 base = [0] * (2 * D - 1)47 for k in range(D):48 base[2 * k] = comb(D - 1, k)49 return polymul(base, [1, D, 1])5051def fill_value(D, q):52 return (q - 1) ** (D - 1) * (q - 1 + D)5354def brute_poly(D, q):55 m = middle(q)56 out = [0] * ((q - 1) * D + 1)57 for v in product(range(q), repeat=D):58 if sum(1 for x in v if x == m) <= 1:59 out[sum(v)] += 160 return out6162# CARRY MATRIX6364def half_width(D):65 return (D - 1) // 26667def coefficient(P, x):68 return P[x] if 0 <= x < len(P) else 06970def assert_window_closed(D, q):71 P = digit_poly(D, q)72 shift = middle(q) * D73 H = half_width(D)74 for c in range(-H, H + 1):75 for s, val in enumerate(P):76 if not val:77 continue78 if (c + shift - s) % q:79 continue80 cp = (c + shift - s) // q81 assert abs(cp) <= H, "carry window escapes at D=%d q=%d" % (D, q)8283def m_full(D, q):84 P = digit_poly(D, q)85 shift = middle(q) * D86 H = half_width(D)87 st = list(range(-H, H + 1))88 return [[coefficient(P, c + shift - q * cp) for c in st] for cp in st]8990def m_even(D, q):91 P = digit_poly(D, q)92 shift = middle(q) * D93 H = half_width(D)94 n = H + 195 N = [[0] * n for _ in range(n)]96 for cp in range(n):97 for c in range(n):98 v = coefficient(P, c + shift - q * cp)99 if c:100 v += coefficient(P, -c + shift - q * cp)101 N[cp][c] = v102 return N103104def m_odd(D, q):105 P = digit_poly(D, q)106 shift = middle(q) * D107 H = half_width(D)108 return [[coefficient(P, c + shift - q * cp) - coefficient(P, -c + shift - q * cp)109 for c in range(1, H + 1)] for cp in range(1, H + 1)]110111# CENSUS BY CARRY TRANSFER112113def census(D, q, top):114 P = digit_poly(D, q)115 shift = middle(q) * D116 H = half_width(D)117 n = 2 * H + 1118 step = []119 for cp in range(-H, H + 1):120 row = []121 for c in range(-H, H + 1):122 v = coefficient(P, c + shift - q * cp)123 if v:124 row.append((c + H, v))125 step.append(row)126 u = [0] * n127 u[H] = 1128 bs = [1]129 m0s = [1]130 for _ in range(top):131 v = [0] * n132 for i in range(n):133 s = 0134 for j, w in step[i]:135 if u[j]:136 s += w * u[j]137 v[i] = s138 u = v139 bs.append(sum(u))140 m0s.append(sum(u[c + H] for c in range(-H, H + 1) if c % q == 0))141 return bs, m0s142143# EXACT DETERMINANT144145def bareiss_det(A):146 n = len(A)147 if n == 0:148 return 1149 M = [row[:] for row in A]150 sign = 1151 prev = 1152 for k in range(n - 1):153 if M[k][k] == 0:154 p = None155 for r in range(k + 1, n):156 if M[r][k] != 0:157 p = r158 break159 if p is None:160 return 0161 M[k], M[p] = M[p], M[k]162 sign = -sign163 for i in range(k + 1, n):164 for j in range(k + 1, n):165 M[i][j] = (M[i][j] * M[k][k] - M[i][k] * M[k][j]) // prev166 prev = M[k][k]167 return sign * M[n - 1][n - 1]168169def fraction_det(A):170 n = len(A)171 if n == 0:172 return Fraction(1)173 M = [[Fraction(x) for x in row] for row in A]174 det = Fraction(1)175 for k in range(n):176 p = None177 for r in range(k, n):178 if M[r][k] != 0:179 p = r180 break181 if p is None:182 return Fraction(0)183 if p != k:184 M[k], M[p] = M[p], M[k]185 det = -det186 det *= M[k][k]187 inv = 1 / M[k][k]188 for j in range(k, n):189 M[k][j] *= inv190 for i in range(k + 1, n):191 f = M[i][k]192 if f:193 for j in range(k, n):194 M[i][j] -= f * M[k][j]195 return det196197# COLLATZ-WIELANDT CERTIFICATES198199def transpose_action(D, q):200 P = digit_poly(D, q)201 shift = middle(q) * D202 H = half_width(D)203 n = 2 * H + 1204 rows = []205 for c in range(-H, H + 1):206 row = []207 for cp in range(-H, H + 1):208 v = coefficient(P, c + shift - q * cp)209 if v:210 row.append((cp + H, v))211 rows.append(row)212 return rows, n213214def certificate_depth(D, q, kmax):215 rows, n = transpose_action(D, q)216 f = fill_value(D, q)217 x = [1] * n218 for K in range(kmax + 1):219 y = [sum(w * x[j] for j, w in row) for row in rows]220 assert all(v > 0 for v in y), "beta vector not positive at D=%d q=%d" % (D, q)221 if all(q * y[i] < f * x[i] for i in range(n)):222 return K, x, y223 x = y224 return None, None, None225226# GF(2) LINEAR ALGEBRA ON BITMASKS227228def bit_rank(rows):229 piv = []230 rk = 0231 for v in rows:232 for p in piv:233 hb = p.bit_length() - 1234 if (v >> hb) & 1:235 v ^= p236 if v:237 piv.append(v)238 piv.sort(key=lambda t: -t.bit_length())239 rk += 1240 return rk241242def echelon(rows):243 B = []244 for v in rows:245 for p in B:246 hb = p.bit_length() - 1247 if (v >> hb) & 1:248 v ^= p249 if v:250 B.append(v)251 B.sort(key=lambda t: -t.bit_length())252 return B253254def reduce_by(v, B):255 for p in B:256 hb = p.bit_length() - 1257 if (v >> hb) & 1:258 v ^= p259 return v260261def nullspace(rows, ncols):262 M = [r for r in rows if r]263 piv = {}264 for c in range(ncols):265 bit = 1 << c266 pr = None267 for r in M:268 if r & bit:269 pr = r270 break271 if pr is None:272 continue273 M = [r for r in M if r is not pr]274 M = [(r ^ pr) if (r & bit) else r for r in M]275 for c2 in list(piv):276 if piv[c2] & bit:277 piv[c2] ^= pr278 piv[c] = pr279 out = []280 for f in [c for c in range(ncols) if c not in piv]:281 v = 1 << f282 for c2, pr in piv.items():283 if (pr >> f) & 1:284 v |= 1 << c2285 out.append(v)286 return out287288# MOD-2 EVEN CORE BY LUCAS289290def m_even_rows_mod2(D):291 assert D % 2 == 1, "the Lucas model is stated for odd D"292 mask = 2 * D - 3293294 def pc(x):295 a = 1 if 0 <= x <= mask and (x & ~mask) == 0 else 0296 b = 1 if 0 <= x - 3 <= mask and ((x - 3) & ~mask) == 0 else 0297 return a ^ b298299 n = (D - 1) // 2 + 1300 rows = []301 for cp in range(n):302 v = 0303 for c in range(n):304 bit = pc(c + D - 3 * cp)305 if c:306 bit ^= pc(-c + D - 3 * cp)307 if bit:308 v |= 1 << c309 rows.append(v)310 return rows, n311312# JACOBSTHAL AND THE TENT313314def jacobsthal(k):315 return (2 ** k - (-1) ** k) // 3316317def trough_set(top):318 T = set()319 k = 0320 while True:321 j = jacobsthal(k)322 T.add(2 * j + 1)323 T.add(2 * j + 3)324 if 2 * j + 1 > 2 * top + 8:325 break326 k += 1327 return sorted(T)328329def tent(D, T):330 return min(abs(D - t) // 2 + 1 for t in T)331332# THE LAYER-2 CLOSED FORM333334def law_e(R):335 b = 0336 while (1 << b) < 3 * R - 1:337 b += 1338 g = abs(2 * R - (1 << (b - 1)) - 1)339 e = 1340 while jacobsthal(e) < (g + 1) // 2:341 e += 1342 k = b - 1 - e343 assert k >= 1, "Law E slot with k < 1 at R=%d" % R344 return b, g, e, k345346def window_top(R, b):347 A = 1 << b348 i0 = max(2, 4 * R - A)349 hi = (6 * R + 2 - A) // 3350 hi_even = hi - (hi % 2)351 assert (hi_even - i0) % 2 == 0, "window parity at R=%d" % R352 return (hi_even - i0) // 2353354def layer2_closed_form(D):355 R = (D - 1) // 2356 b, g, e, k = law_e(R)357 K = window_top(R, b)358 N = jacobsthal(k)359 t = (N - 1) // 2360 deficit = 2 * jacobsthal(e - 1) if k % 2 == 0 else 0361 C = K - deficit362 B = t * (1 << e)363 w = C - B364 ce = 2 * jacobsthal(e - 2) - 1 if e >= 3 else 1365 return K + 1, min(w + 1, ce + 1 - w)366367# THE MOD-4 LIFT368369def even_core_mod4(D):370 P = digit_poly_base3_binomial(D)371 n = (D - 1) // 2 + 1372 rows = []373 for cp in range(n):374 lo = 0375 hi = 0376 for c in range(n):377 v = coefficient(P, c + D - 3 * cp)378 if c:379 v += coefficient(P, -c + D - 3 * cp)380 v %= 4381 if v & 1:382 lo |= 1 << c383 if v & 2:384 hi |= 1 << c385 rows.append((lo, hi))386 return rows, n387388def layers_12(D):389 rows, n = even_core_mod4(D)390 mod2 = [lo for lo, _ in rows]391 V1 = nullspace(mod2, n)392 cols = [0] * n393 for cp, r in enumerate(mod2):394 rr = r395 while rr:396 c = (rr & -rr).bit_length() - 1397 cols[c] |= 1 << cp398 rr &= rr - 1399 image = echelon(cols)400 obstruction = []401 for v in V1:402 o = 0403 for cp in range(n):404 lo, hi = rows[cp]405 s = (bin(lo & v).count("1") + 2 * bin(hi & v).count("1")) % 4406 assert s % 2 == 0, "kernel vector not annihilated mod 2 at D=%d" % D407 if s == 2:408 o |= 1 << cp409 obstruction.append(reduce_by(o, image))410 L = len(V1)411 orows = []412 for cp in range(n):413 r = 0414 for j in range(L):415 if (obstruction[j] >> cp) & 1:416 r |= 1 << j417 if r:418 orows.append(r)419 combos = nullspace(orows, L)420 V2 = []421 for cb in combos:422 v = 0423 for j in range(L):424 if (cb >> j) & 1:425 v ^= V1[j]426 if v:427 V2.append(v)428 return L, len(echelon(V2))429430# CHECK POLYNOMIALS431432def polyphase(P):433 parts = [[], [], []]434 for i, x in enumerate(P):435 parts[i % 3].append(x)436 return parts437438def extraction_image(P, X):439 prod = polymul(P, X)440 top = (len(prod) - 2) // 3 + 1441 return [coefficient(prod, 3 * j + 1) for j in range(max(top, 0))]442443def check_polynomials():444 for q in (3, 5):445 for D in range(2, 9):446 P = digit_poly(D, q)447 B = brute_poly(D, q)448 while len(B) < len(P):449 B.append(0)450 assert P == B[:len(P)] and all(x == 0 for x in B[len(P):]), \451 "digit polynomial disagrees with brute force at q=%d D=%d" % (q, D)452 assert sum(P) == fill_value(D, q), "fill wrong at q=%d D=%d" % (q, D)453 if q == 3:454 assert P == digit_poly_base3_binomial(D), \455 "base-3 factorization wrong at D=%d" % D456 seed = 12345457 for D in range(3, 14, 2):458 R = (D - 1) // 2459 P = digit_poly(D, 3)460 A = m_full(D, 3)461 for c in range(-R, R + 1):462 X = [0] * (2 * R + 1)463 X[R - c] = 1464 img = extraction_image(P, X)465 for cp in range(-R, R + 1):466 assert A[cp + R][c + R] == img[R - cp], \467 "extraction form fails at D=%d c=%d c'=%d" % (D, c, cp)468 P0, P1, P2 = polyphase(P)469 for _ in range(6):470 X = []471 for _ in range(2 * R + 1):472 seed = (seed * 1103515245 + 12345) % (1 << 31)473 X.append(seed % 7 - 3)474 X0, X1, X2 = polyphase(X)475 lhs = extraction_image(P, X)476 parts = (polymul(P1, X0), polymul(P0, X1), [0] + polymul(P2, X2))477 width = max(len(lhs), max(len(p) for p in parts))478 rhs = [0] * width479 for part in parts:480 for j, x in enumerate(part):481 rhs[j] += x482 lhs = lhs + [0] * (width - len(lhs))483 assert lhs == rhs, "polyphase form fails at D=%d" % D484 return "bases 3 and 5, D = 2..8; extraction and polyphase forms, odd D = 3..13"485486# CHECK REDUCTION487488def check_reduction():489 kmax = 6490 for q in (3, 5, 7):491 for D in range(2, 13):492 assert_window_closed(D, q)493 bs, m0s = census(D, q, kmax)494 f = fill_value(D, q)495 sgn = (-1) ** (D - 1)496 for k in range(1, kmax + 1):497 W = q ** (k - 1) * (q * bs[k] - f * bs[k - 1])498 pred = sgn * (D - 1) * q ** (k - 1) * (q * m0s[k - 1] - bs[k - 1])499 assert W == pred, "step identity fails at q=%d D=%d k=%d" % (q, D, k)500 return "bases 3, 5, 7, D = 2..12, k = 1..6"501502# CHECK FOLD503504def check_fold():505 for q in (3, 5):506 for D in range(2, 26):507 full = bareiss_det(m_full(D, q))508 even = bareiss_det(m_even(D, q))509 odd = bareiss_det(m_odd(D, q))510 assert full == even * odd, "fold fails at q=%d D=%d" % (q, D)511 if D <= 9:512 assert Fraction(full) == fraction_det(m_full(D, q)), \513 "Bareiss disagrees with fractions at q=%d D=%d" % (q, D)514 return "bases 3 and 5, D = 2..25, both parities"515516# CHECK EVEN CERTIFICATES517518def check_even_certificates():519 depths = {}520 plan = [521 (5, list(range(2, 41, 2)) + [64, 66, 314, 316]),522 (3, list(range(2, 37, 2))),523 (7, list(range(2, 27, 2)) + [172, 174]),524 (9, list(range(2, 43, 2))),525 (11, list(range(2, 61, 2))),526 ]527 for q, dims in plan:528 for D in dims:529 K, x, y = certificate_depth(D, q, 4 * D + 8)530 assert K is not None, "no certificate found at q=%d D=%d" % (q, D)531 f = fill_value(D, q)532 assert all(q * y[i] < f * x[i] for i in range(len(x))), \533 "certificate does not close at q=%d D=%d" % (q, D)534 depths[(q, D)] = K535 assert depths[(3, 30)] == 50, "base-3 depth at D=30 changed"536 assert depths[(3, 36)] == 72, "base-3 depth at D=36 changed"537 for q in (5, 7, 9, 11):538 assert depths[(q, 2)] == 0, "depth at D=2 changed at q=%d" % q539 assert depths[(q, 4)] == 1, "the K=0 death is not at D=4 at q=%d" % q540 for D in range(4, 15, 2):541 assert depths[(5, D)] == 1, "base-5 staircase moved at D=%d" % D542 for D in list(range(16, 41, 2)) + [64]:543 assert depths[(5, D)] == 2, "base-5 staircase moved at D=%d" % D544 assert depths[(5, 66)] == 3 and depths[(5, 314)] == 3 and depths[(5, 316)] == 4, \545 "a base-5 death dimension moved"546 for D in range(4, 25, 2):547 assert depths[(7, D)] == 1, "base-7 staircase moved at D=%d" % D548 assert depths[(7, 26)] == 2 and depths[(7, 172)] == 2 and depths[(7, 174)] == 3, \549 "a base-7 death dimension moved"550 for D in range(4, 41, 2):551 assert depths[(9, D)] == 1, "base-9 staircase moved at D=%d" % D552 assert depths[(9, 42)] == 2, "the base-9 K=1 death is not at D=42"553 for D in range(4, 59, 2):554 assert depths[(11, D)] == 1, "base-11 staircase moved at D=%d" % D555 assert depths[(11, 60)] == 2, "the base-11 K=1 death is not at D=60"556 dom = ("base 5 even D <= 40 plus 64, 66, 314, 316; base 3 even D <= 36; "557 "base 7 even D <= 26 plus 172, 174; base 9 even D <= 42; base 11 even D <= 60")558 return dom, depths559560# CHECK NO TRANSIENT561562def check_notransient():563 for q in (7, 9, 11):564 for D in range(2, 25, 2):565 assert_window_closed(D, q)566 bs, m0s = census(D, q, 12)567 for L in range(13):568 assert q * m0s[L] - bs[L] > 0, \569 "V dips at q=%d D=%d L=%d" % (q, D, L)570 return "bases 7, 9, 11, even D <= 24, L <= 12, window closure included"571572# THE CROSSING LEVEL573574def k_star(D):575 levels = range(2, 60)576 log_r = sum(math.log(math.cos(math.pi / 3 ** i) / math.cos(2 * math.pi / 3 ** i)) for i in levels)577 s = sum(math.log((D - 2 * math.cos(math.pi / 3 ** i)) / (D - 2))578 - math.log((D + 2 * math.cos(2 * math.pi / 3 ** i)) / (D + 2)) for i in levels)579 return ((D - 1) * log_r + s) / math.log((D + 2) / (D - 2))580581def l_zero(D):582 ks = k_star(D)583 assert abs(ks - round(ks)) > 1e-6, "K* within 1e-6 of an integer at D=%d" % D584 L0 = math.floor(ks)585 return L0 if L0 % 2 else L0 - 1586587# THE ROW SWEEP588589def row_sweep(D, cap):590 P = digit_poly(D, 3)591 assert all(v >= 0 for v in P), "the digit polynomial has a negative coefficient at D=%d" % D592 assert all(P[j] == P[2 * D - j] for j in range(2 * D + 1)), \593 "the digit polynomial is not palindromic at D=%d" % D594 N = m_even(D, 3)595 n = len(N)596 cols = [[(cp, N[cp][c]) for cp in range(n) if N[cp][c]] for c in range(n)]597 r = [(3 if c % 3 == 0 else 0) - 1 for c in range(n)]598 r = [r[c] if c == 0 else 2 * r[c] for c in range(n)]599 last_neg = -1600 for k in range(1, cap + 1):601 r = [sum(w * r[cp] for cp, w in col) for col in cols]602 if r[0] < 0:603 last_neg = k604 if min(r) >= 0:605 assert min(r) > 0, "the row certificate closes without strictness at D=%d" % D606 return k, last_neg607 raise AssertionError("no nonnegative row within the cap at D=%d" % D)608609# CHECK TRANSIENT610611def check_transient(depths):612 table = {}613 for D in range(6, 75, 2):614 assert_window_closed(D, 3)615 K, star = row_sweep(D, D * D // 8 + 200)616 assert star >= 1, "no transient at D=%d" % D617 assert star == l_zero(D), \618 "Lstar is not the greatest odd integer below K* at D=%d: %d vs %d" % (D, star, l_zero(D))619 assert K in (star + 1, star + 2), \620 "certificate depth off the transient at D=%d: K=%d Lstar=%d" % (D, K, star)621 if D <= 36:622 assert K == depths[(3, D)], \623 "the row and column certificates disagree at D=%d: %d vs %d" % (D, K, depths[(3, D)])624 table[D] = (star, K)625 plus_two = sorted(D for D in table if table[D][1] == table[D][0] + 2)626 assert plus_two == [14, 22, 32, 38, 40, 48, 52, 54, 58, 70, 72], \627 "the K_min = Lstar + 2 rows moved: %s" % plus_two628 return "base 3, even D = 6..74, exhaustion by the row certificate", table629630# LEMMA M FINITE WINDOWS631632def valuation_1t(mask):633 v = 0634 while mask and bin(mask).count("1") % 2 == 0:635 n = mask.bit_length()636 x = mask637 s = 1638 while s < n:639 x ^= x << s640 s <<= 1641 mask = x & ((1 << n) - 1)642 v += 1643 return v644645def brute_V(r, d):646 allowed = [e for e in range(d + 1) if e % 3 != r]647 best = -1648 for bits in range(1, 1 << len(allowed)):649 mask = 0650 b = bits651 i = 0652 while b:653 if b & 1:654 mask |= 1 << allowed[i]655 b >>= 1656 i += 1657 v = valuation_1t(mask)658 if v > best:659 best = v660 return best661662def mu_formula(r, d):663 best = -1664 b = 0665 while (1 << b) <= d:666 s0 = (r + (1 << b)) % 3667 if (1 << b) + s0 <= d:668 h = (d - s0 - (1 << b)) // 3 + (1 << b)669 if h > best:670 best = h671 b += 1672 return best673674# CHECK TENT675676def check_tent(deep):677 top = 583 if deep else 301678 T = trough_set(top)679 peaks = {3} | {(1 << (2 * j)) + 1 for j in range(1, 12)}680 for D in range(3, top + 1, 2):681 rows, n = m_even_rows_mod2(D)682 nullity = n - bit_rank(rows)683 pred = tent(D, T)684 assert nullity == pred, \685 "tent law fails at D=%d: nullity=%d tent=%d" % (D, nullity, pred)686 cap = -(-n // 3)687 assert pred <= cap, "tent above the cap at D=%d" % D688 assert (pred == cap) == (D in peaks), "cap equality misplaced at D=%d" % D689 for r in range(3):690 for d in range(2, 13):691 assert brute_V(r, d) == mu_formula(r, d), \692 "valuation window fails at r=%d d=%d" % (r, d)693 return "base 3, odd D = 3..%d, plus the valuation windows r = 0,1,2, d = 2..12" % top694695# CHECK LAYER 2696697def check_layer2(deep):698 top = 401 if deep else 151699 for D in range(5, top + 1, 2):700 L1, L2 = layers_12(D)701 p1, p2 = layer2_closed_form(D)702 assert L1 == p1, "layer-1 window count fails at D=%d: %d vs %d" % (D, L1, p1)703 assert L2 == p2, "layer-2 closed form fails at D=%d: %d vs %d" % (D, L2, p2)704 assert L2 >= 1, "empty layer 2 at D=%d" % D705 return "base 3, odd D = 5..%d" % top706707# CHECK STRICTNESS708709def v2_det(A, prec):710 mod = 1 << prec711 n = len(A)712 M = [[x % mod for x in row] for row in A]713 total = 0714 for k in range(n):715 best = None716 bv = prec717 for r in range(k, n):718 x = M[r][k]719 if x:720 v = (x & -x).bit_length() - 1721 if v < bv:722 best = r723 bv = v724 assert best is not None, "pivot vanished mod 2^%d" % prec725 M[k], M[best] = M[best], M[k]726 total += bv727 inv = pow(M[k][k] >> bv, -1, mod)728 for r in range(k + 1, n):729 x = M[r][k]730 if x:731 f = ((x >> bv) * inv) % mod732 for j in range(k, n):733 M[r][j] = (M[r][j] - f * M[k][j]) % mod734 assert total + 128 < prec, "valuation too close to the working precision"735 return total736737def check_strictness():738 prec = 1024739 E = m_even(7, 3)740 f = fill_value(7, 3)741 pencil = [[f * (i == j) - 3 * E[i][j] for j in range(4)] for i in range(4)]742 delta7 = bareiss_det(pencil)743 assert delta7 == -148506048, "the D=7 pencil determinant changed"744 assert (delta7 & -delta7).bit_length() - 1 == 6, "v2 of the D=7 pencil changed"745 assert delta7 == -(1 << 6) * 2320407, "the odd part of the D=7 pencil changed"746 assert v2_det(E, prec) == 7, "v2 at D=7 changed"747 tight = []748 for D in sorted(set(range(9, 116, 2)) | set(range(13, 152, 6))):749 n = (D + 1) // 2750 val = v2_det(m_even(D, 3), prec)751 assert val <= n, "v2 above n at D=%d: v2=%d n=%d" % (D, val, n)752 if val == n:753 tight.append(D)754 if D % 3 == 1:755 assert val < D - 1, "base-3 strictness fails at D=%d: v2=%d" % (D, val)756 assert tight == [9, 15], "the v2 = n equality set moved: %s" % tight757 for D in range(6, 157, 5):758 threshold = 2 * (D - 1) + (((D + 4) & -(D + 4)).bit_length() - 1)759 val = v2_det(m_even(D, 5), prec)760 assert val < threshold, "base-5 strictness fails at D=%d: v2=%d" % (D, val)761 dom = ("base 3 odd D = 9..115 and class members D = 13..151 plus D = 7 direct, "762 "base 5 class members D = 6..156")763 return dom, delta7764765# CHECK THE INTERVAL CERTIFICATES766767def check_intervals():768 summary, _results, _rows, _rows5 = certify.run()769 return summary770771# MAIN772773def main():774 deep = "--deep" in sys.argv775 results = []776 t0 = time.time()777778 t = time.time()779 dom = check_polynomials()780 results.append(("check_polynomials", dom, time.time() - t))781782 t = time.time()783 dom = check_reduction()784 results.append(("check_reduction", dom, time.time() - t))785786 t = time.time()787 dom = check_fold()788 results.append(("check_fold", dom, time.time() - t))789790 t = time.time()791 dom, depths = check_even_certificates()792 results.append(("check_even_certificates", dom, time.time() - t))793794 t = time.time()795 dom, table = check_transient(depths)796 results.append(("check_transient", dom, time.time() - t))797798 t = time.time()799 dom = check_notransient()800 results.append(("check_notransient", dom, time.time() - t))801802 t = time.time()803 dom = check_tent(deep)804 results.append(("check_tent", dom, time.time() - t))805806 t = time.time()807 dom = check_layer2(deep)808 results.append(("check_layer2", dom, time.time() - t))809810 t = time.time()811 dom, delta7 = check_strictness()812 results.append(("check_strictness", dom, time.time() - t))813814 t = time.time()815 dom = check_intervals()816 results.append(("check_intervals", dom, time.time() - t))817818 for name, dom, secs in results:819 print("%s: PASS (%s, %.1f s)" % (name, dom, secs))820 print("total: %.1f s" % (time.time() - t0))821 print("base 3 transient and certificate depth, even D = 6..74")822 dims = sorted(table)823 for i in range(0, len(dims), 5):824 row = ""825 for D in dims[i:i + 5]:826 star, K = table[D]827 row += " %3d: %4d/%4d" % (D, star, K)828 print(row)829830if __name__ == "__main__":831 main()