verify.py
19.2 kB · python · 539 lines
1#!/usr/bin/env python32# CARPETS34import math5from fractions import Fraction678def sign(n, i, L):9 return 1 if ((n * i) // L) % 2 == 0 else -1101112def integral(m, n):13 L = m * n // math.gcd(m, n)14 total = 015 for i in range(L):16 total += sign(m, i, L) * sign(n, i, L)17 return Fraction(total, L)181920def mean(m):21 total = 022 for i in range(m):23 total += sign(m, i, m)24 return Fraction(total, m)252627def master(m, n):28 g = math.gcd(m, n)29 if (m // g) % 2 == 1 and (n // g) % 2 == 1:30 return Fraction(g * g, m * n)31 return Fraction(0)323334def cov_tree(m, n):35 d = math.gcd(m, n)36 return Fraction(d * d - 1, 4 * m * n)373839def cov_carpet(m, n):40 d = math.gcd(m, n)41 return Fraction((d * d - 1) * (2 * (m - 1) * (n - 1) + d * d - 1), 16 * m * m * n * n)424344def cov_net(m, n):45 d = math.gcd(m, n)46 return Fraction((d * d - 1) * (2 * (m + 1) * (n + 1) + d * d - 1), 16 * m * m * n * n)474849def cov_void(m, n):50 d = math.gcd(m, n)51 return Fraction(d ** 4 - 1, 4 * m * m * n * n)525354def assembled(m, n):55 A = integral(m, n)56 am = Fraction(1, m)57 an = Fraction(1, n)58 mu_m = (1 - am) / 259 mu_n = (1 - an) / 260 nu = (1 - am - an + A) / 461 bar_m = (1 + am) / 262 bar_n = (1 + an) / 263 nub = (1 + am + an + A) / 464 tree = nu - mu_m * mu_n65 carpet = nu * nu - (mu_m * mu_n) ** 266 net = nub * nub - (bar_m * bar_n) ** 267 void = (A * A - (am * an) ** 2) / 468 return tree, carpet, net, void697071def r_tree(m, n):72 d = math.gcd(m, n)73 return (d * d - 1) / math.sqrt((m * m - 1) * (n * n - 1))747576def r_carpet(m, n):77 d = math.gcd(m, n)78 top = (d * d - 1) * (2 * (m - 1) * (n - 1) + d * d - 1)79 bot = (m - 1) * (n - 1) * math.sqrt((3 * m - 1) * (m + 1) * (3 * n - 1) * (n + 1))80 return top / bot818283def r_net(m, n):84 d = math.gcd(m, n)85 top = (d * d - 1) * (2 * (m + 1) * (n + 1) + d * d - 1)86 bot = (m + 1) * (n + 1) * math.sqrt((3 * m + 1) * (m - 1) * (3 * n + 1) * (n - 1))87 return top / bot888990def r_void(m, n):91 d = math.gcd(m, n)92 return (d ** 4 - 1) / math.sqrt((m ** 4 - 1) * (n ** 4 - 1))939495def brute_carpet(m, n):96 L = m * n // math.gcd(m, n)97 chi_m = [(1 - sign(m, i, L)) // 2 for i in range(L)]98 chi_n = [(1 - sign(n, i, L)) // 2 for i in range(L)]99 both = 0100 sm = 0101 sn = 0102 for i in range(L):103 for j in range(L):104 a = chi_m[i] * chi_m[j]105 b = chi_n[i] * chi_n[j]106 both += a * b107 sm += a108 sn += b109 cells = Fraction(1, L * L)110 return both * cells - (sm * cells) * (sn * cells)111112113def brute_net(m, n):114 L = m * n // math.gcd(m, n)115 bar_m = [(1 + sign(m, i, L)) // 2 for i in range(L)]116 bar_n = [(1 + sign(n, i, L)) // 2 for i in range(L)]117 both = 0118 sm = 0119 sn = 0120 for i in range(L):121 for j in range(L):122 a = bar_m[i] * bar_m[j]123 b = bar_n[i] * bar_n[j]124 both += a * b125 sm += a126 sn += b127 cells = Fraction(1, L * L)128 return both * cells - (sm * cells) * (sn * cells)129130131def brute_void(m, n):132 L = m * n // math.gcd(m, n)133 both = 0134 sm = 0135 sn = 0136 for i in range(L):137 for j in range(L):138 a = (1 + sign(m, i, L) * sign(m, j, L)) // 2139 b = (1 + sign(n, i, L) * sign(n, j, L)) // 2140 both += a * b141 sm += a142 sn += b143 cells = Fraction(1, L * L)144 return both * cells - (sm * cells) * (sn * cells)145146147def tri(x):148 r = x % 2149 return 1 - 2 * min(r, 2 - r)150151152def tri_scaled(k, b):153 r = k % (2 * b)154 return b - 2 * min(r, 2 * b - r)155156157def saw(x):158 r = x % 1159 return r - Fraction(1, 2)160161162def franel(m, n):163 L = m * n // math.gcd(m, n)164 h = Fraction(1, L)165 total = Fraction(0)166 for i in range(L):167 b = h * i + h / 2168 total += h * saw(m * b) * saw(n * b) + Fraction(m * n) * h ** 3 / 12169 return total170171172def ray(a, b):173 if b % 2 == 0:174 return Fraction(0)175 return Fraction((-1) ** a, b * b)176177178def is_prime(n):179 if n < 2:180 return False181 k = 2182 while k * k <= n:183 if n % k == 0:184 return False185 k += 1186 return True187188189def mobius_upto(N):190 mu = [1] * (N + 1)191 primes = []192 small = [0] * (N + 1)193 for i in range(2, N + 1):194 if small[i] == 0:195 primes.append(i)196 small[i] = i197 mu[i] = -1198 for p in primes:199 if p > small[i] or i * p > N:200 break201 small[i * p] = p202 mu[i * p] = 0 if p == small[i] else -mu[i]203 return mu204205206def fourier(dim):207 keys = list(range(dim))208 single = []209 pair = []210 triple = 0211 const = 0212 patterns = []213 for mask in range(2 ** dim):214 sigma = [1 - 2 * ((mask >> t) & 1) for t in keys]215 odd = sum(1 for x in sigma if x == -1)216 patterns.append((sigma, 1 if odd <= 1 else 0))217 N = float(2 ** dim)218 for sigma, f in patterns:219 const += f220 const /= N221 for t in keys:222 c = sum(f * sigma[t] for sigma, f in patterns) / N223 single.append(c)224 for t in keys:225 for u in keys:226 if u > t:227 c = sum(f * sigma[t] * sigma[u] for sigma, f in patterns) / N228 pair.append(c)229 if dim == 3:230 triple = sum(f * sigma[0] * sigma[1] * sigma[2] for sigma, f in patterns) / N231 return const, single, pair, triple232233234235def sine_moment(n, a):236 total = 0.0237 for k in range(n):238 cell = math.cos(math.pi * a * k / n) - math.cos(math.pi * a * (k + 1) / n)239 if k % 2 == 1:240 total += cell241 return total / (math.pi * a)242243244def spectrum_direct(S, a, b):245 total = 0.0246 for n in S:247 total += sine_moment(n, a) * sine_moment(n, b)248 return total / len(S)249250251def spectrum_formula(S, a, b):252 L = len(S)253 g = 0254 s1a = sum(n for n in S if a % n == 0)255 s1b = sum(n for n in S if b % n == 0)256 gg = math.gcd(a, b)257 s2 = sum(n * n for n in S if gg % n == 0)258 return (1.0 - s1a / L - s1b / L + s2 / L) / (math.pi ** 2 * a * b)259260261def zeta_em(s, N=20000):262 total = 0.0263 for n in range(1, N + 1):264 total += n ** (-s)265 total += N ** (1 - s) / (s - 1) - 0.5 * N ** (-s) + s * N ** (-s - 1) / 12266 return total267268269def lam(s):270 return (1 - 2.0 ** (-s)) * zeta_em(s)271272def main():273 bad = []274 for m in range(1, 41):275 for n in range(m, 41):276 if integral(m, n) != master(m, n):277 bad.append((m, n))278 assert bad == [], "refined master law: got mismatches %s, want []" % bad[:4]279 naive = 0280 for m in range(1, 41):281 for n in range(m, 41):282 g = math.gcd(m, n)283 if integral(m, n) != Fraction(g * g, m * n):284 naive += 1285 assert naive == 532, "naive law failures to 40: got %d, want 532" % naive286 assert integral(1, 2) == 0, "integral(1,2): got %s, want 0" % integral(1, 2)287 assert integral(2, 6) == Fraction(1, 3), "integral(2,6): got %s, want 1/3" % integral(2, 6)288 print("master law, all pairs 1..40: 0 mismatches, 532 naive failures")289290 odds = list(range(1, 100, 2))291 misses = 0292 for a in range(len(odds)):293 for b in range(a, len(odds)):294 m, n = odds[a], odds[b]295 g = math.gcd(m, n)296 if integral(m, n) != Fraction(g * g, m * n):297 misses += 1298 assert misses == 0, "odd master law to 99: got %d mismatches, want 0" % misses299 for m in odds:300 assert mean(m) == Fraction(1, m), "mean s(%du): got %s, want 1/%d" % (m, mean(m), m)301 print("master law, all odd pairs 1..99: 0 mismatches; mean 1/m on all odd m to 99")302303 for m, n, want in [(3, 9, Fraction(20, 729)), (5, 15, Fraction(68, 1875)),304 (9, 15, Fraction(116, 18225)), (7, 21, Fraction(96, 2401)),305 (3, 5, Fraction(0)), (5, 7, Fraction(0))]:306 got = brute_carpet(m, n)307 assert got == want, "carpet cov(%d,%d): got %s, want %s" % (m, n, got, want)308 assert got == cov_carpet(m, n), "carpet closed form(%d,%d): got %s, want %s" % (309 m, n, cov_carpet(m, n), got)310 gn = brute_net(m, n)311 assert gn == cov_net(m, n), "net cov(%d,%d): got %s, want %s" % (m, n, gn, cov_net(m, n))312 gv = brute_void(m, n)313 assert gv == cov_void(m, n), "void cov(%d,%d): got %s, want %s" % (m, n, gv, cov_void(m, n))314 for m, n in [(3, 9), (5, 15), (9, 15), (7, 21)]:315 direct = float(cov_net(m, n)) / math.sqrt(float(cov_net(m, m)) * float(cov_net(n, n)))316 assert abs(direct - r_net(m, n)) < 1e-12, "net r(%d,%d): got %.12f, want %.12f" % (317 m, n, r_net(m, n), direct)318 print("two-dimensional cell counts on (3,9),(5,15),(9,15),(7,21),(3,5),(5,7): carpet, net and void exact")319320 zeros = 0321 for a in range(1, 20):322 for b in range(a, 20):323 m, n = 2 * a + 1, 2 * b + 1324 tree, carpet, net, void = assembled(m, n)325 for got, want, name in [(tree, cov_tree(m, n), "tree"), (carpet, cov_carpet(m, n), "carpet"),326 (net, cov_net(m, n), "net"), (void, cov_void(m, n), "void")]:327 assert got == want, "%s cov(%d,%d): got %s, want %s" % (name, m, n, got, want)328 coprime = math.gcd(m, n) == 1329 vanish = (tree == 0 and carpet == 0 and net == 0 and void == 0)330 positive = (tree > 0 and carpet > 0 and net > 0 and void > 0)331 assert vanish == coprime, "common zero set at (%d,%d): got %s, want %s" % (332 m, n, vanish, coprime)333 assert vanish or positive, "strict positivity at (%d,%d): got %s, want positive" % (334 m, n, (tree, carpet, net, void))335 if vanish:336 zeros += 1337 assert zeros == 139, "coprime pairs among odd 3..39: got %d, want 139" % zeros338 print("four families, odd pairs 3..39: closed forms exact, all vanish exactly on the 139 coprime pairs")339340 primes = 0341 composites = 0342 worst = None343 for n in range(3, 200, 2):344 score = 0.0345 for m in range(3, n, 2):346 score = max(score, r_carpet(m, n))347 if is_prime(n):348 primes += 1349 assert score == 0.0, "prime %d: got score %r, want 0.0" % (n, score)350 else:351 composites += 1352 assert score > 0.0, "composite %d: got score %r, want positive" % (n, score)353 if worst is None or score < worst[1]:354 worst = (n, score)355 assert primes == 45, "primes in odd 3..199: got %d, want 45" % primes356 assert composites == 54, "composites in odd 3..199: got %d, want 54" % composites357 assert worst[0] == 169, "narrowest composite: got %d, want 169" % worst[0]358 assert abs(worst[1] - 0.0517383422) < 1e-9, "carpet gap at 169: got %.10f, want 0.0517383422" % worst[1]359 tree_score = max(r_tree(m, 169) for m in range(3, 169, 2))360 void_score = max(r_void(m, 169) for m in range(3, 169, 2))361 assert abs(tree_score - 0.0766964989) < 1e-9, "tree score at 169: got %.10f, want 0.0766964989" % tree_score362 assert abs(void_score - 0.0059170562) < 1e-9, "void score at 169: got %.10f, want 0.0059170562" % void_score363 assert r_tree(3, 9) > r_net(3, 9) > r_carpet(3, 9) > r_void(3, 9), (364 "family ordering at (3,9): got %r, want decreasing" % (365 (r_tree(3, 9), r_net(3, 9), r_carpet(3, 9), r_void(3, 9)),))366 for got, want, name in [(r_tree(3, 9), 0.316227766, "tree"), (r_net(3, 9), 0.262950294, "net"),367 (r_carpet(3, 9), 0.219264505, "carpet"), (r_void(3, 9), 0.110431526, "void")]:368 assert abs(got - want) < 1e-9, "%s r(3,9): got %.9f, want %.9f" % (name, got, want)369 print("detector, odd 3..199: 45 primes at exactly 0, 54 composites positive, narrowest 0.0517383 at 169")370371 layers = list(range(1, 56, 2))372373 for n in range(1, 10, 2):374 for p in range(1, 9):375 for q in range(1, 9):376 if math.gcd(p, q) != 1:377 continue378 got = integral(n * p, n * q)379 want = Fraction(1, p * q) if p % 2 == 1 and q % 2 == 1 else Fraction(0)380 assert got == want, "origin ray slope %d/%d at layer %d: got %s, want %s" % (381 q, p, n, got, want)382 print("origin rays: slope q/p carries exactly 1/(pq) when p,q odd and exactly 0 otherwise")383384 for b in [1, 3, 5, 7, 9]:385 got = sum(tri(Fraction(n, b)) for n in layers) / len(layers)386 limit = ray(1, b) if b > 1 else Fraction(1)387 if b == 9:388 assert got > 0 > limit, "b=9 ray at 28 layers: got %s, want a sign disagreeing with %s" % (389 got, limit)390 assert got == Fraction(1, 63), "b=9 ray at 28 layers: got %s, want 1/63" % got391 if b == 7:392 assert got == limit, "b=7 ray at 28 layers: got %s, want %s" % (got, limit)393 if b == 5:394 assert abs(float(got / limit) - 1.4285714) < 1e-6, (395 "b=5 ray at 28 layers: got %.7f of the limit, want 1.4285714" % float(got / limit))396 got2 = sum(tri(Fraction(n, 2)) for n in layers) / len(layers)397 assert got2 == 0, "b=2 ray at 28 layers: got %s, want 0" % got2398 got6 = sum(tri(Fraction(n, 6)) for n in layers) / len(layers)399 assert got6 == Fraction(-1, 42), "b=6 ray at 28 layers: got %s, want -1/42" % got6400 wide = list(range(1, 2222, 2))401 got6w = sum(tri(Fraction(n, 6)) for n in wide) / len(wide)402 assert abs(float(got6w)) < 1e-3, "b=6 ray at 1111 layers: got %s, want under 1e-3" % got6w403 print("slope-one rays at 28 layers: b=7 exact at -1/49, b=9 at +1/63 against a limit -1/81, "404 "b=5 at 10/7 of its limit, b=6 at -1/42 against a limit 0")405406 for L, want in [(28, 0.21014), (101, 0.24516), (501, 0.26400)]:407 stack = list(range(1, 2 * L, 2))408 acc = Fraction(0)409 for m in stack:410 for n in stack:411 acc += cov_carpet(m, n)412 got = float(acc / L)413 print("stack variance at %d layers: L*Var = %.5f" % (L, got))414 assert abs(got - want) < 5e-6, "L*Var at %d layers: got %.5f, want %.5f" % (L, got, want)415416 for m in range(1, 17):417 for n in range(m, 17):418 g = math.gcd(m, n)419 got = franel(m, n)420 want = Fraction(g * g, 12 * m * n)421 assert got == want, "Franel integral (%d,%d): got %s, want %s" % (m, n, got, want)422 assert franel(3, 9) == Fraction(1, 36), "Franel (3,9): got %s, want 1/36" % franel(3, 9)423 print("classical sawtooth sibling, all pairs 1..16: integral is exactly gcd^2/(12mn)")424425 T = 0.0426 top = 1001427 for k in range(1, top + 1, 2):428 for l in range(1, top + 1, 2):429 T += 1.0 / (k * k * l * l * max(k, l))430 assert abs(T - 1.1122336970) < 1e-6, "T truncated at 1001: got %.10f, want 1.1122336970" % T431 print("gcd-sum constant T truncated at 1001: %.9f" % T)432433 mu = mobius_upto(200)434 seen = set()435 checked = 0436 for b in range(2, 201, 2):437 odd_scales = range(1, 2 * b, 2)438 for a in range(1, b):439 if math.gcd(a, b) != 1:440 continue441 stacked = Fraction(sum(tri_scaled(n * a, b) for n in odd_scales), b * b)442 assert stacked == 0, "triangle-wave stack at %d/%d over %d layers: got %s, want 0" % (443 a, b, b, stacked)444 assert stacked == ray(a, b), "ray law at %d/%d: got %s, want %s" % (a, b, ray(a, b), stacked)445 checked += 1446 seen.add(mu[b])447 assert checked == 4081, "reduced fractions with even denominator to 200: got %d, want 4081" % checked448 assert seen == {-1, 0, 1}, "Moebius values on even denominators: got %s, want {-1,0,1}" % sorted(seen)449 print("all 4081 even denominators to 200 carry stacked ray strength exactly 0 "450 "while their Moebius values range over -1, 0, 1")451452 const2, single2, pair2, _ = fourier(2)453 const3, single3, pair3, triple3 = fourier(3)454 assert const2 == 0.75 and single2 == [0.25, 0.25], "plane rule constants: got %r, want 0.75 and 0.25s" % (455 (const2, single2),)456 assert pair2 == [-0.25], "plane rule pair term: got %r, want [-0.25]" % pair2457 assert const3 == 0.5 and single3 == [0.25, 0.25, 0.25], "space rule constants: got %r, want 0.5 and 0.25s" % (458 (const3, single3),)459 assert pair3 == [0.0, 0.0, 0.0], "space rule pair terms: got %r, want zeros" % pair3460 assert triple3 == -0.25, "space rule triple term: got %r, want -0.25" % triple3461 print("fill rule harmonics: the plane carries a pair term -1/4, the space carries none")462463 S = [1, 3, 5, 7, 9]464 for a in range(1, 16):465 for b in range(1, 16):466 got = spectrum_direct(S, a, b)467 if a % 2 == 0 or b % 2 == 0:468 assert abs(got) < 1e-12, "spectrum: even index (%d,%d) got %r want 0" % (a, b, got)469 else:470 want = spectrum_formula(S, a, b)471 assert abs(got - want) < 1e-12 * max(1.0, abs(want)), \472 "spectrum (%d,%d): got %r want %r" % (a, b, got, want)473 got = spectrum_direct(S, 45, 45)474 want = spectrum_formula(S, 45, 45)475 assert abs(got - want) < 1e-14, "spectrum (45,45): got %r want %r" % (got, want)476 print("stack spectrum at L=5: cell-exact integrals match the divisor-sum formula, "477 "all pairs to 15 and (45,45); even indices exactly dark")478479 for L in (2, 3, 4, 5, 6, 8):480 S = list(range(1, 2 * L, 2))481 axis = Fraction(0)482 inter = Fraction(0)483 var = Fraction(0)484 for m in S:485 for n in S:486 d = math.gcd(m, n)487 axis += Fraction(2 * (m - 1) * (n - 1) * (d * d - 1), 16 * m * m * n * n)488 inter += Fraction((d * d - 1) ** 2, 16 * m * m * n * n)489 var += Fraction((d * d - 1) * (2 * (m - 1) * (n - 1) + d * d - 1), 16 * m * m * n * n)490 assert (axis + inter) * L * L == var * L * L and axis + inter == var, \491 "parseval split at L=%d: %r + %r != %r" % (L, axis, inter, var)492 print("variance splits exactly into axis and interior spectral blocks, L = 2..8")493494 A = 1501495 sig2 = [0] * (A + 1)496 for d in range(1, A + 1, 2):497 for mult in range(d, A + 1, 2 * d):498 sig2[mult] += d * d499 for w, use_sq, tol in ((3.0, False, 2e-3), (3.5, True, 2e-3)):500 brute = 0.0501 for a in range(1, A + 1, 2):502 for b in range(1, A + 1, 2):503 g = math.gcd(a, b)504 v = sig2[g]505 if use_sq:506 v *= sig2[g]507 brute += v / float(a * b) ** w508 if not use_sq:509 closed = lam(w) ** 2 * lam(2 * w - 2)510 else:511 closed = lam(w) ** 2 * lam(2 * w - 2) ** 2 * lam(2 * w - 4) / lam(4 * w - 4)512 assert abs(brute / closed - 1) < tol, \513 "zeta quotient w=%r sq=%r: %r vs %r" % (w, use_sq, brute, closed)514 print("zeta quotients: sigma_2 and sigma_2^2 gcd sums match the lambda products at w = 3, 3.5")515516 s = 3.0517 w = 2.0518 brute = 0.0519 for a in range(1, A + 1, 2):520 for b in range(1, A + 1, 2):521 g = math.gcd(a, b)522 sd = 0.0523 d = 1524 while d * d <= g:525 if g % d == 0:526 sd += d ** (2 - s)527 if d * d != g:528 sd += (g // d) ** (2 - s)529 d += 2530 brute += sd / float(a * b) ** w531 closed = lam(w) ** 2 * lam(2 * w + s - 2)532 assert abs(brute / closed - 1) < 2e-3, "weighted stack: %r vs %r" % (brute, closed)533 print("weighted stack at s=3: sigma_{-1} gcd sum matches lambda(w)^2 lambda(2w+1)")534535 print("all green")536537538if __name__ == "__main__":539 main()