stack_levels.py
13.1 kB · python · 347 lines
1from fractions import Fraction2from math import gcd34import numpy as np56# GENERAL LAW78def shadow_mask(base, keep, level):9 span = base ** level10 idx = np.arange(span)11 good = np.ones(span, dtype=bool)12 rest = idx.copy()13 for _ in range(level):14 good &= np.isin(rest % base, list(keep))15 rest //= base16 return good1718def cell_index(total, span, scale):19 step = total // (span * scale)20 return (np.arange(total) // step) % span2122def cross_mean(mask_f, mask_g, m, n):23 span_f = len(mask_f)24 span_g = len(mask_g)25 a = span_f * m26 b = span_g * n27 total = a * b // gcd(a, b)28 hit = mask_f[cell_index(total, span_f, m)] & mask_g[cell_index(total, span_g, n)]29 return Fraction(int(hit.sum()), total)3031def cross_cov(mask_f, mask_g, m, n):32 mu_f = Fraction(int(mask_f.sum()), len(mask_f))33 mu_g = Fraction(int(mask_g.sum()), len(mask_g))34 return cross_mean(mask_f, mask_g, m, n) - mu_f * mu_g3536# PARITY CHECK3738def parity_integral(m, n):39 total = m * n // gcd(m, n)40 i = np.arange(total)41 sm = 1 - 2 * (((m * i) // total) % 2)42 sn = 1 - 2 * (((n * i) // total) % 2)43 return Fraction(int((sm * sn).sum()), total)4445def parity_mean(m):46 i = np.arange(m)47 return Fraction(int((1 - 2 * ((( m * i) // m) % 2)).sum()), m)4849def parity_master(m, n):50 g = gcd(m, n)51 if (m // g) % 2 == 1 and (n // g) % 2 == 1:52 return Fraction(g * g, m * n)53 return Fraction(0)5455def check_parity():56 bad = [(m, n) for m in range(1, 41) for n in range(m, 41)57 if parity_integral(m, n) != parity_master(m, n)]58 assert bad == [], "master integral mismatches %s" % bad[:4]59 for m, n in [(3, 5), (3, 9), (5, 15)]:60 g = gcd(m, n)61 want = Fraction(g * g, m * n)62 got = parity_integral(m, n)63 assert got == want, "parity (%d,%d): got %s want %s" % (m, n, got, want)64 named = ", ".join("(%d,%d) %s" % (m, n, parity_integral(m, n))65 for m, n in [(3, 5), (3, 9), (5, 15)])66 print("parity master integral, all pairs to 40: 0 mismatches; %s" % named)67 odd = [n for n in range(1, 41) if n % 2 == 1]68 tree = []69 for m in odd:70 for n in odd:71 if n < m:72 continue73 g = gcd(m, n)74 cov = (parity_integral(m, n) - parity_mean(m) * parity_mean(n)) / 475 assert cov == Fraction(g * g - 1, 4 * m * n), "tree cov (%d,%d)" % (m, n)76 if g == 1:77 tree.append(cov)78 assert set(tree) == {Fraction(0)}, "tree coprime covariance not zero"79 strip = shadow_mask(2, {1}, 1)80 live = 081 coprime = 082 for m in range(1, 41):83 for n in range(m, 41):84 g = gcd(m, n)85 got = cross_cov(strip, strip, m, n)86 want = Fraction(g * g, 4 * m * n) if (m // g) % 2 and (n // g) % 2 else Fraction(0)87 assert got == want, "odd-strip cov (%d,%d): got %s want %s" % (m, n, got, want)88 if g == 1:89 coprime += 190 live += got != 091 print("tree parity layers: covariance (g^2-1)/(4mn), zero at all %d coprime odd pairs to 40"92 % len(tree))93 print("1-periodic odd strip: covariance g^2/(4mn), %s at (3,5), nonzero at %d of the %d coprime pairs to 40"94 % (cross_cov(strip, strip, 3, 5), live, coprime))9596# CARPET STACK9798def carpet_mask(level):99 return shadow_mask(3, {0, 2}, level)100101def chi3(k):102 r = k % 3103 return 0 if r == 0 else (1 if r == 1 else -1)104105def carpet_law_level1(m, n):106 g = gcd(m, n)107 a, b = m // g, n // g108 return Fraction(2 * chi3(a) * chi3(b), 9 * a * b)109110def carpet_kernel(level, lo=1, hi=41):111 q = 3 ** level112 mask = carpet_mask(level)113 table = {}114 hits = 0115 reduce_fail = 0116 zero_fail = 0117 live = 0118 for m in range(lo, hi):119 for n in range(m, hi):120 g = gcd(m, n)121 a, b = m // g, n // g122 cov = cross_cov(mask, mask, m, n)123 if cov != cross_cov(mask, mask, a, b):124 reduce_fail += 1125 if (cov == 0) != (a % q == 0 or b % q == 0):126 zero_fail += 1127 if g == 1 and cov != 0:128 live += 1129 key = (a % q, b % q)130 val = cov * a * b131 if key in table:132 hits += 1133 if table[key] != val:134 raise AssertionError("residue kernel breaks at %s" % (key,))135 table[key] = val136 return table, hits, reduce_fail, zero_fail, live137138def carpet_2d(level, m, n):139 mask = carpet_mask(level)140 span = len(mask)141 a = span * m142 b = span * n143 total = a * b // gcd(a, b)144 fm = mask[cell_index(total, span, m)]145 fn = mask[cell_index(total, span, n)]146 joint = int((np.outer(fm, fm) & np.outer(fn, fn)).sum())147 return Fraction(joint, total * total) - Fraction(int(fm.sum()) * int(fn.sum()), total * total) ** 2148149def check_carpet():150 mask = carpet_mask(1)151 for m in range(1, 41):152 for n in range(1, 41):153 got = cross_cov(mask, mask, m, n)154 want = carpet_law_level1(m, n)155 assert got == want, "carpet L=1 (%d,%d): got %s want %s" % (m, n, got, want)156 print("carpet L=1: covariance = (2/9) chi(m') chi(n')/(m'n') on all 1600 ordered pairs to 40; "157 "(1,2) %s, (2,5) %s, (1,3) %s"158 % (cross_cov(mask, mask, 1, 2), cross_cov(mask, mask, 2, 5), cross_cov(mask, mask, 1, 3)))159 tables = {}160 for level in (1, 2, 3):161 table, hits, reduce_fail, zero_fail, live = carpet_kernel(level)162 assert reduce_fail == 0, "reduction fails at L=%d" % level163 assert zero_fail == 0, "zero law fails at L=%d" % level164 tables[level] = table165 print("carpet L=%d: 820 pairs to 40, reduction Cov(m,n)=Cov(m/g,n/g) exact, "166 "kernel mod %d agrees on %d repeated residue keys, zero set is 3^%d dividing m' or n', "167 "%d of the 490 coprime pairs carry nonzero covariance"168 % (level, 3 ** level, hits, level, live))169 units = [r for r in range(9) if r % 3]170 rows = " ; ".join("%d: %s" % (u, " ".join(str(tables[2][(min(u, v), max(u, v))]) for v in units))171 for u in units)172 print("carpet L=2 kernel G_2(a,b) = Cov * a b on units mod 9, rows a = %s" % rows)173 sep = tables[2][(1, 1)] * tables[2][(4, 4)] - tables[2][(1, 4)] ** 2174 print("carpet L=2 kernel is not separable: G(1,1)G(4,4) - G(1,4)^2 = %s, so no theta(a)theta(b) form"175 % sep)176 row3 = " ".join(str(tables[3][(1, v)]) for v in range(1, 27) if v % 3)177 print("carpet L=3 kernel row G_3(1,b) over units b mod 27: %s" % row3)178 for level in (1, 2):179 for m, n in [(1, 2), (2, 5), (1, 3), (4, 7)]:180 got = carpet_2d(level, m, n)181 one = cross_mean(carpet_mask(level), carpet_mask(level), m, n)182 mu = Fraction(2, 3) ** (2 * level)183 assert got == (one - mu) * (one + mu), "2d carpet (%d,%d) L=%d" % (m, n, level)184 print("carpet 2D field: cell-counted covariance equals Cov_1D times (mean + (2/3)^(2L)) at "185 "(1,2),(2,5),(1,3),(4,7) and L=1,2, so the 2D and 1D zero sets coincide")186187# LEVELS188189def level_product_mask(level):190 top = carpet_mask(1)191 span = 3 ** level192 out = np.ones(span, dtype=bool)193 for i in range(level):194 out &= top[(np.arange(span) // 3 ** (level - 1 - i)) % 3]195 return out196197def check_levels():198 for level in (1, 2, 3, 4):199 assert (level_product_mask(level) == carpet_mask(level)).all(), "level split at %d" % level200 nested = all(bool(carpet_mask(3)[i]) <= bool(carpet_mask(2)[i % 9]) for i in range(27))201 assert nested, "f_3 not contained in f_2(3x)"202 print("levels: f_L(x) = product of f_1(3^i x) over i < L exact as masks at L = 1,2,3,4; "203 "f_L(nx) <= f_(L-1)(3nx) pointwise, the gap carrying measure (2/3)^(L-1)/3")204 masks = {level: carpet_mask(level) for level in (1, 2, 3)}205 breaches = 0206 unpredicted = 0207 total = 0208 for level in (1, 2, 3):209 for other in (1, 2, 3):210 for m in range(1, 41):211 for n in range(1, 41):212 g = gcd(m, n)213 a, b = m // g, n // g214 cov = cross_cov(masks[level], masks[other], m, n)215 dead = b % 3 ** level == 0 or a % 3 ** other == 0216 total += 1217 if dead and cov != 0:218 breaches += 1219 if not dead and cov == 0:220 unpredicted += 1221 assert breaches == 0 and unpredicted == 0, "cross-level zero law breached"222 print("levels: Cov(f_L(mx), f_L'(nx)) = 0 exactly when 3^L divides n/g or 3^L' divides m/g, "223 "on all %d ordered triples (L,L' in 1..3, m,n to 40); the two levels enter asymmetrically"224 % total)225 print("levels: cross-level examples Cov(f_1(x), f_2(3x)) = %s and Cov(f_2(x), f_1(3x)) = %s, "226 "so level L' at scale 3n is not redundant against level L at scale n"227 % (cross_cov(masks[1], masks[2], 1, 3), cross_cov(masks[2], masks[1], 1, 3)))228229# DESIGNS230231def design_mask(base, removed):232 return shadow_mask(base, set(range(base)) - {removed}, 1)233234def design_law_33(removed_f, removed_g, m, n):235 g = gcd(m, n)236 a, b = m // g, n // g237 if a % 3 == 0 or b % 3 == 0:238 return Fraction(0)239 root = complex(-0.5, 3 ** 0.5 / 2)240 w = root ** ((a * removed_g - b * removed_f) % 3)241 w *= (1 - root ** ((-b) % 3)) * (1 - root ** (a % 3))242 return Fraction(round(2 * w.real), 27 * a * b)243244def design_law_23(removed_g, m, n):245 g = gcd(m, n)246 a, b = m // g, n // g247 if b % 2 == 0 or a % 3 == 0:248 return Fraction(0)249 root = complex(-0.5, 3 ** 0.5 / 2)250 u = root ** ((a * removed_g) % 3) * (1 - root ** (a % 3))251 return Fraction(round(2 * u.real), 18 * a * b)252253def check_designs():254 strip = shadow_mask(2, {1}, 1)255 base3 = {removed: design_mask(3, removed) for removed in (0, 1, 2)}256 seen = set()257 for rf in (0, 1, 2):258 for rg in (0, 1, 2):259 for m in range(1, 41):260 for n in range(1, 41):261 got = cross_cov(base3[rf], base3[rg], m, n)262 want = design_law_33(rf, rg, m, n)263 assert got == want, "base-3 pair (%d,%d) at (%d,%d)" % (rf, rg, m, n)264 g = gcd(m, n)265 if got != 0:266 seen.add(got * 27 * (m // g) * (n // g))267 print("designs, base 3 level 1: Cov = c/(27 m'n') with c in %s, zero exactly when 3 divides m'n', "268 "on all 9 ordered design pairs and 1600 scale pairs to 40"269 % sorted(int(v) for v in seen))270 zero = 0271 for rg in (0, 1, 2):272 for m in range(1, 41):273 for n in range(1, 41):274 got = cross_cov(strip, base3[rg], m, n)275 assert got == design_law_23(rg, m, n), "strip vs base-3 %d at (%d,%d)" % (rg, m, n)276 if rg == 1 and got == 0:277 zero += 1278 print("designs, base 2 against base 3: Cov(p(mx), h_r(nx)) = c/(18 m'n') with c in -3,0,3, "279 "zero when n/g is even or 3 divides m/g")280 print("designs: the odd strip is orthogonal to the Sierpinski shadow at every scale pair, "281 "%d of %d pairs exactly zero, because the centred shadow is even and the centred strip odd "282 "under x -> -x" % (zero, 1600))283284# GRAM285286def jordan2(k):287 out = k * k288 d = 2289 r = k290 while d * d <= r:291 if r % d == 0:292 out = out // (d * d) * (d * d - 1)293 while r % d == 0:294 r //= d295 d += 1296 if r > 1:297 out = out // (r * r) * (r * r - 1)298 return out299300def exact_det(rows):301 size = len(rows)302 work = [list(row) for row in rows]303 det = Fraction(1)304 for col in range(size):305 pivot = next((r for r in range(col, size) if work[r][col] != 0), None)306 if pivot is None:307 return Fraction(0)308 if pivot != col:309 work[col], work[pivot] = work[pivot], work[col]310 det = -det311 det *= work[col][col]312 inv = Fraction(1) / work[col][col]313 for r in range(col + 1, size):314 factor = work[r][col] * inv315 if factor:316 for c in range(col, size):317 work[r][c] -= factor * work[col][c]318 return det319320def check_gram():321 for cap in range(1, 13):322 odds = list(range(1, 2 * cap + 2, 2))323 rows = [[Fraction(gcd(a, b) ** 2, a * b) for b in odds] for a in odds]324 got = exact_det(rows)325 want = Fraction(1)326 for k in odds:327 want *= Fraction(jordan2(k), k * k)328 assert got == want, "gram det at K=%d: got %s want %s" % (cap, got, want)329 assert got > 0, "gram det not positive at K=%d" % cap330 odds = list(range(1, 26, 2))331 final = Fraction(1)332 for k in odds:333 final *= Fraction(jordan2(k), k * k)334 print("gram: det[gcd(m,n)^2/(mn)] over odd m,n <= 2K+1 equals prod J_2(k)/k^2 over the same odds, "335 "exact at K = 1..12, so Smith's factor-closed determinant applies to the odd index set")336 print("gram: the value is prod over odd k <= 2K+1 of prod over p | k of (1 - p^-2); at K = 12 it is "337 "%s, positive, so the odd parity layers are linearly independent in L^2" % final)338339def main():340 check_parity()341 check_carpet()342 check_levels()343 check_designs()344 check_gram()345346if __name__ == "__main__":347 main()