walsh.py
16.5 kB · python · 401 lines
1import math2import sys3import time4from fractions import Fraction5from functools import lru_cache67import numpy as np89# CORNERS1011def corners(code, dim):12 return [i for i in range(2 ** dim) if code >> i & 1]1314def weight(i):15 return bin(i).count("1")1617def profile(code, dim):18 a = [0] * (dim + 1)19 for i in corners(code, dim):20 a[weight(i)] += 121 return a2223def code_of(cs):24 return sum(1 << c for c in cs)2526def mirror(code, dim):27 full = 2 ** dim - 128 return code_of(full ^ c for c in corners(code, dim))2930def complement(code, dim):31 return (1 << 2 ** dim) - 1 - code3233def flip(code, dim, c):34 return code_of(i ^ c for i in corners(code, dim))3536def bits(dim):37 codes = np.arange(2 ** (2 ** dim), dtype=np.int64)38 return (codes[:, None] >> np.arange(2 ** dim)[None, :]) & 13940# WALSH4142def hadamard(dim):43 s = np.arange(2 ** dim)44 pc = np.array([weight(int(v)) for v in (s[:, None] & s[None, :]).ravel()]).reshape(2 ** dim, 2 ** dim)45 return (-1) ** pc4647def spectra(dim):48 return bits(dim) @ hadamard(dim)4950def levels(h, dim):51 lv = np.array([weight(s) for s in range(2 ** dim)])52 return np.stack([h[..., lv == j].sum(axis=-1) for j in range(dim + 1)], axis=-1)5354def square_levels(h, dim):55 return levels(h * h, dim)5657def evaluate(lev, dim, n):58 return sum(int(lev[j]) * n ** (dim - j) for j in range(dim + 1))5960# RENDER6162@lru_cache(maxsize=None)63def parity_index(n, dim):64 grid = np.indices((n,) * dim).reshape(dim, -1) % 265 return sum(grid[k] << (dim - 1 - k) for k in range(dim))6667def histogram(n, dim):68 return np.bincount(parity_index(n, dim), minlength=2 ** dim).astype(np.int64)6970def literal(code, n, dim):71 table = np.array([code >> i & 1 for i in range(2 ** dim)], dtype=bool)72 return int(np.count_nonzero(table[parity_index(n, dim)]))7374def section_expansion():75 bad = [0, 0]76 checks = [0, 0]77 for dim in (1, 2, 3):78 lev = levels(spectra(dim), dim)79 for code in range(1, 2 ** (2 ** dim)):80 w = bin(code).count("1")81 for n in range(1, 10):82 f = literal(code, n, dim)83 want = evaluate(lev[code], dim, n) if n % 2 else w * n ** dim84 checks[n % 2] += 185 bad[n % 2] += 2 ** dim * f != want86 rng = np.random.default_rng(7919)87 sample = rng.choice(np.arange(1, 2 ** 16), size=2000, replace=False)88 lev4 = levels(spectra(4), 4)89 for code in sample:90 code = int(code)91 w = bin(code).count("1")92 for n in range(1, 10):93 f = literal(code, n, 4)94 want = evaluate(lev4[code], 4, n) if n % 2 else w * n ** 495 checks[n % 2] += 196 bad[n % 2] += 16 * f != want97 b4 = bits(4)98 for n in range(1, 10):99 f = b4 @ histogram(n, 4)100 if n % 2:101 want = sum(lev4[:, j] * n ** (4 - j) for j in range(5))102 else:103 want = b4.sum(axis=1) * n ** 4104 checks[n % 2] += len(f)105 bad[n % 2] += int(np.count_nonzero(16 * f != want))106 print(f"expansion: 2^dim fill(N) = sum_j 2^dim W_j N^(dim-j) at odd N 1..9, {checks[1]} checks, {bad[1]} mismatches; w N^dim at even N 2..8, {checks[0]} checks, {bad[0]} mismatches")107 print(" literal grids: all 273 nonempty codes at dim 1..3 and 2000 seeded codes at dim 4, sides 1..9; literal corner histogram x all 65536 codes at dim 4")108 for dim, code in ((1, 1), (2, 7), (2, 11), (3, 23), (3, 232)):109 lev = levels(spectra(dim), dim)[code]110 ws = ", ".join(str(Fraction(int(x), 2 ** dim)) for x in lev)111 print(f" dim {dim} code {code}: W = ({ws}), 2^dim fill = {poly_text(lev, dim)}, fill(3) = {literal(code, 3, dim)}")112113def poly_text(lev, dim):114 terms = []115 for j in range(dim + 1):116 c = int(lev[j])117 if c == 0:118 continue119 p = dim - j120 mono = "" if p == 0 else ("N" if p == 1 else f"N^{p}")121 coef = str(abs(c)) if (abs(c) != 1 or p == 0) else ""122 body = coef + ("*" if coef and mono else "") + mono123 terms.append(("- " if c < 0 else "+ ") + body)124 s = " ".join(terms)125 return s[2:] if s.startswith("+ ") else "-" + s[2:]126127# PROFILE128129def krawtchouk(j, w, dim):130 return sum((-1) ** i * math.comb(w, i) * math.comb(dim - w, j - i) for i in range(j + 1))131132def section_profile():133 bad = total = 0134 for dim in (1, 2, 3, 4):135 lev = levels(spectra(dim), dim)136 b = bits(dim)137 wt = np.array([weight(i) for i in range(2 ** dim)])138 a = np.stack([b[:, wt == w].sum(axis=1) for w in range(dim + 1)], axis=1)139 K = np.array([[krawtchouk(j, w, dim) for w in range(dim + 1)] for j in range(dim + 1)])140 bad += int(np.count_nonzero(a @ K.T != lev))141 total += lev.size142 assert np.array_equal(K @ K, 2 ** dim * np.eye(dim + 1, dtype=int))143 print(f"profile: 2^dim W_j = sum_w a_w K_j(w) on every code at dim 1..4, {total} level sums, {bad} mismatches; K^2 = 2^dim I at dim 1..4")144145# MIRROR146147def section_mirror():148 for dim in (1, 2, 3, 4):149 lev = levels(spectra(dim), dim)150 sign = np.array([(-1) ** (dim - j) for j in range(dim + 1)])151 mbad = 0152 swap, crit, selfdual, other = set(), set(), set(), set()153 for code in range(1, 2 ** (2 ** dim)):154 neg = sign * lev[code]155 mbad += not np.array_equal(neg, (-1) ** dim * lev[mirror(code, dim)])156 void = -lev[code].copy()157 void[0] += 2 ** dim158 if np.array_equal(neg, (-1) ** dim * void):159 swap.add(code)160 if np.array_equal(neg, -((-1) ** dim) * void):161 other.add(code)162 a = profile(code, dim)163 if all(a[j] + a[dim - j] == math.comb(dim, j) for j in range(dim + 1)):164 crit.add(code)165 if complement(code, dim) == mirror(code, dim):166 selfdual.add(code)167 formula = 1168 for j in range(dim + 1):169 m = math.comb(dim, j)170 if 2 * j < dim:171 formula *= math.comb(2 * m, m)172 elif 2 * j == dim:173 formula *= math.comb(m, m // 2) if m % 2 == 0 else 0174 nonself = sorted(swap - selfdual)175 print(f"mirror dim {dim}: fill_F(-N) = (-1)^dim fill_F'(N) fails on {mbad} codes; swap {len(swap)}, criterion {len(crit)}, equal {swap == crit}, formula {formula}, self-dual {len(selfdual)} all inside {selfdual <= swap}, other sign {len(other)}")176 if dim <= 3:177 print(f" swapping, not self-dual: {len(nonself)}, least {nonself[:6]}")178 if dim == 2:179 print(f" dim 2 swapping codes {sorted(swap)}; code 7 swaps: {7 in swap}; code 3 swaps: {3 in swap}")180 if dim == 3:181 print(f" dim 3 code 1 swaps: {1 in swap}; code 23 swaps: {23 in swap}; code 27 swaps: {27 in swap}, self-dual: {27 in selfdual}, corners {[format(c, '03b') for c in corners(27, 3)]}, profile {profile(27, 3)}")182183# ROOTS184185def k_coefficients(code, dim):186 a = profile(code, dim)187 c = [0] * (dim + 1)188 for j, aj in enumerate(a):189 poly = np.poly1d([1.0])190 for _ in range(dim - j):191 poly = poly * np.poly1d([1.0, 0.0])192 for _ in range(j):193 poly = poly * np.poly1d([1.0, -1.0])194 c = np.polyadd(c, aj * poly.coeffs)195 return np.poly1d(c)196197def section_roots():198 worst = 0.0199 seen = set()200 count = 0201 for dim in (1, 2, 3):202 lev = levels(spectra(dim), dim)203 for code in range(1, 2 ** (2 ** dim)):204 key = (dim, tuple(profile(code, dim)))205 if key in seen:206 continue207 seen.add(key)208 count += 1209 P = k_coefficients(code, dim)210 R = np.poly1d([float(lev[code][j]) for j in range(dim, -1, -1)])211 scale = float(np.abs(lev[code]).sum())212 for r in P.roots:213 if abs(r - 0.5) < 1e-9:214 continue215 mu = 1 / (2 * r - 1)216 worst = max(worst, abs(R(mu)) / (scale * max(1.0, abs(mu)) ** dim))217 grid = np.linspace(-0.999, 0.999, 2001)218 assert np.all(R(grid) > 0)219 lev = levels(spectra(2), 2)[7]220 R = np.poly1d([float(lev[j]) for j in range(2, -1, -1)])221 print(f"roots: {count} profiles at dim 1..3, every root r != 1/2 of P_F(n) sends mu = 1/(2r-1) to a zero of R(mu) = sum W_j mu^j, worst scaled residual {worst:.1e}; R > 0 on (-1, 1) at every profile")222 print(f" dim 2 code 7: R roots {sorted(np.round(R.roots, 12))}, P_F roots {sorted(np.round(k_coefficients(7, 2).roots, 12))}")223224# DRIFT225226def downset(code, dim):227 cs = set(corners(code, dim))228 return all((c & ~(1 << b)) in cs for c in cs for b in range(dim) if c >> b & 1)229230def bichromatic(code, dim):231 cs = set(corners(code, dim))232 return sum(1 for c in range(2 ** dim) for b in range(dim) if not c >> b & 1 and ((c in cs) != ((c | 1 << b) in cs)))233234def section_drift():235 bad = total = 0236 dbad = off = above = 0237 dcount = []238 for dim in (1, 2, 3, 4):239 lev = levels(spectra(dim), dim)240 dcount.append(0)241 for code in range(1, 2 ** (2 ** dim)):242 a = profile(code, dim)243 w = sum(a)244 mean = Fraction(sum(j * aj for j, aj in enumerate(a)), w)245 drift = Fraction(dim, 2) - mean246 total += 1247 bad += drift != Fraction(int(lev[code][1]), 2 * int(lev[code][0]))248 u = bichromatic(code, dim)249 above += int(lev[code][1]) > u250 if downset(code, dim):251 dcount[-1] += 1252 dbad += int(lev[code][1]) != u253 else:254 off += int(lev[code][1]) != u255 print(f"drift: dim/2 - mean = W_1/(2 W_0) on all {total} nonempty codes at dim 1..4, {bad} mismatches")256 print(f" I = 2^-dim sum_x s(f, x) = U/2^(dim-1): 2^dim W_1 > U on {above} codes; 2^dim W_1 = U on the nonempty down-sets {dcount} at dim 1..4, {dbad} mismatches; off down-sets it fails on {off} codes")257258# STABILITY259260def section_stability():261 bad = checks = 0262 for dim in (1, 2, 3):263 h = spectra(dim)264 sq = square_levels(h, dim)265 for n in range(1, 10, 2):266 hist = histogram(n, dim)267 for code in range(1, 2 ** (2 ** dim)):268 fills = [int(sum(hist[i] for i in corners(flip(code, dim, c), dim))) for c in range(2 ** dim)]269 inside = sum(fills[c] for c in corners(code, dim))270 checks += 2271 bad += 2 ** dim * inside != evaluate(sq[code], dim, n)272 bad += sum(fills) != bin(code).count("1") * n ** dim273 polys = sorted({poly_text(levels(spectra(3), 3)[flip(23, 3, c)], 3) for c in range(8)})274 print(f"stability: 2^dim sum_(c in F) fill_(F+c)(N) = sum_j 2^(2 dim) W^j N^(dim-j) and sum over all c = w N^dim, every code at dim 1..3, odd sides 1..9; {checks} checks, {bad} mismatches")275 print(f" dim 3 code 23 flip orbit, 8 fills of {len(polys)} polynomials in N (times 8): {polys}")276277# ENTROPY278279def section_entropy():280 for dim, code in ((1, 1), (1, 2), (2, 7), (3, 23)):281 lev = levels(spectra(dim), dim)[code]282 w = bin(code).count("1")283 inf = math.log2(2 ** dim / w)284 row = []285 for n in (2, 3, 9, 27, 81, 243, 729):286 ratio = Fraction(evaluate(lev, dim, n) if n % 2 else w * n ** dim, 2 ** dim * n ** dim)287 row.append(f"{n}: {-math.log2(ratio):.9f}")288 slope = (int(lev[1]) / int(lev[0])) / math.log(2)289 n = 729290 ratio = Fraction(evaluate(lev, dim, n), 2 ** dim * n ** dim)291 print(f"entropy dim {dim} code {code}: bits per level {', '.join(row)}, infinity {inf:.9f}; N (limit - cost) at 729 {n * (inf + math.log2(ratio)):.6f} against (W_1/W_0)/ln 2 = {slope:.6f}")292293# THRESHOLD294295def partial(n, t):296 return sum(math.comb(n, j) for j in range(t + 1))297298def level_code(dim, allowed):299 return code_of(i for i in range(2 ** dim) if weight(i) in allowed)300301def section_threshold(top=2048):302 bad = 0303 for dim in range(1, 9):304 H = hadamard(dim)305 for t in range(dim + 1):306 code = level_code(dim, set(range(t + 1)))307 vec = np.array([code >> i & 1 for i in range(2 ** dim)], dtype=np.int64)308 l = levels(vec @ H, dim)309 bad += int(l[0]) != partial(dim, t)310 bad += int(l[1]) != (t + 1) * math.comb(dim, t + 1)311 print(f"threshold: at most t odd has 2^dim W_0 = S(dim, t), 2^dim W_1 = (t+1) C(dim, t+1) at dim 1..8, every t, {bad} mismatches")312 maj = []313 for dim in range(1, 9):314 r = Fraction(partial(dim, dim // 2), 2 ** dim)315 want = Fraction(1, 2) if dim % 2 else Fraction(1, 2) + Fraction(math.comb(dim, dim // 2), 2 ** (dim + 1))316 maj.append(f"{r}{'' if r == want else ' MISMATCH'}")317 print(f" majority, at most dim/2 odd, ratio at dim 1..8: {', '.join(maj)}")318 lev4 = levels(spectra(4), 4)319 one = level_code(4, {0, 1})320 two = level_code(4, {2, 3, 4})321 print(f" dim 4: at most one odd code {one} ratio {Fraction(int(lev4[one][0]), 16)}, 2^4 fill = {poly_text(lev4[one], 4)}; at least two odd code {two} ratio {Fraction(int(lev4[two][0]), 16)}, is the mirror of at most two odd: {two == mirror(level_code(4, {0, 1, 2}), 4)}")322 for dim, t in ((2, 1), (3, 1), (4, 2)):323 code = level_code(dim, set(range(t + 1)))324 print(f" dim {dim} at most {t} odd, code {code}: 2^dim fill = {poly_text(levels(spectra(dim), dim)[code], dim)}")325 groups = {}326 row = [1]327 for n in range(1, top + 1):328 row = [1] + [row[i] + row[i + 1] for i in range(len(row) - 1)] + [1]329 s = 0330 for t in range(n):331 s += row[t]332 v = (s & -s).bit_length() - 1333 groups.setdefault((s >> v, n - v), []).append((n, t))334 assert partial(274, 52) == 8 * partial(271, 51)335 half = groups.pop((1, 1))336 odd_majority = half == [(n, (n - 1) // 2) for n in range(1, top + 1, 2)]337 shared = {k: m for k, m in groups.items() if len(m) > 1}338 sizes = sorted({len(m) for m in shared.values()})339 print(f" coincidences of W_0 = S(dim, t)/2^dim among 0 <= t < dim <= {top}: 1/2 is held by exactly the odd majorities: {odd_majority}; {len(shared)} further shared values, group sizes {sizes}:")340 for k, m in sorted(shared.items(), key=lambda kv: kv[1][0]):341 odd, e = k342 if odd == 1:343 label = f"2^-{e}"344 elif odd == (1 << e) - 1:345 label = f"1 - 2^-{e}"346 else:347 label = f"odd/2^{e}, odd of {odd.bit_length()} bits"348 print(f" {label}: {m}")349350# FLAT351352def section_flat():353 counts = []354 for dim in (1, 2, 3, 4):355 lev = levels(spectra(dim), dim)356 free = set(np.nonzero(np.all(lev[:, 1:dim] == 0, axis=1))[0].tolist()) if dim > 1 else set(range(2 ** (2 ** dim)))357 crit = set()358 formula_bad = 0359 for code in range(2 ** (2 ** dim)):360 a = profile(code, dim)361 ev = {Fraction(a[w], math.comb(dim, w)) for w in range(0, dim + 1, 2)}362 od = {Fraction(a[w], math.comb(dim, w)) for w in range(1, dim + 1, 2)}363 if len(ev) <= 1 and len(od) <= 1:364 crit.add(code)365 alpha = ev.pop() if ev else Fraction(0)366 beta = od.pop() if od else Fraction(0)367 l = lev[code]368 formula_bad += Fraction(int(l[0]), 2 ** dim) != (alpha + beta) / 2369 formula_bad += Fraction(int(l[dim]), 2 ** dim) != (alpha - beta) / 2 if dim > 0 else 0370 counts.append(len(free))371 print(f"flat dim {dim}: W_j = 0 for 0 < j < dim on {len(free)} of {2 ** (2 ** dim)} codes, level-parity criterion {len(crit)}, equal {free == crit}, W_0 = (alpha+beta)/2 and W_dim = (alpha-beta)/2 fails {formula_bad}")372 if dim == 2:373 print(f" dim 2 lower-order-free codes {sorted(free)}; code 11 fill: 2^2 fill = {poly_text(lev[11], 2)}")374 if dim == 4:375 h = spectra(4)376 g = -2 * h377 g[:, 0] += 16378 bent = set(np.nonzero(np.all(np.abs(g) == 4, axis=1))[0].tolist())379 x = code_of(c for c in range(16) if ((c >> 3 & 1) & (c >> 2 & 1)) ^ ((c >> 1 & 1) & (c & 1)))380 print(f" dim 4: {len(bent)} bent codes, {len(bent & free)} of them lower-order-free; x1 x2 + x3 x4 is code {x}, bent {x in bent}, profile {profile(x, 4)}, 2^4 fill = {poly_text(lev[x], 4)}")381382# MAIN383384SECTIONS = {385 "expansion": section_expansion,386 "profile": section_profile,387 "mirror": section_mirror,388 "roots": section_roots,389 "drift": section_drift,390 "stability": section_stability,391 "entropy": section_entropy,392 "threshold": section_threshold,393 "flat": section_flat,394}395396if __name__ == "__main__":397 verbs = sys.argv[1:] or list(SECTIONS)398 for verb in verbs:399 start = time.perf_counter()400 SECTIONS[verb]()401 print(f" [{verb} {time.perf_counter() - start:.1f}s]")