unequal_split.py

1import math2import subprocess3import sys4import time5from bisect import bisect_right6from fractions import Fraction7from math import comb89import numpy as np1011K = 512HALF = Fraction(1, 2)13THIRD = Fraction(1, 3)14CHILDREN = (15    (HALF, (Fraction(0), Fraction(0))),16    (THIRD, (Fraction(2, 3), Fraction(0))),17    (THIRD, (Fraction(2, 3), Fraction(1, 3))),18    (THIRD, (Fraction(2, 3), Fraction(2, 3))),19    (THIRD, (Fraction(0), Fraction(2, 3))),20    (THIRD, (Fraction(1, 3), Fraction(2, 3))),21)22LN2, LN3 = math.log(2), math.log(3)23NUMER, DENOM = 65, 4124HEIGHT = 60.025TERMS = 1926UMAX = 300.027WINDOW = 10.028WINDOW_STARTS = (10.0, 20.0, 40.0, 80.0, 160.0, 290.0)29SPEC_LO = 50.030SPEC_STEP = 0.00231NBINS = 4032OMEGA_MAX = 600.033PEAKS = 1034TEXTBOOK = 0.7675115443 + 45.55415979j35TOLERANCE = 1e-436RIPPLE_BINS = 2437RIPPLE_PERIODS = 438RIPPLE_CENTRES = (((0.0, 0.0), 0.5, 2 / 3), ((1.0, 0.0), 1 / 3, 1 / 3), ((1.0, 1.0), 1 / 3, 1 / 3))39STEP = 2e-440T0 = time.time()4142def clock(label):43    print(f"[{time.time() - T0:6.1f}s] {label}")4445def gp(script):46    return subprocess.run(["gp", "-q"], input="default(realprecision, 60);\n" + script, capture_output=True, text=True, check=True).stdout4748def number(text):49    return float(text.strip().replace(" E", "E"))5051def dimension(k):52    return number(gp(f"print(solve(s = 0, 4, 2^-s + {k}*3^-s - 1))"))5354def left_edge(k):55    return number(gp(f"print(solve(s = -2, 4, {k}*3^-s - 1 - 2^-s))"))5657def overlap(a, b):58    (ra, (xa, ya)), (rb, (xb, yb)) = a, b59    return xa < xb + rb and xb < xa + ra and ya < yb + rb and yb < ya + ra6061def letter_grid():62    grid = [["." for _ in range(6)] for _ in range(6)]63    for index, (r, (x, y)) in enumerate(CHILDREN):64        n = int(r * 6)65        for i in range(n):66            for j in range(n):67                grid[int(y * 6) + j][int(x * 6) + i] = "H" if index == 0 else str(index)68    return ["".join(row) for row in reversed(grid)]6970def words(level):71    cells = [(Fraction(1), Fraction(0), Fraction(0), 0, 0)]72    for _ in range(level):73        cells = [(size * r, x + size * tx, y + size * ty, a + (r == HALF), b + (r == THIRD)) for size, x, y, a, b in cells for r, (tx, ty) in CHILDREN]74    return cells7576def render(level):77    side = 6**level78    image = np.zeros((side, side), dtype=np.int32)79    for size, x, y, _, _ in words(level):80        n, i, j = int(size * side), int(x * side), int(y * side)81        image[j:j + n, i:i + n] += 182    return image8384def census(depth):85    counts = {}86    frontier = [(Fraction(1), Fraction(0), Fraction(0), 0, 0)]87    for _ in range(depth + 1):88        nxt = []89        for size, x, y, a, b in frontier:90            counts[(a, b)] = counts.get((a, b), 0) + 191            if a + b < depth:92                nxt.extend((size * r, x + size * tx, y + size * ty, a + (r == HALF), b + (r == THIRD)) for r, (tx, ty) in CHILDREN)93        frontier = nxt94    return counts9596def patch():97    print("the letter on the 6-grid, H the 1/2 child, 1..5 the 1/3 children, . the gap")98    for row in letter_grid():99        print("  " + row)100    inside = all(0 <= x and x + r <= 1 and 0 <= y and y + r <= 1 for r, (x, y) in CHILDREN)101    disjoint = all(not overlap(CHILDREN[i], CHILDREN[j]) for i in range(6) for j in range(i + 1, 6))102    print(f"children inside the closed unit square: {inside}   open images pairwise disjoint (open set condition): {disjoint}")103    area = sum(r * r for r, _ in CHILDREN)104    print(f"area of the children {area} = 1/4 + {K}/9; gap {1 - area}")105    out = gp(f"D = solve(s = 0, 4, 2^-s + {K}*3^-s - 1); print(D); print(log(29)/log(6)); print(log(2)/log(3))").split()106    print(f"dimension D solving 2^-s + {K}*3^-s = 1: {out[0][:22]}")107    print(f"level render reading log(29)/log(6) = {out[1][:12]}   log 2/log 3 = {out[2][:12]}")108    for level in range(1, 5):109        image = render(level)110        fill = int((image > 0).sum())111        print(f"  level {level}: side 6^{level} = {6**level}, fill {fill} = (9 + 4*{K})^{level}: {fill == (9 + 4 * K) ** level}   max cover {image.max()}")112    tile = (render(1) > 0).astype(np.int32)113    kron = int((np.kron(tile, tile) != (render(2) > 0)).sum())114    print(f"level 2 render against the Kronecker square of the level 1 tile: {kron} of {36 * 36} cells differ")115    depth = 7116    counts = census(depth)117    good = all(counts[(a, b)] == comb(a + b, a) * K**b for (a, b) in counts)118    print(f"cell census to word length {depth}: {len(counts)} sizes 2^-a 3^-b, every count C(a+b, a) {K}^b: {good}")119    clock("patch")120121def moran(s, k):122    return 1 - np.exp(-s * LN2) - k * np.exp(-s * LN3)123124def moran_prime(s, k):125    return LN2 * np.exp(-s * LN2) + k * LN3 * np.exp(-s * LN3)126127def seeds(k, height):128    out = gp(f"r = polroots(1 - x^{DENOM} - {k}*x^{NUMER}); for(j = 1, #r, z = -{DENOM}*log(r[j])/log(2); print(real(z), \" \", imag(z)))")129    rows = np.array([[float(t) for t in line.split()] for line in out.strip().splitlines()])130    z = rows[:, 0] + 1j * rows[:, 1]131    period = 2 * math.pi * DENOM / LN2132    z = np.concatenate([z + 1j * period * m for m in (-1, 0, 1)])133    return z[np.abs(z.imag) <= height + 2]134135def newton(z, k):136    with np.errstate(all="ignore"):137        for _ in range(60):138            z = z - moran(z, k) / moran_prime(z, k)139    return z140141def roots_in_box(k, lo, hi, height):142    re = np.linspace(lo, hi, 9)143    im = np.linspace(-height, height, int(8 * height) + 1)144    grid = (re[:, None] + 1j * im[None, :]).ravel()145    lattice = seeds(k, height)146    z = newton(np.concatenate([lattice, grid]), k)147    with np.errstate(all="ignore"):148        good = np.isfinite(z) & (np.abs(moran(z, k)) < 1e-12) & (z.real >= lo) & (z.real <= hi) & (np.abs(z.imag) <= height)149    found = []150    for w in sorted(z[good], key=lambda w: (w.imag, w.real)):151        if all(abs(w - q) > 1e-8 for q in found):152            found.append(w)153    polished = newton(lattice, k)154    with np.errstate(all="ignore"):155        hit = np.isfinite(polished) & (np.abs(moran(polished, k)) < 1e-12) & (np.abs(polished.imag) <= height) & (polished.real >= lo) & (polished.real <= hi)156    from_lattice = sum(1 for w in found if np.min(np.abs(polished[hit] - w), initial=1.0) < 1e-8)157    return np.array(found), from_lattice158159def edge_floor(k, lo, hi, height):160    x = np.linspace(lo, hi, 4001)161    return float(np.abs(moran(x + 1j * height, k)).min())162163def winding(k, lo, hi, height, step):164    corners = [lo - 1j * height, hi - 1j * height, hi + 1j * height, lo + 1j * height]165    lipschitz = LN2 * 2 ** (-lo) + k * LN3 * 3 ** (-lo)166    turns, margin = 0.0, math.inf167    for a, b in zip(corners, corners[1:] + corners[:1]):168        n = int(math.ceil(abs(b - a) / step))169        path = a + (b - a) * np.arange(n + 1) / n170        values = moran(path, k)171        margin = min(margin, float((np.abs(values[:-1]) / (lipschitz * abs(b - a) / n)).min()))172        turns += float(np.angle(values[1:] / values[:-1]).sum())173    return int(round(turns / (2 * math.pi))), margin174175def box(k):176    d, edge = dimension(k), left_edge(k)177    lo, hi = edge - 0.25, d + 0.25178    height = HEIGHT179    while edge_floor(k, lo, hi, height) < 0.02:180        height += 0.5181    found, from_lattice = roots_in_box(k, lo, hi, height)182    count, margin = winding(k, lo, hi, height, STEP)183    return d, edge, lo, hi, height, found, from_lattice, count, margin184185def offsets(found):186    upper = found[found.imag > 1e-9]187    out = []188    for base in (2, 3, 6):189        omega = 2 * math.pi / math.log(base)190        ratio = upper.imag / omega191        out.append(math.floor(1000 * float(np.abs(ratio - np.round(ratio)).max())) / 1000)192    return out193194def convergents():195    script = "\n".join([196        f"k = {K}; D = solve(s = 0, 4, 2^-s + k*3^-s - 1); P = 2^-D; Q = k*3^-D; mu = P*log(2) + Q*log(3);",197        "mf(z) = 1 - 2^-z - k*3^-z;",198        "mg(z) = log(2)*2^-z + k*log(3)*3^-z;",199        f"c = contfrac(log(3)/log(2), , {TERMS}); p0 = 1; q0 = 0; p1 = c[1]; q1 = 1;",200        "for(j = 2, #c - 1, p2 = c[j]*p1 + p0; q2 = c[j]*q1 + q0; p0 = p1; q0 = q1; p1 = p2; q1 = q2; w = D + I*2*Pi*q1/log(2); for(i = 1, 80, w = w - mf(w)/mg(w)); th = 2*Pi*(q1*log(3)/log(2) - p1); pr = P*Q*log(2)^2*th^2/(2*mu^3); if(q1 > 1, print(q1, \" \", p1, \" \", imag(w), \" \", D - real(w), \" \", pr, \" \", abs(mf(w)))))",201    ])202    rows = [line.replace(" E", "E").split() for line in gp(script).strip().splitlines()]203    print(f"roots near D + 2 pi i q/ln 2 for the convergents p/q of log2(3), k = {K}: D - Re w against P Q (ln 2)^2 theta^2/(2 f'(D)^3), theta = 2 pi (q log2 3 - p)")204    for q, p, im, gap, pred, res in rows:205        print(f"  q {int(q):>7}  p {int(p):>7}  Im w {float(im):>14.6f}  D - Re w {float(gap):.6e}  law {float(pred):.6e}  ratio {float(gap) / float(pred):.6f}  |f(w)| {float(res):.0e}")206207def poles():208    print(f"complex dimensions: the zeros of 1 - 2^-s - k 3^-s in the box Re [D_l - 0.25, D + 0.25], Im [-T, T], T >= {HEIGHT}")209    print(f"seeds from PARI polroots of 1 - x^{DENOM} - k x^{NUMER} ({NUMER}/{DENOM} a convergent of log2 3), then Newton; the count certified by the winding number")210    print("   k        D        D_l      T   roots  winding  margin  from-lattice  Re min      Re max     offset 2pi/ln2 ln3 ln6, floored")211    worst = 0.0212    for k in range(1, 28):213        d, edge, lo, hi, height, found, from_lattice, count, margin = box(k)214        o = offsets(found)215        print(f"  {k:>2}  {d:.6f}  {edge:+.6f}  {height:5.1f}  {found.size:>4}  {count:>5}  {margin:7.1f}  {from_lattice:>5}       {found.real.min():+.6f}  {found.real.max():+.6f}   {o[0]:.3f} {o[1]:.3f} {o[2]:.3f}")216        worst = max(worst, found.real.min() - edge)217        if k == K:218            keep = found219        if k == 1:220            textbook = found[np.argmin(np.abs(found - TEXTBOOK))]221    print(f"largest gap from D_l to the lowest real part over k = 1..27, rounded up: {math.ceil(worst * 10000) / 10000:.4f}")222    z = complex(TEXTBOOK)223    for _ in range(60):224        z -= (1 - 2**-z - 2 ** (-485 * z / 306)) / (LN2 * 2**-z + 485 / 306 * LN2 * 2 ** (-485 * z / 306))225    print(f"k = 1, the 2-3 nonlattice equation: the printed root .7675115443 + 45.55415979 i; the root of its lattice approximant 1 - 2^-s - 2^(-485 s/306) {z.real:.10f} + {z.imag:.8f} i, at {abs(z - TEXTBOOK):.1e}; the true root {textbook.real:.10f} + {textbook.imag:.8f} i, at {abs(textbook - TEXTBOOK):.1e}")226    print(f"the first complex dimensions of the patch, k = {K}, Im >= 0:")227    for w in keep[keep.imag > -1e-9]:228        print(f"  {w.real:+.6f} {w.imag:+.6f} i")229    clock("poles")230    convergents()231    clock("convergents")232233def breakpoints(k, top):234    rows = []235    for b in range(int(top / LN3) + 1):236        for a in range(int((top - b * LN3) / LN2) + 1):237            rows.append((a * LN2 + b * LN3, a, b))238    rows.sort()239    exact = [2**a * 3**b for _, a, b in rows]240    ordered = all(x < y for x, y in zip(exact, exact[1:]))241    weights = [comb(a + b, a) * k**b for _, a, b in rows]242    totals, running = [], 0243    for w in weights:244        running += w245        totals.append(running)246    return np.array([u for u, _, _ in rows]), exact, totals, ordered247248def cells_at_least(exact, totals, x):249    i = bisect_right(exact, x)250    return totals[i - 1] if i else 0251252def renewal_check(k, exact, totals):253    return all(totals[j] == 1 + cells_at_least(exact, totals, Fraction(v, 2)) + k * cells_at_least(exact, totals, Fraction(v, 3)) for j, v in enumerate(exact))254255def normalized(u_breaks, log_totals, d, grid):256    i = np.searchsorted(u_breaks, grid, side="right") - 1257    return np.exp(log_totals[i] - d * grid)258259def fold(u, g, period, nbins):260    phase = np.floor((u % period) / period * nbins).astype(int)261    means = np.bincount(phase, weights=g, minlength=nbins) / np.maximum(np.bincount(phase, minlength=nbins), 1)262    centred = g - g.mean()263    return 1 - float(((g - means[phase]) ** 2).sum() / (centred**2).sum())264265def spectrum(u, g, top):266    y = (g - g.mean()) * np.blackman(g.size)267    pad = 1 << 21268    power = np.abs(np.fft.rfft(y, pad)) ** 2269    omega = 2 * np.pi * np.fft.rfftfreq(pad, u[1] - u[0])270    keep = (omega > 2.0) & (omega < OMEGA_MAX)271    w, pw = omega[keep], power[keep]272    peaks = np.where((pw[1:-1] > pw[:-2]) & (pw[1:-1] >= pw[2:]))[0] + 1273    peaks = peaks[np.argsort(pw[peaks])[::-1][:top]]274    amplitude = np.sqrt(pw[peaks]) / (np.blackman(g.size).sum() / 2)275    return w[peaks], amplitude276277def roots_to(k, height):278    found, _ = roots_in_box(k, left_edge(k) - 0.25, dimension(k) + 0.25, height)279    return found[found.imag > 0]280281def count():282    d = dimension(K)283    p, q = 2**-d, K * 3**-d284    mu = p * LN2 + q * LN3285    limit = 1 / (d * mu)286    u_breaks, exact, totals, ordered = breakpoints(K, UMAX)287    print(f"cells of size 2^-a 3^-b >= r = e^-U, U <= {UMAX:.0f}: {len(exact)} sizes, sorted exactly by 2^a 3^b: {ordered}")288    print(f"N(e^-U) = sum of C(a+b, a) {K}^b; N({UMAX:.0f}) has {len(str(totals[-1]))} digits")289    clock("breakpoints")290    print(f"renewal identity N(U) = 1 + N(U - ln 2) + {K} N(U - ln 3) exact at every breakpoint: {renewal_check(K, exact, totals)}")291    clock("renewal")292    log_totals = np.array([math.log(t) for t in totals])293    print(f"limit N(r) r^D -> C = 1/(D f'(D)), f'(D) = 2^-D ln 2 + {K} 3^-D ln 3 = {mu:.9f}: {limit:.9f}")294    print(f"   window U     min F/C     max F/C     mean F/C - 1    swing      swing sqrt(U)   carpet swing")295    dc = math.log(8) / LN3296    for start in WINDOW_STARTS:297        grid = np.linspace(start, start + WINDOW, 200001)298        f = normalized(u_breaks, log_totals, d, grid) / limit299        m = np.floor(grid / LN3)300        carpet = (8 ** (m + 1) - 1) / 7 * np.exp(-dc * grid)301        cswing = (carpet.max() - carpet.min()) / carpet.mean()302        swing = f.max() - f.min()303        print(f"  [{start:5.0f},{start + WINDOW:5.0f}]  {f.min():.6f}   {f.max():.6f}   {f.mean() - 1:+.6f}      {swing:.6f}   {swing * math.sqrt(start):.6f}       {cswing:.6f}")304    clock("windows")305    u = np.arange(SPEC_LO, UMAX, SPEC_STEP)306    g = np.log(normalized(u_breaks, log_totals, d, u))307    m = np.floor(u / LN3)308    gc = np.log((8 ** (m + 1) - 1) / 7) - dc * u309    print(f"folded variance of g(u) = ln N(e^-u) - D u on [{SPEC_LO:.0f}, {UMAX:.0f}], {NBINS} bins:   ln 2    ln 3    ln 6")310    print(f"  patch                                               {fold(u, g, LN2, NBINS):.3f}   {fold(u, g, LN3, NBINS):.3f}   {fold(u, g, math.log(6), NBINS):.3f}")311    print(f"  carpet control                                      {fold(u, gc, LN2, NBINS):.3f}   {fold(u, gc, LN3, NBINS):.3f}   {fold(u, gc, math.log(6), NBINS):.3f}")312    roots = roots_to(K, OMEGA_MAX + 5)313    print(f"Blackman periodogram of g on [{SPEC_LO:.0f}, {UMAX:.0f}], the {PEAKS} highest peaks in ({2.0}, {OMEGA_MAX:.0f}), each against the nearest complex dimension:")314    w, a = spectrum(u, g, PEAKS)315    taper = np.blackman(u.size)316    for wi, ai in sorted(zip(w, a)):317        near = roots[np.argmin(np.abs(roots.imag - wi))]318        decay = float((taper * np.exp(-(d - near.real) * u)).sum() / taper.sum())319        predicted = 2 * decay / (abs(near * moran_prime(near, K)) * limit)320        print(f"  peak {wi:8.3f}  amplitude {ai:.3e}  root {near.real:.6f} + {near.imag:.3f} i  D - Re {d - near.real:.2e}  offset {wi - near.imag:+.3f}  residue amplitude {predicted:.3e}")321    wc, ac = spectrum(u, gc, 6)322    print("carpet control peaks, as multiples of 2 pi/ln 3:", " ".join(f"{x / (2 * math.pi / LN3):.3f}" for x in sorted(wc)))323    clock("count")324325def ball_mass(cx, cy, radius, d):326    ratios = np.array([float(r) for r, _ in CHILDREN])327    shifts = np.array([[float(x), float(y)] for _, (x, y) in CHILDREN])328    weights = ratios**d329    x, y, size, mass = np.zeros(1), np.zeros(1), np.ones(1), np.ones(1)330    inside, r2 = 0.0, radius * radius331    while size.size:332        far = np.maximum(np.abs(x - cx), np.abs(x + size - cx)) ** 2 + np.maximum(np.abs(y - cy), np.abs(y + size - cy)) ** 2333        nx = np.clip(cx, x, x + size) - cx334        ny = np.clip(cy, y, y + size) - cy335        near = nx * nx + ny * ny336        full = far <= r2 * (1 - 1e-12)337        inside += float(mass[full].sum())338        open_ = ~full & (near < r2 * (1 + 1e-12))339        x, y, size, mass = x[open_], y[open_], size[open_], mass[open_]340        if size.size == 0 or size.max() < TOLERANCE * radius:341            return inside, inside + float(mass.sum())342        x = (x[:, None] + size[:, None] * shifts[None, :, 0]).ravel()343        y = (y[:, None] + size[:, None] * shifts[None, :, 1]).ravel()344        mass = (mass[:, None] * weights[None, :]).ravel()345        size = (size[:, None] * ratios[None, :]).ravel()346    return inside, inside347348def ripple():349    d = dimension(K)350    print(f"the spin mass M(r) = mu(B(c, r)) of the natural measure, weights 2^-D and 3^-D, D = {d:.9f}, enclosed to relative {TOLERANCE:g} by cells")351    print("   centre    map ratio  window   identity M(r) = ratio^D M(r/ratio)   ripple swing   bar     swing/bar  fold ln 2   fold ln 3   fold ln 6")352    for (cx, cy), ratio, window in RIPPLE_CENTRES:353        period = math.log(1 / ratio)354        steps = np.arange(RIPPLE_PERIODS * RIPPLE_BINS + RIPPLE_BINS) * period / RIPPLE_BINS355        radii = window * 0.999 * np.exp(-steps)356        bounds = np.array([ball_mass(cx, cy, r, d) for r in radii])357        mid = bounds.mean(axis=1)358        bar = float(np.max(np.log(bounds[:, 1] / bounds[:, 0])))359        lo, hi = bounds[RIPPLE_BINS:, 0], bounds[RIPPLE_BINS:, 1]360        scaled_lo, scaled_hi = ratio**d * bounds[:-RIPPLE_BINS, 0], ratio**d * bounds[:-RIPPLE_BINS, 1]361        meets = bool(np.all((lo <= scaled_hi * (1 + 1e-9)) & (scaled_lo <= hi * (1 + 1e-9))))362        g = np.log(mid) - d * np.log(radii)363        logr = np.log(radii)364        folds = [fold(logr, g, math.log(b), RIPPLE_BINS) for b in (2, 3, 6)]365        print(f"   ({cx},{cy})   1/{round(1 / ratio)}       {window:.4f}   bands meet at all {lo.size} shifted radii: {meets}          {g.max() - g.min():.5f}   {bar:.1e}   {math.floor((g.max() - g.min()) / bar):>6}     {folds[0]:.3f}       {folds[1]:.3f}       {folds[2]:.3f}")366    clock("ripple")367368VERBS = {"patch": patch, "poles": poles, "count": count, "ripple": ripple}369370if __name__ == "__main__":371    VERBS[sys.argv[1] if len(sys.argv) > 1 else "patch"]()