stack_dilations.py

1from fractions import Fraction2from math import cos, gcd, log, pi, sqrt34import numpy as np56C = 2.0 * sqrt(2.0) / pi78def sine_coefficient(n):9    return C / n if n % 2 == 1 else 0.01011def dilate_coefficient(n, k):12    total = 0.013    for j in range(n):14        total += ((-1) ** j) * (cos(k * pi * j / n) - cos(k * pi * (j + 1) / n)) / (k * pi)15    return sqrt(2.0) * total1617def check_dilation_shift(scales, harmonics):18    worst = 0.019    for n in range(1, scales + 1):20        for k in range(1, harmonics + 1):21            want = sine_coefficient(k // n) if k % n == 0 else 0.022            worst = max(worst, abs(dilate_coefficient(n, k) - want))23    print("shift: <s(nx), e_k> equals a_{k/n} when n divides k and 0 otherwise, max deviation "24          "%.2e over n <= %d and k <= %d, so the dilate s(nx) carries the Dirichlet series n^-s S(s)"25          % (worst, scales, harmonics))2627def grid_integral(m, n):28    lcm = m * n // gcd(m, n)29    total = 030    for k in range(lcm):31        a = (m * (2 * k + 1)) // (2 * lcm)32        b = (n * (2 * k + 1)) // (2 * lcm)33        total += 1 if (a + b) % 2 == 0 else -134    return Fraction(total, lcm)3536def symbol_entry(m, n):37    g = gcd(m, n)38    if (m // g) % 2 == 1 and (n // g) % 2 == 1:39        return Fraction(g * g, m * n)40    return Fraction(0)4142def symbol_series(m, n, terms):43    g = gcd(m, n)44    mm, nn = m // g, n // g45    total = 0.046    for j in range(1, terms + 1):47        total += sine_coefficient(j * nn) * sine_coefficient(j * mm)48    return total4950def check_gram_two_ways(cap, terms):51    pairs = 052    worst = 0.053    for m in range(1, cap + 1):54        for n in range(m, cap + 1):55            exact = grid_integral(m, n)56            assert exact == symbol_entry(m, n), "gram entry (%d,%d)" % (m, n)57            g = gcd(m, n)58            slack = C * C / ((m // g) * (n // g) * 2.0 * terms)59            gap = abs(symbol_series(m, n, terms) - float(exact))60            assert gap <= slack + 1e-12, "symbol series (%d,%d)" % (m, n)61            worst = max(worst, gap)62            pairs += 163    print("gram: the L^2(0,1) integral of s(mx)s(nx) on the lcm grid equals the symbol sum "64          "sum_j a_{jn'} a_{jm'} at all %d pairs m <= n <= %d, exact in rationals, reading %s at "65          "(3,5), %s at (3,9) and %s at (2,3)"66          % (pairs, cap, grid_integral(3, 5), grid_integral(3, 9), grid_integral(2, 3)))67    print("gram: the truncated symbol sum at %d terms sits within its own tail bound of the closed "68          "form at every pair, largest gap %.2e" % (terms, worst))6970def check_blocks(cap):71    off = 072    for m in range(1, cap + 1):73        for n in range(1, cap + 1):74            e = symbol_entry(m, n)75            am, an = (m & -m).bit_length() - 1, (n & -n).bit_length() - 176            if am != an:77                assert e == 0, "block leak (%d,%d)" % (m, n)78                off += 179            else:80                mo, no = m >> am, n >> an81                assert e == Fraction(gcd(mo, no) ** 2, mo * no), "block value (%d,%d)" % (m, n)82    print("blocks: the Gram entry vanishes whenever v_2(m) differs from v_2(n), %d ordered pairs to "83          "%d, and the block at v_2 = a is the odd Gram itself, so the full dilation system's Gram "84          "is a countable direct sum of copies of the odd one" % (off, cap))8586def exact_det(rows):87    a = [row[:] for row in rows]88    n = len(a)89    det = Fraction(1)90    for i in range(n):91        p = next(r for r in range(i, n) if a[r][i] != 0)92        if p != i:93            a[i], a[p] = a[p], a[i]94            det = -det95        det *= a[i][i]96        inv = Fraction(1) / a[i][i]97        for r in range(i + 1, n):98            f = a[r][i] * inv99            if f:100                for c in range(i, n):101                    a[r][c] -= f * a[i][c]102    return det103104def smith_product(odds):105    want = Fraction(1)106    for k in odds:107        num, den, t, p = 1, 1, k, 3108        if t % 2 == 0:109            t //= 2110        while p * p <= t:111            if t % p == 0:112                num *= p * p - 1113                den *= p * p114                while t % p == 0:115                    t //= p116            p += 2117        if t > 1:118            num *= t * t - 1119            den *= t * t120        want *= Fraction(num, den)121    return want122123def check_determinant(cap):124    for k in range(1, cap + 1):125        odds = list(range(1, 2 * k, 2))126        rows = [[Fraction(gcd(a, b) ** 2, a * b) for b in odds] for a in odds]127        assert exact_det(rows) == smith_product(odds), "det at K=%d" % k128    odds = list(range(1, 26, 2))129    print("det: the odd Gram determinant equals prod over odd k of prod over p | k of (1 - p^-2) at "130          "every K = 1..%d, the value over the %d odd scales n <= 25 being %s, the number "131          "lab/py/stack-levels prints" % (cap, len(odds), smith_product(odds)))132133def check_inverse(cap, primecap):134    ok = 0135    for n in range(1, cap + 1):136        total = Fraction(0)137        for d in range(1, n + 1):138            if n % d == 0 and d % 2 == 1 and (n // d) % 2 == 1:139                mu, t, p, sq = 1, d, 3, False140                while p * p <= t:141                    if t % p == 0:142                        t //= p143                        if t % p == 0:144                            sq = True145                            break146                        mu = -mu147                    p += 2148                if sq:149                    continue150                if t > 1:151                    mu = -mu152                total += Fraction(mu, d) * Fraction(1, n // d)153        assert total == (1 if n == 1 else 0), "inverse at n=%d" % n154        ok += 1155    sieve = bytearray([1]) * (primecap + 1)156    sieve[0] = sieve[1] = 0157    for p in range(2, int(primecap ** 0.5) + 1):158        if sieve[p]:159            sieve[p * p::p] = bytearray(len(sieve[p * p::p]))160    tail = sum(C / p for p in range(3, primecap + 1, 2) if sieve[p])161    print("inverse: the Dirichlet inverse of a_n/a_1 is mu(n)/n on the odd n and 0 on the even, "162          "checked at every n <= %d, so 1/S(s) = 1/(a_1 (1 - 2^-1-s) zeta(1 + s)) and the Mobius "163          "square wave, coefficients mu(n) a_n, has symbol a_1^2/S(s)" % ok)164    print("inverse: a_n/a_1 and its inverse are totally multiplicative, and the Riesz test sum over "165          "the primes reads %.4f already at p <= %d, a divergent sum by Mertens, so neither system "166          "is a Riesz basis; the verdict is invariant under scaling phi by a_1" % (tail, primecap))167168def spectrum(caps):169    print("spectrum: K, top odd scale N, lambda_max, lambda_min, condition number, "170          "lambda_max/(log N)^2, lambda_max/(log log N)^2")171    for k in caps:172        odds = np.arange(1, 2 * k, 2, dtype=np.float64)173        g = np.gcd(odds.astype(np.int64)[:, None], odds.astype(np.int64)[None, :]).astype(np.float64)174        mat = g * g / (odds[:, None] * odds[None, :])175        ev = np.linalg.eigvalsh(mat)176        top, bot = ev[-1], ev[0]177        n = 2 * k - 1178        print("spectrum: %5d %6d %10.5f %12.3e %12.3e %8.4f %8.4f"179              % (k, n, top, bot, top / bot, top / (log(n) ** 2), top / (log(log(n)) ** 2)))180181def main():182    check_dilation_shift(8, 60)183    check_gram_two_ways(12, 20000)184    check_blocks(64)185    check_determinant(13)186    check_inverse(400, 100000)187    spectrum([25, 50, 100, 150, 200])188189if __name__ == "__main__":190    main()