spun_stack.py
13.1 kB · python · 425 lines
1from fractions import Fraction2from math import gcd, isqrt, lcm3from sympy import Catalan as SCatalan4from sympy import N as N_5from sympy import Poly, Rational, Symbol, cos, cyclotomic_poly, minimal_polynomial, sin6from sympy import pi as spi7from sympy import sqrt as ssqrt8from sympy import sympify9from sympy import zeta as szeta1011pi = spi1213X = Symbol("x")14CYC = 36015PHI360 = Poly(cyclotomic_poly(CYC, X), X, domain="QQ")16SCHEDULE_N = 3017NEAR_BOX = 618STACK_N = 5019COUNT_NS = [50, 200, 800]20SPOT_DEGREES = [0, 7, 30, 45, 60, 90, 97, 113, 143, 180, 270, 271]21BASE_BOX = 1022ODD_N = 5523INCREMENTS = [24 (Fraction(0), "0"),25 (Fraction(18), "18"),26 (Fraction(45, 2), "22.5"),27 (Fraction(30), "30"),28 (Fraction(45), "45"),29 (Fraction(60), "60"),30 (Fraction(135, 2), "67.5"),31 (Fraction(36), "36"),32 (Fraction(15), "15"),33 (Fraction(45, 4), "11.25"),34 (Fraction(10), "10"),35 (90 * (ssqrt(2) - 1), "90 (sqrt2 - 1)"),36]37BASE_DEPTH = 83839def norm(z):40 return z[0] * z[0] + z[1] * z[1]4142def gmul(z, w):43 return (z[0] * w[0] - z[1] * w[1], z[0] * w[1] + z[1] * w[0])4445def gconj(z):46 return (z[0], -z[1])4748def gdivides(w, z):49 p = gmul(z, gconj(w))50 n = norm(w)51 return p[0] % n == 0 and p[1] % n == 05253def gquo(z, w):54 p = gmul(z, gconj(w))55 n = norm(w)56 return (p[0] // n, p[1] // n)5758def gnearest(z, w):59 p = gmul(z, gconj(w))60 n = norm(w)61 return ((2 * p[0] + n) // (2 * n), (2 * p[1] + n) // (2 * n))6263def ggcd(z, w):64 while w != (0, 0):65 q = gnearest(z, w)66 z, w = w, (z[0] - gmul(q, w)[0], z[1] - gmul(q, w)[1])67 return z6869def canon(z):70 if z == (0, 0):71 return z72 for c in (z, (-z[1], z[0]), (-z[0], -z[1]), (z[1], -z[0])):73 if c[0] > 0 and c[1] >= 0:74 return c75 return z7677def classes_up_to(n):78 out = []79 for a in range(1, isqrt(n) + 1):80 for b in range(0, isqrt(n - a * a) + 1):81 out.append((a, b))82 return sorted(out, key=lambda z: (norm(z), z))8384def circle_classes_jacobi(t):85 s = 086 j = 087 while 4 * j + 1 <= t:88 s += t // (4 * j + 1) - t // (4 * j + 3)89 j += 190 return s9192def circle_classes_direct(t):93 return len(classes_up_to(t))9495def gauss_primes_dividing(z):96 n = norm(z)97 out = []98 m = n99 q = 2100 while q * q <= m:101 if m % q == 0:102 while m % q == 0:103 m //= q104 out.append(q)105 q += 1106 if m > 1:107 out.append(m)108 primes = []109 for q in out:110 if q == 2:111 primes.append((1, 1))112 elif q % 4 == 3:113 primes.append((q, 0))114 else:115 r = 2116 a = None117 while a is None:118 if pow(r, (q - 1) // 2, q) == q - 1:119 a = pow(r, (q - 1) // 4, q)120 r += 1121 p = ggcd((q, 0), (a, 1))122 primes.append(canon(p))123 primes.append(canon(gconj(p)))124 return [p for p in primes if gdivides(p, z)]125126def gauss_totient(z):127 v = norm(z)128 if v == 1:129 return 1130 for p in gauss_primes_dividing(z):131 v = v // norm(p) * (norm(p) - 1)132 return v133134def totient_sum(n_max):135 return sum(gauss_totient(z) for z in classes_up_to(n_max))136137def cyc_reduce(terms):138 acc = {}139 for e, c in terms:140 k = (e % CYC,)141 acc[k] = acc.get(k, 0) + c142 acc = {k: v for k, v in acc.items() if v}143 if not acc:144 return Poly(0, X, domain="QQ")145 return Poly(acc, X, domain="QQ").rem(PHI360)146147RATIONAL_CACHE = {}148149def rational_cos_sin(d):150 d = d % CYC151 if d not in RATIONAL_CACHE:152 c = cyc_reduce([(d, 1), (-d, 1)])153 s = cyc_reduce([(90 + d, -1), (90 - d, 1)])154 RATIONAL_CACHE[d] = (c.degree() <= 0, s.degree() <= 0)155 return RATIONAL_CACHE[d]156157def rational_angle_degrees():158 return [d for d in range(CYC) if all(rational_cos_sin(d))]159160def spot_check_degrees():161 out = []162 for d in SPOT_DEGREES:163 c = minimal_polynomial(cos(pi * Rational(d, 180)), X, polys=True).degree()164 s = minimal_polynomial(sin(pi * Rational(d, 180)), X, polys=True).degree()165 out.append((d, c == 1, s == 1, rational_cos_sin(d)))166 return out167168def primes(k):169 out = []170 n = 2171 while len(out) < k:172 if all(n % p for p in out if p * p <= n):173 out.append(n)174 n += 1175 return out176177def prime_schedule(k):178 return list(zip(range(1, k + 1), [0] + primes(k - 1)))179180def dead_spin_pairs(schedule):181 hits = []182 for i in range(len(schedule)):183 for j in range(i + 1, len(schedule)):184 m, a = schedule[i]185 n, b = schedule[j]186 if all(rational_cos_sin((a - b) % CYC)):187 hits.append((m, a, n, b))188 return hits189190def shared_witness(m, a, n, b):191 g = gcd(m, n)192 q = ((a - b) % CYC) // 90193 v = [(1, 0), (0, 1), (-1, 0), (0, -1)][q]194 return (Fraction(v[0], g), Fraction(v[1], g))195196def rot(d):197 from math import cos as fc, radians, sin as fs198 return fc(radians(d)), fs(radians(d))199200def near_miss(schedule, box):201 best = None202 for i in range(len(schedule)):203 for j in range(i + 1, len(schedule)):204 m, a = schedule[i]205 n, b = schedule[j]206 if all(rational_cos_sin((a - b) % CYC)):207 continue208 c, s = rot(a - b)209 for u in range(-box, box + 1):210 for v in range(-box, box + 1):211 if u == 0 and v == 0:212 continue213 x = (c * u - s * v) * n / m214 y = (s * u + c * v) * n / m215 dx = abs(x - round(x))216 dy = abs(y - round(y))217 d = max(dx, dy)218 if best is None or d < best[0]:219 best = (d, m, a, n, b, u, v)220 return best221222def layer_nodes(z):223 a, b = z224 n = norm(z)225 out = set()226 span = a + b + 1227 for m in range(-span, span + 1):228 for k in range(-span, span + 1):229 x = Fraction(m * a + k * b, n)230 y = Fraction(k * a - m * b, n)231 out.add((x - int(x // 1), y - int(y // 1)))232 return out233234def literal_stack(n_max):235 hits = {}236 for z in classes_up_to(n_max):237 for p in layer_nodes(z):238 hits[p] = hits.get(p, 0) + 1239 return hits240241def reduced_denominator(node):242 x, y = node243 d = lcm(x.denominator, y.denominator)244 u = (x.numerator * (d // x.denominator), y.numerator * (d // y.denominator))245 g = ggcd(u, (d, 0))246 return canon(gquo((d, 0), g))247248def closed_brightness(n_max, d):249 return circle_classes_jacobi(n_max // norm(d))250251def pythagorean_hits(bound, reach):252 seen = set()253 for a in range(-reach, reach + 1):254 for b in range(-reach, reach + 1):255 if a == 0 and b == 0:256 continue257 n = a * a + b * b258 seen.add((Fraction(a * a - b * b, n), Fraction(2 * a * b, n)))259 direct = set()260 for r in range(1, bound + 1):261 for p in range(-r, r + 1):262 q2 = r * r - p * p263 s = isqrt(q2)264 if s * s == q2:265 for t in (s, -s):266 direct.add((Fraction(p, r), Fraction(t, r)))267 return len(direct), direct <= seen268269def base_c_overlap(c, box):270 best = None271 count = 0272 for m in range(-box, box + 1):273 for k in range(-box, box + 1):274 if m == 0 and k == 0:275 continue276 x = c[0] * m - c[1] * k277 y = c[1] * m + c[0] * k278 d = max(abs(x - round(x)), abs(y - round(y)))279 if d < 1e-12:280 count += 1281 if best is None or d < best:282 best = d283 return count, best284285def base_layer(c, k):286 z = (1, 0)287 for _ in range(k):288 z = gmul(z, c)289 return layer_nodes(canon(z))290291def unit_square_shares(g, a):292 c, s = rot(a)293 count = 0294 for k in range(-4 * g - 4, 4 * g + 5):295 for l in range(-4 * g - 4, 4 * g + 5):296 if k == 0 and l == 0:297 continue298 x = (c * k - s * l) / g299 y = (s * k + c * l) / g300 if 1e-9 < x < 1 - 1e-9 and 1e-9 < y < 1 - 1e-9:301 count += 1302 return count303304def is_whole_turn(ratio, d):305 v = d * ratio306 if isinstance(v, Fraction):307 return v.denominator == 1308 return bool(sympify(v).is_integer)309310def increment_period(ratio, layers):311 for q in range(1, layers + 1):312 if is_whole_turn(ratio, q):313 return q314 return None315316def increment_classes(theta_deg, layers):317 ratio = theta_deg / 90318 q = increment_period(ratio, layers)319 if q is None:320 sizes = [1] * layers321 else:322 sizes = [len(range(r, layers, q)) for r in range(q)]323 closed = sum(c * (c - 1) // 2 for c in sizes)324 pairwise = sum(325 1326 for j in range(layers)327 for k in range(j + 1, layers)328 if is_whole_turn(ratio, j - k)329 )330 niven = None331 if isinstance(theta_deg, Fraction) and theta_deg.denominator == 1:332 niven = sum(333 1334 for j in range(layers)335 for k in range(j + 1, layers)336 if all(rational_cos_sin(int(theta_deg) * (j - k)))337 )338 return q, sizes, closed, pairwise, niven339340def share_count_spread(g):341 return sorted({unit_square_shares(g, a) for a in range(1, 90)})342343def base_depth_check(c, depth):344 layers = [base_layer(c, k) for k in range(depth + 1)]345 nested = all(layers[k] <= layers[k + 1] for k in range(depth))346 bad = 0347 for p in layers[depth]:348 d = min(k for k in range(depth + 1) if p in layers[k])349 b = sum(1 for k in range(depth + 1) if p in layers[k])350 if b != depth + 1 - d:351 bad += 1352 return nested, len(layers[depth]), bad353354def main():355 rats = rational_angle_degrees()356 print("rational rotation degrees", *rats)357 print("rational rotation count", len(rats))358 spots = spot_check_degrees()359 print("spot degrees", len(spots))360 print("spot minpoly agrees with cyclotomic", all((c, s) == r for _, c, s, r in spots))361362 hits, covered = pythagorean_hits(60, 12)363 print("rational unit-circle points denominator <= 60", hits)364 print("all are w^2/N(w), Gaussian w in box 12", covered)365366 sched = prime_schedule(SCHEDULE_N)367 print("schedule layers", len(sched))368 print("schedule angles", *[a for _, a in sched])369 pairs = dead_spin_pairs(sched)370 print("schedule pairs", len(sched) * (len(sched) - 1) // 2)371 print("sharing pairs", len(pairs))372 for m, a, n, b in pairs:373 g = gcd(m, n)374 ok = (a - b) % 90 == 0 and n % g == 0 and m % g == 0375 w = shared_witness(m, a, n, b)376 print("share", m, a, n, b, "gcd", g, "exact", ok, "witness", str(w[0]), str(w[1]), "open square nodes", unit_square_shares(g, a))377 print("all sharing pairs congruent mod 90", all((a - b) % 90 == 0 for _, a, _, b in pairs))378 print("shared node window", "open unit square, origin excluded")379 print("shared lattice density per unit area", "g^2")380 for g in (2, 3, 7):381 print("share counts over whole degrees 1..89 at g", g, *share_count_spread(g))382 odds = len(range(1, ODD_N + 1, 2))383 print("increment schedule odd scales 1 to", ODD_N, "layers", odds)384 for theta, name in INCREMENTS:385 q, sizes, closed, pairwise, niven = increment_classes(theta, odds)386 print(387 "increment", name,388 "classes", "none" if q is None else q,389 "sizes", *(["all 1"] if q is None else sizes),390 "pairs", closed,391 "pairwise equal", closed == pairwise,392 "niven equal", "skipped" if niven is None else closed == niven,393 )394 d, m, a, n, b, u, v = near_miss(sched, NEAR_BOX)395 print("near miss box", NEAR_BOX)396 print("near miss min distance", "%.6f" % d, "at", m, a, n, b, u, v)397398 stack = literal_stack(STACK_N)399 print("stack norm bound", STACK_N)400 print("stack layers", len(classes_up_to(STACK_N)))401 print("stack nodes", len(stack))402 print("gaussian totient sum", totient_sum(STACK_N))403 print("nodes equal totient sum", len(stack) == totient_sum(STACK_N))404 bad = 0405 for p, b in stack.items():406 if closed_brightness(STACK_N, reduced_denominator(p)) != b:407 bad += 1408 print("brightness comparisons", len(stack))409 print("brightness mismatches", bad)410 print("brightest node", max(stack.values()), "at 0", stack[(Fraction(0), Fraction(0))])411 print("circle count jacobi equals direct", all(circle_classes_jacobi(t) == circle_classes_direct(t) for t in range(0, 401)))412 print("circle counts t=1..12", *[circle_classes_jacobi(t) for t in range(1, 13)])413 for n in COUNT_NS:414 s = totient_sum(n)415 print("node count N", n, "=", s, "ratio N^2", "%.6f" % (s / (n * n)))416 print("pi / (8 zeta(2) Catalan)", "%.6f" % float(N_(spi / (8 * szeta(2) * SCatalan))))417418 for c, name, depth in (((1, 1), "1+i", BASE_DEPTH), ((2, 1), "2+i", 4)):419 nested, size, bad = base_depth_check(c, depth)420 print("base", name, "depth", depth, "nested", nested, "deepest layer nodes", size, "address mismatches", bad)421 for c, name in (((1.0, 1.0), "1+i"), ((1.5, 0.5), "3/2+i/2"), ((2 ** 0.5 * __import__("math").cos(1.0), 2 ** 0.5 * __import__("math").sin(1.0)), "sqrt2 e^i")):422 count, best = base_c_overlap(c, BASE_BOX)423 print("base", name, "box", BASE_BOX, "overlaps", count, "min distance", "%.6g" % best)424425main()