complex_dimensions.py
14.1 kB · python · 378 lines
1import math2import time3from fractions import Fraction45import numpy as np6from scipy.special import gammaln78LATTICE = [(3, (0, 2)), (5, (0, 2, 4)), (15, (0, 4, 10, 14)), (15, (0, 2, 4, 10, 12, 14))]9RATIOS = (1 / 3, 1 / 5)10CONTROL_MAPS = ((1 / 3, 0.0), (1 / 5, 4 / 5))11CONTROL_SEED = Fraction(7, 15)12APERIODIC = ((3, (0, 2)), (5, (0, 4)))13POLE_RANGE = 4014REMAX, IMMAX = 3.0, 40.015BOUNDARY = 20000016UMIN, UMAX, NU = math.log(1 / 0.03), math.log(1e6), 300017SAFETY = 50.018OMEGA_GRID = np.linspace(0.5, 30.0, 24001)19NBINS = 4020SAMPLES = 600021WINDOWS = [6.0, 15.0, 24.0, 33.0, 42.0, 51.0, 60.0]22TAIL = 16023LN3, LN5 = math.log(3), math.log(5)242526def label(base, digits):27 return f"base {base}, digits {{{','.join(str(d) for d in digits)}}}"282930def gap_seeds(base, digits):31 ds = sorted(digits)32 return [Fraction(b - a - 1, base) for a, b in zip(ds, ds[1:]) if b - a > 1]333435def dimension(base, digits):36 return math.log(len(digits)) / math.log(base)373839def solve_moran(ratios):40 d = 0.541 for _ in range(200):42 f = sum(r**d for r in ratios) - 143 fp = sum(math.log(r) * r**d for r in ratios)44 d -= f / fp45 return d464748def lattice_poles(base, digits):49 k = len(digits)50 d = dimension(base, digits)51 omega = 2 * math.pi / math.log(base)52 s = d + 1j * omega * np.arange(-POLE_RANGE, POLE_RANGE + 1)53 denominator = np.abs(1 - k * np.exp(-s * math.log(base)))54 seeds = np.array([float(g) for g in gap_seeds(base, digits)])55 numerator = np.abs(np.exp(np.outer(s, np.log(seeds))).sum(axis=1))56 d1 = sum(gap_seeds(base, digits))57 print(f" {label(base, digits)}")58 print(f" d = {d:.6f} ({d:.9f}) omega = 2*pi/ln({base}) = {omega:.6f} ({omega:.9f})")59 print(f" {s.size} predicted poles m = -{POLE_RANGE}..{POLE_RANGE}: max |1 - k n^(-s)| = {denominator.max():.1e}")60 print(f" min |D(s)| = {numerator.min():.6f} D(1) = {d1} = 1 - k/n: {d1 == 1 - Fraction(k, base)}")616263def moran(s):64 return 1 - sum(np.exp(s * math.log(r)) for r in RATIOS)656667def moran_prime(s):68 return -sum(math.log(r) * np.exp(s * math.log(r)) for r in RATIOS)697071def winding_number():72 corners = [-REMAX - 1j * IMMAX, REMAX - 1j * IMMAX, REMAX + 1j * IMMAX, -REMAX + 1j * IMMAX]73 steps = np.arange(BOUNDARY) / BOUNDARY74 path = np.concatenate([a + (b - a) * steps for a, b in zip(corners, corners[1:] + corners[:1])])75 phase = np.angle(moran(path))76 turns = np.diff(np.concatenate([phase, phase[:1]]))77 turns = (turns + np.pi) % (2 * np.pi) - np.pi78 return int(round(turns.sum() / (2 * np.pi)))798081def control_poles():82 re = np.linspace(-REMAX, REMAX, 61)83 im = np.linspace(-IMMAX, IMMAX, 121)84 s = (re[:, None] + 1j * im[None, :]).ravel()85 alive = np.ones(s.size, dtype=bool)86 with np.errstate(all="ignore"):87 for _ in range(200):88 s = s - moran(s) / moran_prime(s)89 alive &= np.isfinite(s) & (np.abs(s.real) <= 40) & (np.abs(s.imag) <= 4000)90 s = np.where(alive, s, np.nan)91 good = alive & (np.abs(moran(s)) < 1e-13) & (np.abs(s.imag) <= IMMAX) & (np.abs(s.real) <= REMAX)92 roots = []93 for z in sorted(s[good], key=lambda z: (z.imag, z.real)):94 if all(abs(z - q) > 1e-7 for q in roots):95 roots.append(z)96 roots = np.array(roots)97 gaps = np.diff(roots.imag)98 d = solve_moran(RATIOS)99 print(f" two-ratio control, ratios 1/3 and 1/5, d = {d:.9f} solving 3^(-d) + 5^(-d) = 1")100 print(f" roots in Re [-{REMAX:.0f},{REMAX:.0f}], Im [-{IMMAX:.0f},{IMMAX:.0f}]: {roots.size} winding number of 1 - 3^(-s) - 5^(-s) over the box: {winding_number()}")101 print(f" Re range {roots.real.min():.6f} .. {roots.real.max():.6f} Im gaps {gaps.min():.6f} .. {gaps.max():.6f}")102 for base in (3, 5, 15):103 omega = 2 * math.pi / math.log(base)104 ratio = roots.imag / omega105 print(f" worst offset from spacing 2*pi/ln{base}: {np.abs(ratio - np.round(ratio)).max():.2f} of a step")106 print(f" ln3/ln5 = {LN3 / LN5:.12f}")107108109def compose(schedule):110 base, digits = 1, [0]111 for b, ds in schedule:112 digits = [x * b + d for x in digits for d in ds]113 base *= b114 return base, sorted(digits)115116117def integer_cover(schedule, levels):118 lefts, den = [0], 1119 for i in range(levels):120 b, ds = schedule[i % len(schedule)]121 lefts = [x * b + d for x in lefts for d in ds]122 den *= b123 return sorted(Fraction(x, den) for x in lefts)124125126def composition():127 for schedule in [((3, (0, 2)), (5, (0, 4))), ((3, (0, 2)), (5, (0, 2, 4)))]:128 base, digits = compose(schedule)129 print(f" {' then '.join(label(b, ds) for b, ds in schedule)} -> {label(base, digits)}")130 alternating = integer_cover(((3, (0, 2)), (5, (0, 4))), 8)131 product = integer_cover(((15, (0, 4, 10, 14)),), 4)132 print(f" 8 alternating levels and 4 base-15 levels: {len(alternating)} and {len(product)} intervals, identical as exact fractions: {alternating == product}")133134135def level_cover(base, digits, floor):136 lefts = np.zeros(1)137 length = 1.0138 while length > floor:139 length /= base140 lefts = (lefts[:, None] + np.array(digits, dtype=float) * length).ravel()141 return np.sort(lefts), np.full(lefts.size, length)142143144def schedule_cover(schedule, floor):145 lefts = np.zeros(1)146 length = 1.0147 used = 0148 lnk, lnn = 0.0, 0.0149 while length > floor:150 base, digits = schedule[used % len(schedule)]151 used += 1152 length /= base153 lnk += math.log(len(digits))154 lnn += math.log(base)155 lefts = (lefts[:, None] + np.array(digits, dtype=float) * length).ravel()156 return np.sort(lefts), np.full(lefts.size, length), used, lnk / lnn157158159def ratio_cover(maps, floor):160 lefts, lens = np.zeros(1), np.ones(1)161 done_l, done_s = [], []162 while lefts.size:163 small = lens <= floor164 done_l.append(lefts[small])165 done_s.append(lens[small])166 lefts, lens = lefts[~small], lens[~small]167 if not lefts.size:168 break169 lefts = np.concatenate([lefts + c * lens for _, c in maps])170 lens = np.concatenate([r * lens for r, _ in maps])171 lefts, lens = np.concatenate(done_l), np.concatenate(done_s)172 order = np.argsort(lefts)173 return lefts[order], lens[order]174175176def thue_morse(n):177 return [bin(i).count("1") % 2 for i in range(n)]178179180def box_count(lefts, rights, eps):181 lo = np.floor(lefts / eps).astype(np.int64)182 hi = np.floor(rights / eps).astype(np.int64)183 return int((hi - lo + 1).sum() - np.count_nonzero(hi[:-1] == lo[1:]))184185186def detrended(lefts, rights, d):187 u = np.linspace(UMIN, UMAX, NU)188 g = np.array([math.log(box_count(lefts, rights, math.exp(-x))) for x in u]) - d * u189 slope, intercept = np.polyfit(u, g, 1)190 return u, g - (slope * u + intercept)191192193def power_at(w, dt, y):194 return abs(np.exp(-1j * w * dt) @ y) ** 2195196197def golden(f, lo, hi):198 phi = (math.sqrt(5) - 1) / 2199 a, b = lo, hi200 c, d = b - phi * (b - a), a + phi * (b - a)201 fc, fd = f(c), f(d)202 for _ in range(80):203 if fc < fd:204 b, d, fd = d, c, fc205 c = b - phi * (b - a)206 fc = f(c)207 else:208 a, c, fc = c, d, fd209 d = a + phi * (b - a)210 fd = f(d)211 return (a + b) / 2212213214def periodogram_peak(u, g):215 j = np.arange(g.size)216 window = 0.42 - 0.5 * np.cos(2 * np.pi * j / (g.size - 1)) + 0.08 * np.cos(4 * np.pi * j / (g.size - 1))217 y = g * window218 y = y - y.mean()219 dt = u - u[0]220 power = np.concatenate([np.abs(np.exp(-1j * np.outer(chunk, dt)) @ y) ** 2 for chunk in np.array_split(OMEGA_GRID, 48)])221 i = int(np.argmax(power))222 lo, hi = OMEGA_GRID[max(i - 1, 0)], OMEGA_GRID[min(i + 1, OMEGA_GRID.size - 1)]223 return golden(lambda w: -power_at(w, dt, y), lo, hi)224225226def folded_variance(u, g, period, nbins):227 phase = ((u - u[0]) / period) % 1.0228 index = np.minimum((phase * nbins).astype(int), nbins - 1)229 sums = np.bincount(index, g, nbins)230 counts = np.bincount(index, minlength=nbins)231 profile = np.where(counts > 0, sums / np.maximum(counts, 1), 0.0)232 return 1 - (g - profile[index]).var() / g.var()233234235def box_report(name, lefts, lens, d, base):236 u, g = detrended(lefts, lefts + lens, d)237 omega = periodogram_peak(u, g)238 fold = [folded_variance(u, g, math.log(p), NBINS) for p in (3, 5, 15)]239 print(f" {name}: {lefts.size} intervals, d = {d:.6f}")240 if base:241 predicted = 2 * math.pi / math.log(base)242 err = abs(omega - predicted) / predicted * 100243 print(f" predicted omega {predicted:.4f} periodogram peak {omega:.4f} error {err:.3f}% within 1%: {err < 1}")244 else:245 print(f" no lattice prediction periodogram peak {omega:.4f}")246 print(f" folding variance explained at ln3 {fold[0]:.3f} ln5 {fold[1]:.3f} ln15 {fold[2]:.3f}")247248249def box_counts():250 floor = math.exp(-UMAX) / SAFETY251 print(f" u in [{UMIN:.4f},{UMAX:.4f}], {NU} points, covers refined below eps/{SAFETY:.0f}, {NBINS} folding bins")252 for base, digits in LATTICE:253 lefts, lens = level_cover(base, digits, floor)254 box_report(label(base, digits), lefts, lens, dimension(base, digits), base)255 lefts, lens = ratio_cover(CONTROL_MAPS, floor)256 box_report("two-ratio control x/3, (x+4)/5", lefts, lens, solve_moran(RATIOS), None)257 schedule = [APERIODIC[bit] for bit in thue_morse(64)]258 lefts, lens, used, d = schedule_cover(schedule, floor)259 box_report(f"aperiodic control, Thue-Morse schedule of {label(3, (0, 2))} and {label(5, (0, 4))}, {used} levels", lefts, lens, d, None)260261262def tube_digit(base, k, seeds, t):263 total = 0.0264 for g in seeds:265 a0 = 0 if g <= t else math.ceil(math.log(g / t) / math.log(base))266 total += t * (k**a0 - 1) / (k - 1) + g * (k / base) ** a0 * base / (base - k)267 return total268269270def brute_digit(base, digits, level, t):271 lefts = np.zeros(1, dtype=np.int64)272 for _ in range(level):273 lefts = (lefts[:, None] * base + np.array(sorted(digits))).ravel()274 den = base**level275 gaps = [Fraction(int(b - a - 1), den) for a, b in zip(lefts, lefts[1:]) if b - a > 1]276 cut = Fraction(t)277 covered = sum(g if g <= cut else cut for g in gaps)278 return float(covered), lefts.size / den279280281A, B = np.meshgrid(np.arange(TAIL), np.arange(TAIL), indexing="ij")282LOG_COUNT = gammaln(A + B + 1) - gammaln(A + 1) - gammaln(B + 1)283LOG_WEIGHT = LOG_COUNT - A * LN3 - B * LN5284LOG_SCALE = A * LN3 + B * LN5285286287def tube_control(t):288 g = float(CONTROL_SEED)289 small = LOG_SCALE >= math.log(g / t)290 return t * np.exp(LOG_COUNT[~small]).sum() + g * np.exp(LOG_WEIGHT[small]).sum()291292293def content_curve(measure, d, ulo, uhi):294 u = np.linspace(ulo, uhi, SAMPLES)295 return u, np.array([measure(2 * math.exp(-x)) * math.exp(-x) ** (d - 1) for x in u])296297298def window_swings(measure, d):299 u, m = content_curve(measure, d, WINDOWS[0], WINDOWS[-1])300 swings = []301 for lo, hi in zip(WINDOWS, WINDOWS[1:]):302 vals = m[(u >= lo) & (u <= hi)]303 swings.append((vals.max() - vals.min()) / vals.mean() * 100)304 return swings305306307def periodicity_defect(measure, d, period):308 worst = 0.0309 for u in np.linspace(40.0, 50.0, 600):310 a = measure(2 * math.exp(-u)) * math.exp(-u) ** (d - 1)311 b = measure(2 * math.exp(-u - period)) * math.exp(-u - period) ** (d - 1)312 worst = max(worst, abs(a - b) / abs(a))313 return worst314315316def swing_line(swings):317 return " swing of M(eps) per window u in " + ", ".join(f"[{lo:.0f},{hi:.0f}] {s:.4f}%" for (lo, hi), s in zip(zip(WINDOWS, WINDOWS[1:]), swings))318319320def tube_lattice(base, digits):321 k = len(digits)322 d = dimension(base, digits)323 seeds = [float(g) for g in gap_seeds(base, digits)]324 measure = lambda t: tube_digit(base, k, seeds, t)325 print(f" {label(base, digits)}: total gap length {sum(seeds) * base / (base - k):.12f}")326 brute, tail = brute_digit(base, digits, 6, 0.01)327 closed = measure(0.01)328 print(f" level 6, t = 0.01: brute {brute:.9f} closed form {closed:.9f} excess {closed - brute:.6e} <= cover length {tail:.6e}: {-1e-12 <= closed - brute <= tail * (1 + 1e-9) + 1e-12}")329 print(swing_line(window_swings(measure, d)))330 defects = [periodicity_defect(measure, d, math.log(p)) for p in (3, 5, 15)]331 print(f" max relative defect of M(eps) = M(eps/p) on u in [40,50]: p=3 {defects[0]:.1e} p=5 {defects[1]:.1e} p=15 {defects[2]:.1e}")332333334def tube_two_ratio():335 d = solve_moran(RATIOS)336 total = CONTROL_SEED / (1 - sum(Fraction(1, int(round(1 / r))) for r in RATIOS))337 print(f" two-ratio control, gap seed {CONTROL_SEED}, d = {d:.9f}, total gap length {total}, eps at u = 60: {math.exp(-60):.1e}")338 swings = window_swings(tube_control, d)339 print(swing_line(swings))340 print(f" swing decays {swings[0]:.2f}% -> {swings[-1]:.2f}%, monotone: {all(a > b for a, b in zip(swings, swings[1:]))}")341342343def cantor_profile():344 d = dimension(3, (0, 2))345 profile = lambda t: 2 ** (1 - d) * (t ** (d - 1) + t**d)346 tstar = (1 - d) / d347 lo, hi = profile(tstar), profile(1.0)348 seeds = [float(g) for g in gap_seeds(3, (0, 2))]349 u, m = content_curve(lambda t: tube_digit(3, 2, seeds, t), d, 50.0, 60.0)350 print(f" Cantor limit profile 2^(1-d) (t^(d-1) + t^d), t in [1/3, 1)")351 print(f" minimum {lo:.9f} at t* = (1-d)/d = {tstar:.6f} maximum {hi:.9f} at the ends swing {(hi - lo) / lo * 100:.2f}%")352 print(f" measured on u in [50,60]: min {m.min():.9f} max {m.max():.9f} min gap {abs(m.min() - lo):.1e} max gap {abs(m.max() - hi):.1e}")353354355def main():356 start = time.time()357 print("POLES: zeros of 1 - k n^(-s) at s = d + 2 pi i m/ln(n), the numerator D(s) there, and the two-ratio Moran roots")358 for base, digits in LATTICE:359 lattice_poles(base, digits)360 control_poles()361 print()362 print("COMPOSITION: alternating one base per level multiplies into the product base")363 composition()364 print()365 print("BOX COUNT: g(u) = ln N(e^-u) - d u, Blackman periodogram on a direct DFT grid, and period folding")366 box_counts()367 print()368 print("TUBE: V(eps) = sum over gaps of min(gap, 2 eps) in closed form, M(eps) = eps^(d-1) V(eps)")369 for base, digits in LATTICE:370 tube_lattice(base, digits)371 tube_two_ratio()372 cantor_profile()373 print()374 print(f"wall {time.time() - start:.1f} s")375376377if __name__ == "__main__":378 main()