weighted.py
27.5 kB · python · 730 lines
1from decimal import Decimal, ROUND_CEILING, ROUND_FLOOR, getcontext2from fractions import Fraction3from math import comb, e, gcd, isqrt, log, pi4import bisect56getcontext().prec = 8078# THE OBJECT910Q = 311D = 212CELLS = ((0, 0), (2, 0), (0, 2))13WEIGHT = (Fraction(3, 8), Fraction(3, 8), Fraction(1, 4))14LEVEL = 1015BINS = 2416FOLD_LO, FOLD_HI = 4, 817MASS_A = (20, 40, 60, 80, 100)18MASS_MAX = 130.019SPECTRUM_LEVELS = (6, 8, 10)20GRIPENBERG_LEN = 821NORM_LEN = 142223PROBS = (24 ("1/3,1/3,1/3", (Fraction(1, 3), Fraction(1, 3), Fraction(1, 3))),25 ("2/5,2/5,1/5", (Fraction(2, 5), Fraction(2, 5), Fraction(1, 5))),26 ("3/8,3/8,1/4", (Fraction(3, 8), Fraction(3, 8), Fraction(1, 4))),27 ("1/2,1/4,1/4", (Fraction(1, 2), Fraction(1, 4), Fraction(1, 4))),28)2930RENEWAL = (31 ("1/2,1/2,1/2", (Fraction(1, 2), Fraction(1, 2), Fraction(1, 2))),32 ("1/2,1/2,1/4", (Fraction(1, 2), Fraction(1, 2), Fraction(1, 4))),33 ("1/2,1/2,1/3", (Fraction(1, 2), Fraction(1, 2), Fraction(1, 3))),34 ("1/3,1/3,1/3", (Fraction(1, 3), Fraction(1, 3), Fraction(1, 3))),35 ("1/2,1/4,1/4", (Fraction(1, 2), Fraction(1, 4), Fraction(1, 4))),36 ("2/5,2/5,1/5", (Fraction(2, 5), Fraction(2, 5), Fraction(1, 5))),37 ("3/8,3/8,1/4", (Fraction(3, 8), Fraction(3, 8), Fraction(1, 4))),38)3940CONTROLS = ("1/3,1/3,1/3", "1/2,1/4,1/4", "3/8,3/8,1/4")4142NAPKIN_LENGTH = {43 "1/3,1/3,1/3": ("0.652441", "0.022937", "0.000000", "0.0208 0.0094 0.0047 0.0021"),44 "2/5,2/5,1/5": ("0.536956", "0.019130", "0.066973", "0.0224 0.0091 0.0043 0.0022"),45 "3/8,3/8,1/4": ("0.580858", "0.020478", "0.044604", "0.0205 0.0092 0.0041 0.0022"),46 "1/2,1/4,1/4": ("0.427372", "0.014471", "0.153880", "0.0299 0.0153 0.0077 0.0039"),47}4849NAPKIN_MASS = {50 "1/2,1/2,1/2": ("1.5849625007", "1.09812 1.09852 1.09812 1.09812 1.09803"),51 "1/2,1/2,1/4": ("1.2715533032", "0.88098 0.88130 0.88098 0.88098 0.88091"),52 "1/2,1/2,1/3": ("1.3646005647", "0.20259 0.14175 0.11707 0.10183 0.09030"),53 "1/3,1/3,1/3": ("1.0000000000", "1.09825 1.09825 1.09825 1.09825 1.09825"),54 "1/2,1/4,1/4": ("1.0000000000", "0.69284 0.69309 0.69284 0.69284 0.69278"),55 "2/5,2/5,1/5": ("1.0000000000", "0.23297 0.19792 0.18623 0.17623 0.16776"),56 "3/8,3/8,1/4": ("1.0000000000", "0.23133 0.16019 0.13671 0.11972 0.10734"),57}5859NAPKIN_RUNG0 = {60 "1/3,1/3,1/3": ("[3, 3, 3]", "3", "1.000000000", "1.000000000", "-1.000000000"),61 "2/5,2/5,1/5": ("[18/5, 18/5, 9/5]", "18/5", "0.834043767", "1.464973521", "-1.165956233"),62 "3/8,3/8,1/4": ("[27/8, 27/8, 9/4]", "27/8", "0.892789261", "1.261859507", "-1.107210739"),63 "1/2,1/4,1/4": ("[9/2, 9/4, 9/4]", "9/2", "0.630929754", "1.261859507", "-1.369070246"),64}6566NAPKIN_RUNG1 = {67 "hat (1/2, 1, 1/2)": ("0.5000000", "0.5000000", "1.0000000", "1.0000000"),68 "D4": ("0.6830127", "0.7105812", "0.4929285", "0.5500157"),69}7071FAILURES = []7273def safe_log2(x, digits, mode):74 v = -(Decimal(x.numerator) / Decimal(x.denominator)).ln() / Decimal(2).ln()75 return str(v.quantize(Decimal(1).scaleb(-digits), rounding=mode))7677def check(tag, got, want):78 if got != want:79 FAILURES.append("%s reads %s wants %s" % (tag, got, want))80 return got8182# THE ARITHMETIC CLASS, EXACT8384def factor(n):85 out, p = {}, 286 while p * p <= n:87 while n % p == 0:88 out[p] = out.get(p, 0) + 189 n //= p90 p += 1 if p == 2 else 291 if n > 1:92 out[n] = out.get(n, 0) + 193 return out9495def exponent_rows(w):96 primes = set()97 for x in w:98 primes |= set(factor(x.numerator)) | set(factor(x.denominator))99 primes = sorted(primes)100 rows = []101 for x in w:102 up, dn = factor(x.numerator), factor(x.denominator)103 rows.append(tuple(up.get(p, 0) - dn.get(p, 0) for p in primes))104 return primes, rows105106def rank(rows):107 m = [[Fraction(v) for v in r] for r in rows]108 top = 0109 for col in range(len(m[0])):110 p = next((k for k in range(top, len(m)) if m[k][col] != 0), None)111 if p is None:112 continue113 m[top], m[p] = m[p], m[top]114 d = m[top][col]115 m[top] = [v / d for v in m[top]]116 for k in range(len(m)):117 if k != top and m[k][col] != 0:118 f = m[k][col]119 m[k] = [a - f * b for a, b in zip(m[k], m[top])]120 top += 1121 return top122123def arithmetic_class(w):124 primes, rows = exponent_rows(w)125 r = rank(rows)126 if r != 1:127 return "nonlattice", None128 seed = next(v for v in rows if any(v))129 g = 0130 for v in seed:131 g = gcd(g, abs(v))132 unit = tuple(v // g for v in seed)133 pos = next(i for i, b in enumerate(unit) if b)134 mult = 0135 for v in rows:136 mult = gcd(mult, abs(v[pos] // unit[pos]))137 span = Fraction(1)138 for p, ex in zip(primes, unit):139 span *= Fraction(p) ** (ex * mult)140 return "lattice", (1 / span if span < 1 else span)141142# THE MASS SIDE143144def classes(w):145 out = {}146 for x in w:147 out[x] = out.get(x, 0) + 1148 return sorted(out.items(), key=lambda t: (-t[1], -t[0]))149150def log_weights(w):151 return [((Decimal(x.denominator) / Decimal(x.numerator)).ln(), g) for x, g in classes(w)]152153def dirichlet_root(lw):154 lo, hi = 1e-9, 40.0155 for _ in range(200):156 mid = 0.5 * (lo + hi)157 if sum(g * e ** (-mid * float(l)) for l, g in lw) > 1:158 lo = mid159 else:160 hi = mid161 return 0.5 * (lo + hi)162163def stopping_table(lw, amax):164 rows = []165 caps = [int(amax / float(l)) + 2 for l, _ in lw]166167 def walk(i, acc, ks):168 if i == len(lw):169 total, rest = 1, sum(ks)170 for k, (_, g) in zip(ks, lw):171 total *= comb(rest, k) * g ** k172 rest -= k173 rows.append((float(acc), total))174 return175 for k in range(caps[i] + 1):176 nxt = acc + k * lw[i][0]177 if float(nxt) > amax:178 break179 walk(i + 1, nxt, ks + [k])180181 walk(0, Decimal(0), [])182 rows.sort()183 ws, cum, acc = [], [], 0184 for x, c in rows:185 acc += c186 if ws and ws[-1] == x:187 cum[-1] = acc188 else:189 ws.append(x)190 cum.append(acc)191 return ws, cum192193def oscillation(ws, cum, delta, window, points=3000):194 out = []195 for a in MASS_A:196 vals = []197 for i in range(points):198 x = a + i * window / (points - 1)199 k = bisect.bisect_right(ws, x)200 if k:201 vals.append(log(cum[k - 1]) - delta * x)202 out.append(max(vals) - min(vals))203 return out204205# THE LENGTH SIDE206207def integer_weights(w):208 d = 1209 for x in w:210 d = d * x.denominator // gcd(d, x.denominator)211 return [int(x * d) for x in w], d212213def mass_within(level, nums):214 pts = [(0, 0, 1)]215 for _ in range(level):216 nxt = []217 for (i, j, m) in pts:218 for k, (a, b) in enumerate(CELLS):219 nxt.append((Q * i + a, Q * j + b, m * nums[k]))220 pts = nxt221 side = Q ** level222 top = isqrt(2 * side * side) + 2223 hist = [0] * (top + 2)224 for (i, j, m) in pts:225 s = i * i + j * j226 r = isqrt(s)227 if r * r < s:228 r += 1229 hist[r] += m230 out, acc = [0] * (top + 1), 0231 for r in range(top + 1):232 acc += hist[r]233 out[r] = acc234 return out235236def fold(rs, g, bins=BINS):237 b = [[] for _ in range(bins)]238 for r, v in zip(rs, g):239 u = log(r, Q)240 b[int((u - int(u)) * bins) % bins].append(v)241 return [sum(x) / len(x) if x else None for x in b]242243def bar(rs, g, cut):244 f1, f2 = fold(rs[:cut], g[:cut]), fold(rs[cut:], g[cut:])245 shared = [abs(a - b) for a, b in zip(f1, f2) if a is not None and b is not None]246 return max(shared), len(shared)247248def ripple(w):249 nums, den = integer_weights(w)250 M = mass_within(LEVEL, nums)251 prev = mass_within(LEVEL - 1, nums)252 lim = len(prev) - 1253 exact = all(M[r] == nums[0] * prev[min(r, lim)] for r in range(1, 2 * Q ** (LEVEL - 1)))254 alpha = -log(float(w[0]), Q)255 rs = list(range(Q ** FOLD_LO, Q ** FOLD_HI + 1))256 g = [log(M[r]) - alpha * log(r) for r in rs]257 f = fold(rs, g)258 f = [x - sum(f) / len(f) for x in f]259 drift, shared = bar(rs, g, len(rs) // 2)260 period, _ = bar(rs, g, Q ** ((FOLD_LO + FOLD_HI) // 2) - Q ** FOLD_LO)261 decade = []262 for k in range(FOLD_LO, FOLD_HI):263 decade.append(max(abs((log(M[Q * r]) - alpha * log(Q * r)) - (log(M[r]) - alpha * log(r)))264 for r in range(Q ** k, Q ** (k + 1))))265 return alpha, f, max(f) - min(f), drift, shared, period, decade, exact266267# SAFE ROOTS268269def floor_root(num, den, n, digits=7):270 scale = 10 ** digits271 lo = int(e ** ((log(num) - log(den)) / n) * scale) - 4272 while (lo + 1) ** n * den <= num * scale ** n:273 lo += 1274 while lo ** n * den > num * scale ** n:275 lo -= 1276 return "%d.%0*d" % (lo // scale, digits, lo % scale)277278def ceil_root(num, den, n, digits=7):279 scale = 10 ** digits280 hi = int(e ** ((log(num) - log(den)) / n) * scale) + 4281 while (hi - 1) ** n * den >= num * scale ** n:282 hi -= 1283 while hi ** n * den < num * scale ** n:284 hi += 1285 return "%d.%0*d" % (hi // scale, digits, hi % scale)286287# EXACT ARITHMETIC IN Q(sqrt ROOT)288289SQRT_PREC = 10 ** 20290SQRT_CACHE = {}291ZERO = (Fraction(0), Fraction(0))292293def rat_sqrt_down(x):294 if x <= 0:295 return Fraction(0)296 n, d = x.numerator, x.denominator297 return Fraction(isqrt(n * d * SQRT_PREC * SQRT_PREC), d * SQRT_PREC)298299def rat_sqrt_up(x):300 if x <= 0:301 return Fraction(0)302 n, d = x.numerator, x.denominator303 return Fraction(isqrt(n * d * SQRT_PREC * SQRT_PREC) + 1, d * SQRT_PREC)304305def root_bounds(root):306 if root not in SQRT_CACHE:307 SQRT_CACHE[root] = (rat_sqrt_down(Fraction(root)), rat_sqrt_up(Fraction(root)))308 return SQRT_CACHE[root]309310def sadd(x, y):311 return (x[0] + y[0], x[1] + y[1])312313def ssub(x, y):314 return (x[0] - y[0], x[1] - y[1])315316def smul(x, y, root):317 return (x[0] * y[0] + root * x[1] * y[1], x[0] * y[1] + x[1] * y[0])318319def sbounds(x, root):320 a, b = x321 if b == 0:322 return a, a323 lo, hi = root_bounds(root)324 if b > 0:325 return a + b * lo, a + b * hi326 return a + b * hi, a + b * lo327328def mat_mul(A, B, root):329 n = len(A)330 out = []331 for i in range(n):332 row = []333 for j in range(n):334 acc = ZERO335 for k in range(n):336 acc = sadd(acc, smul(A[i][k], B[k][j], root))337 row.append(acc)338 out.append(tuple(row))339 return tuple(out)340341def rho_lower(P, root):342 if len(P) == 1:343 lo, hi = sbounds(P[0][0], root)344 return lo if lo >= 0 else (-hi if hi <= 0 else Fraction(0))345 t = sadd(P[0][0], P[1][1])346 dd = ssub(smul(P[0][0], P[1][1], root), smul(P[0][1], P[1][0], root))347 disc = ssub(smul(t, t, root), (4 * dd[0], 4 * dd[1]))348 tlo, thi = sbounds(t, root)349 tabs = tlo if tlo >= 0 else (-thi if thi <= 0 else Fraction(0))350 dlo, dhi = sbounds(disc, root)351 if dlo >= 0:352 return (tabs + rat_sqrt_down(dlo)) / 2353 if dhi <= 0:354 return rat_sqrt_down(max(Fraction(0), sbounds(dd, root)[0]))355 return tabs / 2356357def norm_upper(P, root):358 if len(P) == 1:359 lo, hi = sbounds(P[0][0], root)360 return max(abs(lo), abs(hi))361 G = [[ZERO, ZERO], [ZERO, ZERO]]362 for i in range(2):363 for j in range(2):364 acc = ZERO365 for k in range(2):366 acc = sadd(acc, smul(P[k][i], P[k][j], root))367 G[i][j] = acc368 t = sadd(G[0][0], G[1][1])369 dd = ssub(smul(G[0][0], G[1][1], root), smul(G[0][1], G[1][0], root))370 thi = sbounds(t, root)[1]371 dlo = max(Fraction(0), sbounds(dd, root)[0])372 disc = max(Fraction(0), thi * thi - 4 * dlo)373 return rat_sqrt_up((thi + rat_sqrt_up(disc)) / 2)374375def gripenberg(mats, root, maxlen):376 words = [((k,), m) for k, m in enumerate(mats)]377 best, arg = "0.0000000", ()378 for n in range(1, maxlen + 1):379 for word, P in words:380 r = rho_lower(P, root)381 if r <= 0:382 continue383 s = floor_root(r.numerator, r.denominator, n)384 if float(s) > float(best):385 best, arg = s, word386 if n == maxlen:387 break388 words = [(word + (k,), mat_mul(m, P, root)) for word, P in words for k, m in enumerate(mats)]389 return best, arg390391def norm_bracket(mats, root, maxlen):392 prods, best, arg = list(mats), None, 0393 for n in range(1, maxlen + 1):394 u = max(norm_upper(P, root) for P in prods)395 s = ceil_root(u.numerator, u.denominator, n)396 if best is None or float(s) < float(best):397 best, arg = s, n398 if n == maxlen:399 break400 prods = [mat_mul(m, P, root) for P in prods for m in mats]401 return best, arg402403def dl_matrices(c):404 N = len(c) - 1405 cc = lambda k: c[k] if 0 <= k <= N else ZERO406 T0 = tuple(tuple(cc(2 * i - j - 1) for j in range(1, N + 1)) for i in range(1, N + 1))407 T1 = tuple(tuple(cc(2 * i - j) for j in range(1, N + 1)) for i in range(1, N + 1))408 return T0, T1409410def restrict(T, root):411 N = len(T)412 TB = [[ssub(T[i][j], T[i][j + 1]) for j in range(N - 1)] for i in range(N)]413 C = []414 acc = [ZERO] * (N - 1)415 for i in range(N - 1):416 acc = [sadd(acc[j], TB[i][j]) for j in range(N - 1)]417 C.append(tuple(acc))418 ok = all(sadd(acc[j], TB[N - 1][j]) == ZERO for j in range(N - 1))419 return tuple(C), ok420421# THE MULTIFRACTAL SPECTRUM422423def partition(w, s):424 return sum(float(x) ** s for x in w)425426def tau(w, s):427 return log(partition(w, s), Q)428429def alpha_of(w, s):430 z = partition(w, s)431 return -sum(float(x) ** s * log(float(x)) for x in w) / (z * log(Q))432433def legendre(w, a):434 lo, hi = -300.0, 300.0435 for _ in range(300):436 mid = 0.5 * (lo + hi)437 if alpha_of(w, mid) > a:438 lo = mid439 else:440 hi = mid441 s = 0.5 * (lo + hi)442 return tau(w, s) + s * a443444def coarse(w, level):445 cls = classes(w)446 out = {}447448 def walk(i, rest, ks):449 if i == len(cls) - 1:450 ks = ks + [rest]451 mass, cnt, left = Fraction(1), 1, level452 for k, (x, g) in zip(ks, cls):453 mass *= x ** k454 cnt *= comb(left, k) * g ** k455 left -= k456 hit = out.setdefault(mass, [0, []])457 hit[0] += cnt458 hit[1].append(tuple(ks))459 return460 for k in range(rest + 1):461 walk(i + 1, rest - k, ks + [k])462463 walk(0, level, [])464 return cls, out465466def lnfac(n):467 if n < 1:468 return 0.0469 return 0.5 * log(2 * pi * n) + n * log(n) - n + 1.0 / (12 * n) - 1.0 / (360 * n ** 3)470471def stirling_count(level, ks, cls):472 v = lnfac(level)473 for k, (_, g) in zip(ks, cls):474 v += -lnfac(k) + k * log(g)475 return v476477# THE RUN478479print("WEIGHTED DESIGNS")480print()481print("== THE OBJECT ==")482print("base q %d, dimension D %d, cells %s, level %d" % (Q, D, " ".join(str(c) for c in CELLS), LEVEL))483print("weights %s, a probability vector over a common denominator" % ", ".join(str(x) for x in WEIGHT))484nums, den = integer_weights(WEIGHT)485print("integer masses %s over %d, so every level-%d cell mass is an integer over %d^%d"486 % (nums, den, LEVEL, den, LEVEL))487check("weights sum to one", sum(WEIGHT), Fraction(1))488print()489490print("== THE ARITHMETIC CLASS ==")491print("positive rational weights only: a weight has a prime-exponent vector only there, and the logs of the primes are Q-independent")492print("under that hypothesis the group generated by the log-weights is cyclic iff the prime-exponent matrix has rank 1")493CLASSOF = {}494for name, w in RENEWAL:495 check("%s is a vector of rationals" % name, all(isinstance(x, Fraction) for x in w), True)496 primes, rows = exponent_rows(w)497 kind, span = arithmetic_class(w)498 CLASSOF[name] = kind499 print(" %-12s primes %-10s rows %-28s rank %d %s%s"500 % (name, primes, rows, rank(rows), kind,501 "" if span is None else ", span ln %s" % span))502print()503504print("== THE MASS SIDE ==")505print("delta is the root of sum_f w_f^s = 1; N(a) counts words of mass at least e^(-a)")506OSC = {}507for name, w in RENEWAL:508 lw = log_weights(w)509 delta = dirichlet_root(lw)510 ws, cum = stopping_table(lw, MASS_MAX)511 osc = oscillation(ws, cum, delta, log(Q))512 OSC[name] = osc513 spread = (max(osc) - min(osc)) / max(osc)514 falling = all(x > y for x, y in zip(osc, osc[1:])) and osc[-1] / osc[0] < 0.8515 trend = "flat" if spread < 0.01 else ("decaying" if falling else "mixed")516 want = "flat" if CLASSOF[name] == "lattice" else "decaying"517 check("trend %s" % name, trend, want)518 check("delta %s" % name, "%.10f" % delta, NAPKIN_MASS[name][0])519 check("oscillation %s" % name, " ".join("%.5f" % v for v in osc), NAPKIN_MASS[name][1])520 print(" %-12s delta %.10f %-11s osc of ln N(a) - delta a over a window of ln %d at a = %s: %s %s"521 % (name, delta, CLASSOF[name], Q, ",".join(str(a) for a in MASS_A),522 " ".join("%.5f" % v for v in osc), trend))523print("for a probability vector s = 1 solves the Dirichlet equation, so delta = 1 identically;")524print("what weights move on the mass side is the arithmetic class alone")525for name, w in PROBS:526 check("delta one %s" % name, "%.10f" % dirichlet_root(log_weights(w)), "1.0000000000")527print()528529print("== THE LENGTH SIDE ==")530print("M(r) is the mass of cells within radius r of the corner fixed point, detrended by")531print("alpha = -log_q w_0 and folded into %d bins of log_%d r over whole periods of R in %d^%d..%d^%d"532 % (BINS, Q, Q, FOLD_LO, Q, FOLD_HI))533BASE, RIP = None, {}534for name, w in PROBS:535 alpha, f, swing, drift, shared, period, decade, exact = ripple(w)536 if BASE is None:537 BASE, gap = f, 0.0538 else:539 gap = max(abs(a - b) for a, b in zip(f, BASE))540 RIP[name] = (swing, drift, gap)541 check("swing %s" % name, "%.6f" % swing, NAPKIN_LENGTH[name][0])542 check("drift %s" % name, "%.6f" % drift, NAPKIN_LENGTH[name][1])543 check("0/1 gap %s" % name, "%.6f" % gap, NAPKIN_LENGTH[name][2])544 check("decade %s" % name, " ".join("%.4f" % v for v in decade), NAPKIN_LENGTH[name][3])545 check("ripple kept %s" % name, swing > 10 * drift, True)546 check("ripple kept on whole periods %s" % name, swing > 10 * period, True)547 check("cross level %s" % name, exact, True)548 print(" %-12s alpha %.9f swing %.6f on drift bar %.6f gap against equal weights %.6f"549 % (name, alpha, swing, drift, gap))550 print(" %-12s log_%d periodicity residual by decade: %s cross-level identity exact %s"551 % ("", Q, " ".join("%.4f" % v for v in decade), exact))552 print(" %-12s drift bar on %d shared bins, %.6f on a whole-period split of the window"553 % ("", shared, period))554print("equal weights give every cell mass 1, so M(r) is the cell count and the fold is the 0/1 ripple")555print()556557print("== THE LOCAL DIMENSIONS ==")558for name, w in PROBS:559 amin = -log(float(max(w)), Q)560 amax = -log(float(min(w)), Q)561 ainf = -sum(float(x) * log(float(x)) for x in w) / log(Q)562 print(" %-12s alpha_min %.9f alpha_max %.9f alpha at s = 1 %.9f" % (name, amin, amax, ainf))563print()564565print("== THE JSR LADDER ==")566print("rung 0, the designs: supp phi sits in the unit cell, so the Daubechies-Lagarias matrices")567print("are 1 x 1, T_f = [q^D w_f], and the joint spectral radius is q^D max_f w_f in closed form")568for name, w in PROBS:569 T = [Q ** D * x for x in w]570 jsr = max(T)571 mats = [((( x, Fraction(0)),),) for x in T]572 lo, word = gripenberg(mats, 0, GRIPENBERG_LEN)573 hi, at = norm_bracket(mats, 0, GRIPENBERG_LEN)574 holder = -D - log(float(max(w)), Q)575 print(" %-12s T_f = %-28s JSR %-6s solver [%s, %s] alpha_Holder(phi) %.9f"576 % (name, "[" + ", ".join(str(x) for x in T) + "]", str(jsr), lo, hi, holder))577 check("rung 0 matrices %s" % name, "[" + ", ".join(str(x) for x in T) + "]", NAPKIN_RUNG0[name][0])578 check("rung 0 jsr %s" % name, str(jsr), NAPKIN_RUNG0[name][1])579 check("rung 0 alpha_min %s" % name, "%.9f" % -log(float(max(w)), Q), NAPKIN_RUNG0[name][2])580 check("rung 0 alpha_max %s" % name, "%.9f" % -log(float(min(w)), Q), NAPKIN_RUNG0[name][3])581 check("rung 0 Holder %s" % name, "%.9f" % holder, NAPKIN_RUNG0[name][4])582 check("rung 0 lower %s" % name, float(lo) <= float(jsr), True)583 check("rung 0 upper %s" % name, float(hi) >= float(jsr), True)584 check("rung 0 collapse %s" % name, float(hi) - float(lo) < 1e-6, True)585print()586print("rung 1, the first overlap at q = 2: the same solver on T_0 = (c_(2i-j-1)), T_1 = (c_(2i-j)),")587print("restricted to sum v_i = 0, Gripenberg words to length %d and norm upper bounds to length %d"588 % (GRIPENBERG_LEN, NORM_LEN))589MASKS = (590 ("hat (1/2, 1, 1/2)", 0,591 ((Fraction(1, 2), Fraction(0)), (Fraction(1), Fraction(0)), (Fraction(1, 2), Fraction(0)))),592 ("D4", 3,593 ((Fraction(1, 4), Fraction(1, 4)), (Fraction(3, 4), Fraction(1, 4)),594 (Fraction(3, 4), Fraction(-1, 4)), (Fraction(1, 4), Fraction(-1, 4)))),595)596SQ3 = (Fraction(1, 4), Fraction(1, 4))597D4_JSR = "0.6830127"598D4_ALPHA = "0.5500157"599print("the hat mask carries alpha = 1 exactly; D4's lower end is the closed form (1 + sqrt 3)/4 and")600print("its Holder exponent 2 - log_2(1 + sqrt 3) = %s must sit inside the printed bracket" % D4_ALPHA)601for name, root, c in MASKS:602 even = ZERO603 odd = ZERO604 for i, x in enumerate(c):605 if i % 2 == 0:606 even = sadd(even, x)607 else:608 odd = sadd(odd, x)609 check("mask %s even" % name, even, (Fraction(1), Fraction(0)))610 check("mask %s odd" % name, odd, (Fraction(1), Fraction(0)))611 T0, T1 = dl_matrices(c)612 V0, ok0 = restrict(T0, root)613 V1, ok1 = restrict(T1, root)614 check("invariant subspace %s" % name, ok0 and ok1, True)615 lo, word = gripenberg([V0, V1], root, GRIPENBERG_LEN)616 hi, at = norm_bracket([V0, V1], root, NORM_LEN)617 alo = safe_log2(Fraction(hi), 7, ROUND_FLOOR)618 ahi = safe_log2(Fraction(lo), 7, ROUND_CEILING)619 check("rung 1 lower %s" % name, lo, NAPKIN_RUNG1[name][0])620 check("rung 1 upper %s" % name, hi, NAPKIN_RUNG1[name][1])621 check("rung 1 alpha lower %s" % name, alo, NAPKIN_RUNG1[name][2])622 check("rung 1 alpha upper %s" % name, ahi, NAPKIN_RUNG1[name][3])623 check("rung 1 bracket order %s" % name, float(lo) <= float(hi), True)624 print(" %-18s dim %d JSR in [%s, %s], the lower at word %s, the upper at length %d"625 % (name, len(V0), lo, hi, "".join(str(k) for k in word), at))626 print(" %-18s alpha = -log_2 JSR in [%s, %s]" % ("", alo, ahi))627 if name == "D4":628 check("D4 lower is (1 + sqrt 3)/4", lo,629 floor_root(sbounds(SQ3, 3)[0].numerator, sbounds(SQ3, 3)[0].denominator, 1))630 check("D4 lower reads", lo, D4_JSR)631 closed = safe_log2(sbounds(SQ3, 3)[1], 7, ROUND_CEILING)632 check("D4 closed form", closed, D4_ALPHA)633 check("D4 closed form inside the bracket", float(alo) <= float(closed) <= float(ahi), True)634print()635636print("== THE CONTROLS ==")637print("pre-registered: equal weights reproduce the 0/1 design exactly; log-commensurable weights")638print("keep every ripple; only an irrational log ratio may move a mass observable")639print(" %-12s %-11s %-9s %-9s %-8s %-9s %-9s %s"640 % ("weights", "class", "swing", "drift", "0/1 gap", "osc first", "osc last", "verdict"))641moved = 0642for name in CONTROLS:643 swing, drift, gap = RIP[name]644 osc = OSC[name]645 kept = swing > 10 * drift646 flat = osc[-1] / osc[0] > 0.9647 if not flat:648 moved += 1649 verdict = "length kept, mass %s" % ("flat" if flat else "moved")650 print(" %-12s %-11s %-9.6f %-9.6f %-8.6f %-9.5f %-9.5f %s"651 % (name, CLASSOF[name], swing, drift, gap, osc[0], osc[-1], verdict))652 check("control ripple %s" % name, kept, True)653check("equal weights reproduce the 0/1 fold", "%.6f" % RIP["1/3,1/3,1/3"][2], "0.000000")654check("exactly one control moves the mass side", moved, 1)655print("length observable moves under none of the three; mass observable moves under one of the three")656print()657658print("== THE MULTIFRACTAL SPECTRUM ==")659print("Moran self-similar measure, ratios all 1/q, open set condition: tau(s) solves")660print("sum_f w_f^s (1/q)^tau(s) = 1, hence tau(s) = log_q sum_f w_f^s, and f(alpha) = inf_s (alpha s + tau(s))")661print(" s sum_f w_f^s tau(s) alpha(s) f(alpha(s))")662for s in (-2, -1, 0, 1, 2, 3):663 zr = sum(x ** s for x in WEIGHT)664 a = alpha_of(WEIGHT, s)665 print(" %2d %-18s %-13.9f %-13.9f %.9f" % (s, str(zr), tau(WEIGHT, s), a, tau(WEIGHT, s) + s * a))666check("tau at 0 is the box dimension", "%.12f" % tau(WEIGHT, 0), "%.12f" % log(len(WEIGHT), Q))667check("tau at 1 vanishes", "%.12f" % tau(WEIGHT, 1), "%.12f" % 0.0)668a1 = alpha_of(WEIGHT, 1)669check("f is tangent to the diagonal at s = 1", "%.9f" % legendre(WEIGHT, a1), "%.9f" % a1)670amin = -log(float(max(WEIGHT)), Q)671amax = -log(float(min(WEIGHT)), Q)672print("alpha ranges over [%.9f, %.9f]; f(alpha_min) = log_q #argmax w = %.9f, f(alpha_max) = %.9f"673 % (amin, amax, log(sum(1 for x in WEIGHT if x == max(WEIGHT)), Q),674 log(sum(1 for x in WEIGHT if x == min(WEIGHT)), Q)))675print()676print("the level-L box partition carries the moments exactly: sum_i mu_i^s = (sum_f w_f^s)^L")677for level in SPECTRUM_LEVELS:678 cls, groups = coarse(WEIGHT, level)679 for s in (-2, -1, 0, 1, 2, 3):680 got = sum(cnt * mass ** s for mass, (cnt, _) in groups.items())681 check("partition identity L %d s %d" % (level, s), got, sum(x ** s for x in WEIGHT) ** level)682 print(" level %2d %d distinct box masses, %d boxes, moments exact at s = -2..3"683 % (level, len(groups), sum(c for c, _ in groups.values())))684print()685print("the coarse-grained spectrum f_L(alpha) = log_q N(alpha) / L against the Legendre transform")686print("f_L <= f at every level and every achievable alpha, since N_i mu_i^s <= q^(L tau(s)); the gap")687print("is the Stirling volume term of the multinomial count and falls like log L / L")688print(" level alpha f_L(alpha) f(alpha) deficit max deficit max |f_L - Stirling|")689DEF = {}690for level in SPECTRUM_LEVELS:691 cls, groups = coarse(WEIGHT, level)692 worst, stir, mid = 0.0, 0.0, None693 for mass, (cnt, ks) in sorted(groups.items()):694 a = -log(float(mass)) / (level * log(Q))695 fl = log(cnt) / (level * log(Q))696 fx = legendre(WEIGHT, a)697 check("f_L under the transform L %d" % level, fl <= fx + 1e-12, True)698 worst = max(worst, fx - fl)699 if len(ks) == 1:700 stir = max(stir, abs(fl - stirling_count(level, ks[0], cls) / (level * log(Q))))701 if ks[0][0] * 2 == level:702 mid = (a, fl, fx)703 if len(ks) == 1 and 0 in ks[0]:704 check("endpoint exact L %d" % level, "%.9f" % abs(fx - fl), "0.000000000")705 DEF[level] = worst706 print(" %5d %-13.9f %-13.9f %-13.9f %-10.6f %-12.6f %.6f"707 % (level, mid[0], mid[1], mid[2], mid[2] - mid[1], worst, stir))708 check("Stirling agreement L %d" % level, stir < 1e-3, True)709check("the deficit falls with the level", DEF[6] > DEF[8] > DEF[10], True)710print()711print("the whole coarse-grained band at level %d, printed row by row" % SPECTRUM_LEVELS[-1])712print(" k box mass boxes alpha f_L(alpha) f(alpha)")713cls, groups = coarse(WEIGHT, SPECTRUM_LEVELS[-1])714level = SPECTRUM_LEVELS[-1]715for mass, (cnt, ks) in sorted(groups.items(), reverse=True):716 a = -log(float(mass)) / (level * log(Q))717 print(" %2d %-25s %-8d %-13.9f %-13.9f %.9f"718 % (ks[0][0], str(mass), cnt, a, log(cnt) / (level * log(Q)), legendre(WEIGHT, a)))719top = max(log(c) / (level * log(Q)) for c, _ in groups.values())720print("the band's top %.9f rises to tau(0) = log_q |F| = %.9f as the level grows"721 % (top, tau(WEIGHT, 0)))722check("the band's top sits under the box dimension", top < tau(WEIGHT, 0), True)723print()724725print("== THE VERDICT ==")726if FAILURES:727 for line in FAILURES:728 print("FAIL " + line)729 raise SystemExit(1)730print("every assertion holds")