leaning_stack.py
18.8 kB · python · 477 lines
1import cmath2import math3import os4from fractions import Fraction56import numpy as np7from PIL import Image89N_LINEAR = 3010N_QUAD = 6011B_MAX = 1212DRIFTS = [Fraction(1, 2), Fraction(1, 3), Fraction(1, 4), Fraction(2, 5), Fraction(1, 6), Fraction(3, 8), Fraction(5, 9)]13LINEAR_DRIFT = Fraction(2, 7)14GAUSS_PRIMES = [5, 7, 11, 13]15ORIGIN_D = 6016WEYL_NS = [1000, 10000]17PAIR_MAX = 4018SHARE_MAX = 1219SHARE_D = 1620WIDE_B = 2021WIDE_D = 1622WIDE_PAIR = 12023DISTINCT_N = 2524RENDER_N = 3025PANEL = 50026DIGITS = 502728# EXACT IRRATIONAL DRIFTS2930def sqrt_fraction(k):31 scale = 10 ** DIGITS32 return Fraction(math.isqrt(k * scale * scale), scale)3334def silver():35 return sqrt_fraction(2) - 13637def golden():38 return (sqrt_fraction(5) - 1) / 23940# THE LEAN4142def phase_linear(n, drift):43 return n * drift4445def phase_quadratic(n, drift):46 return n * n * drift4748def lit(n, x, theta):49 return (n * x - theta).denominator == 15051def brightness_literal(top, x, drift, phase):52 return sum(1 for n in range(1, top + 1) if lit(n, x, phase(n, drift)))5354def reduced_points(top_b):55 out = [(0, 1)]56 for b in range(2, top_b + 1):57 for a in range(1, b):58 if math.gcd(a, b) == 1:59 out.append((a, b))60 return out6162# THE QUADRATIC CONGRUENCE6364def valuation(m, p):65 e = 066 while m % p == 0:67 m //= p68 e += 169 return e7071def prime_factors(m):72 out = []73 q = m74 p = 275 while p * p <= q:76 if q % p == 0:77 out.append(p)78 while q % p == 0:79 q //= p80 p += 181 if q > 1:82 out.append(q)83 return out8485def local_classes(p, a, b, c, d):86 beta = valuation(b, p)87 delta = valuation(d, p)88 modulus = p ** max(beta, (delta + 1) // 2)89 if beta >= 1 and beta <= delta < 2 * beta:90 k = 2 * beta - delta91 unit = (a * (d // p ** delta) * pow(c * (b // p ** beta), -1, p ** k)) % p ** k92 return modulus, [0, (p ** (delta - beta) * unit) % modulus]93 return modulus, [0]9495def crt_merge(state, modulus, classes):96 period, residues = state97 step = period * modulus98 inv = pow(period % modulus, -1, modulus)99 out = []100 for r in residues:101 for s in classes:102 out.append((r + period * ((s - r) * inv % modulus)) % step)103 return step, sorted(out)104105def solution_classes(a, b, c, d):106 state = (1, [0])107 for p in prime_factors(b * d):108 state = crt_merge(state, *local_classes(p, a, b, c, d))109 return state110111def count_class(top, period, residue):112 if residue == 0:113 return top // period114 return 0 if residue > top else (top - residue) // period + 1115116def brightness_form(top, a, b, c, d):117 period, residues = solution_classes(a, b, c, d)118 return sum(count_class(top, period, r) for r in residues)119120def first_lit(a, b, c, d):121 period, residues = solution_classes(a, b, c, d)122 return min(period if r == 0 else r for r in residues)123124def square_root_kernel(d):125 out = 1126 for p in prime_factors(d):127 out *= p ** ((valuation(d, p) + 1) // 2)128 return out129130def middle_primes(b, d):131 return [p for p in prime_factors(b) if valuation(b, p) <= valuation(d, p) < 2 * valuation(b, p)]132133def period_form(b, d):134 k = square_root_kernel(d)135 return b * k // math.gcd(b, k)136137def halving_rule(b, d):138 period = period_form(b, d)139 beta = valuation(b, 2)140 return period // 2 if beta >= 1 and valuation(d, 2) == 2 * beta - 1 else period141142# THE SOLUTION SET READ LITERALLY143144def literal_mask(a, b, c, d):145 m = b * d146 n = np.arange(m, dtype=np.int64)147 return (n * (a * d - b * c * n)) % m == 0148149def divisors(m):150 out = []151 k = 1152 while k * k <= m:153 if m % k == 0:154 out.append(k)155 out.append(m // k)156 k += 1157 return sorted(set(out))158159def minimal_period(mask):160 m = mask.size161 count = int(mask.sum())162 for q in divisors(m):163 if count % (m // q):164 continue165 if np.array_equal(mask, np.roll(mask, q)):166 return q167 return m168169def formula_set(a, b, c, d):170 period, residues = solution_classes(a, b, c, d)171 return {r + period * k for r in residues for k in range(b * d // period)}172173def shares_point(m, n, drift):174 return (m * n // math.gcd(m, n) * (m - n) * drift).denominator == 1175176def layer_points(n, drift):177 theta = (n * n * drift) % 1178 return {(k + theta) / n for k in range(n)}179180# GAUSS SUMS181182def legendre(a, p):183 a %= p184 if a == 0:185 return 0186 return 1 if pow(a, (p - 1) // 2, p) == 1 else -1187188def phase_census(c, p):189 out = {}190 for n in range(p):191 j = c * n * n % p192 out[j] = out.get(j, 0) + 1193 return out194195def gauss_sum(coefficient, m):196 return sum(cmath.exp(2j * math.pi * (coefficient * n * n % m) / m) for n in range(m))197198def epsilon(p):199 return 1 if p % 4 == 1 else 1j200201def quadratic_sum(a_coeff, b_coeff, m):202 return sum(cmath.exp(2j * math.pi * ((a_coeff * n * n + b_coeff * n) % m) / m) for n in range(m))203204# SECTIONS205206def linear_lean():207 print("LINEAR LEAN t_n = delta, phase n delta, N =", N_LINEAR)208 drift = LINEAR_DRIFT209 bad = 0210 for a, b in reduced_points(N_LINEAR):211 x = drift + Fraction(a, b)212 if brightness_literal(N_LINEAR, x, drift, phase_linear) != N_LINEAR // b:213 bad += 1214 off = [drift + Fraction(1, q) for q in (31, 37, 41, 43)]215 blank = sum(1 for x in off if brightness_literal(N_LINEAR, x, drift, phase_linear) != 0)216 print(" delta =", drift, " nodes tested", len(reduced_points(N_LINEAR)), " brightness != floor(N/b):", bad, " nonzero off the translate:", blank)217 for name, value in (("sqrt 2 - 1", math.sqrt(2) - 1), ("phi - 1", (math.sqrt(5) - 1) / 2)):218 bad = 0219 for a, b in reduced_points(N_LINEAR):220 x = value + a / b221 hit = sum(1 for n in range(1, N_LINEAR + 1) if abs(n * x - n * value - round(n * x - n * value)) < 1e-12)222 if hit != N_LINEAR // b:223 bad += 1224 closest = min(min(abs(n * (a / b) - n * value - round(n * (a / b) - n * value)) for n in range(1, N_LINEAR + 1)) for a, b in reduced_points(N_LINEAR)[1:])225 print(" delta = %s brightness != floor(N/b) at 1e-12: %d closest approach at a rational x: %.6f" % (name, bad, closest))226227def quadratic_lean():228 print("QUADRATIC LEAN t_n = n delta, phase n^2 delta, N =", N_QUAD, " all reduced a/b with b <=", B_MAX)229 points = reduced_points(B_MAX)230 for drift in DRIFTS:231 c, d = drift.numerator, drift.denominator232 bad = sum(1 for a, b in points if brightness_literal(N_QUAD, Fraction(a, b), drift, phase_quadratic) != brightness_form(N_QUAD, a, b, c, d))233 set_bad = 0234 period_bad = 0235 halved = 0236 for a, b in points:237 mask = literal_mask(a, b, c, d)238 if formula_set(a, b, c, d) != set(int(v) for v in np.nonzero(mask)[0]):239 set_bad += 1240 least = minimal_period(mask)241 if least != halving_rule(b, d):242 period_bad += 1243 if least < period_form(b, d):244 halved += 1245 top = max(brightness_form(N_QUAD, a, b, c, d) for a, b in points)246 best = [Fraction(a, b) for a, b in points if brightness_form(N_QUAD, a, b, c, d) == top]247 origin = brightness_form(N_QUAD, 0, 1, c, d)248 print(" delta = %s points %d brightness mismatches %d solution-set mismatches %d minimal-period mismatches %d minimal period below lcm(b, d*) at %d points" % (drift, len(points), bad, set_bad, period_bad, halved))249 print(" top brightness %d at %s origin %d origin is top: %s" % (top, ", ".join(str(v) for v in best), origin, "yes" if origin == top else "no"))250251def minimal_period_law():252 print("THE MINIMAL PERIOD lcm(b, d*) is a period; the least one halves it exactly when v_2(b) >= 1 and v_2(d) = 2 v_2(b) - 1")253 points = reduced_points(WIDE_B)254 drifts = [Fraction(c, e) for e in range(1, WIDE_B + 1) for c in range(1, e + 1) if math.gcd(c, e) == 1]255 tuples = 0256 halved = 0257 bad = 0258 for drift in drifts:259 c, d = drift.numerator, drift.denominator260 for a, b in points:261 tuples += 1262 least = minimal_period(literal_mask(a, b, c, d))263 if least < period_form(b, d):264 halved += 1265 if least != halving_rule(b, d):266 bad += 1267 print(" tuples %d over b, d <= %d minimal period below lcm(b, d*) at %d breaches of the rule %d" % (tuples, WIDE_B, halved, bad))268 for a, b, c, d in ((1, 2, 1, 2), (1, 6, 1, 2), (1, 4, 3, 8)):269 mask = literal_mask(a, b, c, d)270 print(" x = %d/%d delta = %d/%d lcm(b, d*) = %d classes %s minimal period %d" % (a, b, c, d, period_form(b, d), sorted(solution_classes(a, b, c, d)[1]), minimal_period(mask)))271272def lit_set():273 print("THE LIT SET first layer lighting a/b, and its dependence on the numerator")274 points = reduced_points(B_MAX)275 for drift in DRIFTS:276 c, d = drift.numerator, drift.denominator277 bad = 0278 for a, b in points:279 reach = first_lit(a, b, c, d)280 literal = next((n for n in range(1, N_QUAD + 1) if lit(n, Fraction(a, b), phase_quadratic(n, drift))), None)281 if literal != (reach if reach <= N_QUAD else None):282 bad += 1283 generic = sum(1 for a, b in points if first_lit(a, b, c, d) == period_form(b, d))284 print(" delta = %s first-lit mismatches %d points first lit at lcm(b, d*) itself %d of %d" % (drift, bad, generic, len(points)))285 print(" numerator dependence, delta = 1/4:")286 for a in (1, 3):287 period, residues = solution_classes(a, 4, 1, 4)288 print(" x = %d/4 classes mod %d: %s first lit n = %d B_61 = %d" % (a, period, residues, first_lit(a, 4, 1, 4), brightness_form(61, a, 4, 1, 4)))289 split = [(a, b, Fraction(c, d)) for b in range(2, B_MAX + 1) for d in (2, 4, 8, 9) for c in range(1, d) if math.gcd(c, d) == 1 for a in range(1, b) if math.gcd(a, b) == 1 and first_lit(a, b, c, d) != first_lit(1, b, c, d)]290 print(" reduced (a/b, delta) pairs with b <= %d and d in 2, 4, 8, 9 whose first lit layer moves with the numerator: %d" % (B_MAX, len(split)))291292def zeta_em(s, terms=400000):293 k = np.arange(1, terms + 1, dtype=np.float64)294 total = float(np.sum(k ** (-float(s))))295 m = float(terms)296 return total + m ** (1 - s) / (s - 1) - 0.5 * m ** (-s) + s * m ** (-s - 1) / 12.0297298def inner_roots(limit):299 inner = np.ones(limit + 1, dtype=np.int64)300 k = 2301 while k * k <= limit:302 inner[k * k :: k * k] = k303 k += 1304 return inner305306def origin_law():307 print("THE ORIGIN x = 0, brightness floor(N/d*) with d* the least k with d | k^2")308 bad = 0309 for d in range(1, ORIGIN_D + 1):310 for c in range(1, d + 1):311 if math.gcd(c, d) != 1:312 continue313 drift = Fraction(c, d)314 if brightness_literal(N_QUAD, Fraction(0), drift, phase_quadratic) != N_QUAD // square_root_kernel(d):315 bad += 1316 kernels = [square_root_kernel(d) for d in range(1, 31)]317 print(" drifts tested to d =", ORIGIN_D, " mismatches", bad)318 print(" d* for d = 1..30:", ", ".join(str(v) for v in kernels))319 limit = 10 ** 6320 d = np.arange(1, limit + 1, dtype=np.float64)321 star = d / inner_roots(limit)[1:].astype(np.float64)322 breaches = sum(1 for k in range(1, 1001) if int(star[k - 1]) != square_root_kernel(k))323 print(" d/A000188(d) against the local form to d = 1000: breaches", breaches)324 for s in (1, 2):325 partials = [float(np.sum(1.0 / (star[:cut] * d[:cut] ** s))) for cut in (10 ** 4, 10 ** 5, 10 ** 6)]326 target = zeta_em(2 * s + 1) * zeta_em(s + 1) / zeta_em(2 * s + 2)327 print(" s = %d partial sums %s zeta(%d) zeta(%d)/zeta(%d) = %.6f" % (s, ", ".join("%.6f" % v for v in partials), 2 * s + 1, s + 1, 2 * s + 2, target))328329def gauss_sums():330 print("GAUSS SUMS the phase census of the quadratic lean at prime drift denominator, and the centred twist S(c,p) = (c/p) S(1,p)")331 for p in GAUSS_PRIMES:332 residue = next(c for c in range(1, p) if legendre(c, p) == 1)333 nonresidue = next(c for c in range(1, p) if legendre(c, p) == -1)334 for c in (residue, nonresidue):335 census = phase_census(c, p)336 inverse = pow(c, -1, p)337 bad = sum(1 for j in range(p) if census.get(j, 0) != (1 if j == 0 else 1 + legendre(j * inverse, p)))338 value = gauss_sum(c, p)339 claim = legendre(c, p) * epsilon(p) * math.sqrt(p)340 twist = sum(1 for j in range(p) if census.get(j, 0) - 1 != legendre(c, p) * (phase_census(1, p).get(j, 0) - 1))341 print(" p = %2d c = %2d (c/p) = %2d distinct phases %d census breaches %d centred twist breaches %d S(c,p) = %.6f%+.6fi (c/p) eps_p sqrt p = %.6f%+.6fi" % (p, c, legendre(c, p), len(census), bad, twist, value.real, value.imag, claim.real, claim.imag))342 print(" the count of solutions as a Fourier sum, R = (1/m) sum_h sum_n e(h Q(n)/m)")343 for a, b, c, d in ((1, 1, 1, 5), (1, 2, 1, 3), (1, 4, 1, 4), (1, 6, 1, 6), (5, 12, 3, 8)):344 m = b * d345 total = sum(quadratic_sum(h * b * c, -h * a * d, m) for h in range(m)) / m346 period, residues = solution_classes(a, b, c, d)347 exact = len(residues) * (m // period)348 print(" a/b = %d/%d delta = %d/%d m = %2d Fourier count %.6f%+.6fi classes mod %d %s R = %d" % (a, b, c, d, m, total.real, total.imag, period, residues, exact))349350def irrational_lean():351 print("IRRATIONAL DRIFT no two layers share a point, brightness at most 1 everywhere")352 for name, drift in (("sqrt 2 - 1", silver()), ("phi - 1", golden())):353 worst = None354 for n in range(1, PAIR_MAX + 1):355 for m in range(1, n):356 span = m * n // math.gcd(m, n)357 t = drift * (m - n) * span358 frac = t - math.floor(t)359 gap = Fraction(1, span) * min(frac, 1 - frac)360 if worst is None or gap < worst[0]:361 worst = (gap, m, n)362 print(" delta = %s closest two layers to N = %d come: %.3e at (m, n) = (%d, %d)" % (name, PAIR_MAX, float(worst[0]), worst[1], worst[2]))363 row = []364 for top in WEYL_NS:365 q = drift.denominator366 p = drift.numerator367 values = sorted((n * n * p) % q for n in range(1, top + 1))368 star = 0.0369 for i, r in enumerate(values, start=1):370 u = r / q371 star = max(star, i / top - u, u - (i - 1) / top)372 row.append("N = %5d D*_N = %.6f 1/sqrt N = %.6f ratio %.4f" % (top, star, 1 / math.sqrt(top), star * math.sqrt(top)))373 for line in row:374 print(" ", line)375376# ADVERSARIAL PASS377378def sharing_law():379 print("SHARING layers m and n share a lit point iff lcm(m, n)(m - n) delta is an integer")380 drifts = [Fraction(c, d) for d in range(1, SHARE_D + 1) for c in range(1, d + 1) if math.gcd(c, d) == 1]381 pairs = 0382 sharing = 0383 bad = 0384 loose = 0385 slack = 0386 witness = None387 for drift in drifts:388 for n in range(2, SHARE_MAX + 1):389 for m in range(1, n):390 pairs += 1391 truth = len(layer_points(m, drift) & layer_points(n, drift)) > 0392 if truth:393 sharing += 1394 if truth != shares_point(m, n, drift):395 bad += 1396 if not truth and (Fraction(n) * (m - n) * drift).denominator == 1:397 loose += 1398 if not truth and (Fraction(m * n) * (m - n) * drift).denominator == 1:399 slack += 1400 if witness is None:401 witness = (m, n, drift)402 print(" pairs %d over m < n <= %d and every reduced c/d with d <= %d sharing %d criterion failures %d" % (pairs, SHARE_MAX, SHARE_D, sharing, bad))403 print(" the rational reading is vacuous at rational drift, every one of the %d non-sharing pairs having n (m - n) delta rational; n (m - n) delta in Z holds at %d of them and m n (m - n) delta in Z at %d, first at (m, n, delta) = (%d, %d, %s), so both are weaker than the criterion" % (pairs - sharing, loose, slack, witness[0], witness[1], witness[2]))404405def adversarial():406 print("ADVERSARIAL the closed form against literal stacking on the widest domain this study runs")407 points = reduced_points(WIDE_B)408 drifts = [Fraction(c, d) for d in range(1, WIDE_D + 1) for c in range(1, d + 1) if math.gcd(c, d) == 1]409 bad = 0410 tested = 0411 for drift in drifts:412 c, d = drift.numerator, drift.denominator413 for a, b in points:414 tested += 1415 if brightness_literal(N_QUAD, Fraction(a, b), drift, phase_quadratic) != brightness_form(N_QUAD, a, b, c, d):416 bad += 1417 print(" points %d drifts %d pairs %d brightness mismatches at N = %d: %d" % (len(points), len(drifts), tested, N_QUAD, bad))418 for name, drift in (("sqrt 2 - 1", silver()), ("phi - 1", golden())):419 seen = set()420 for n in range(1, DISTINCT_N + 1):421 theta = (n * n * drift) % 1422 for k in range(n):423 seen.add((k + theta) / n)424 worst = None425 for n in range(1, WIDE_PAIR + 1):426 for m in range(1, n):427 span = m * n // math.gcd(m, n)428 t = drift * (m - n) * span429 frac = t - math.floor(t)430 gap = Fraction(1, span) * min(frac, 1 - frac)431 if worst is None or gap < worst[0]:432 worst = (gap, m, n)433 print(" delta = %s lit points of layers 1..%d: %d distinct of %d drawn closest pair to N = %d: %.3e at (%d, %d)" % (name, DISTINCT_N, len(seen), DISTINCT_N * (DISTINCT_N + 1) // 2, WIDE_PAIR, float(worst[0]), worst[1], worst[2]))434435# RENDER436437def lean_axis(top, drift, resolution):438 counts = np.zeros(resolution, dtype=np.int64)439 for n in range(1, top + 1):440 theta = float((n * n * drift) % 1)441 for k in range(n):442 counts[int((k + theta) / n * resolution) % resolution] += 1443 return counts444445def panel_grey(drift):446 axis = lean_axis(RENDER_N, drift, PANEL)447 field = np.outer(axis, axis)448 levels = np.unique(field)449 rank = np.searchsorted(levels, field).astype(np.float64)450 grey = np.full(field.shape, 255.0)451 lit_mask = field > 0452 grey[lit_mask] = 190 - 190 * (rank[lit_mask] - 1) / max(levels.size - 2, 1)453 return grey.astype(np.uint8), int(field.max())454455def render():456 left, peak_left = panel_grey(Fraction(1, 5))457 right, peak_right = panel_grey(silver())458 sheet = np.full((PANEL, 2 * PANEL + 8), 255, dtype=np.uint8)459 sheet[:, :PANEL] = left460 sheet[:, PANEL + 8 :] = right461 path = "research/lab/py/leaning-stack/leaning-stack.png"462 Image.fromarray(sheet, mode="L").save(path, optimize=True)463 print("RENDER leaning-stack.png N = %d panels delta = 1/5 and delta = sqrt 2 - 1 peaks %d and %d %d bytes" % (RENDER_N, peak_left, peak_right, os.path.getsize(path)))464465def main():466 linear_lean()467 quadratic_lean()468 minimal_period_law()469 lit_set()470 origin_law()471 gauss_sums()472 sharing_law()473 adversarial()474 irrational_lean()475 render()476477main()