sides.py
17.4 kB · python · 441 lines
1import math2import sys3import time4from fractions import Fraction5from math import gcd, isqrt, log67import numpy as np89A030979_URL = "https://raw.githubusercontent.com/oeis/oeisdata/main/seq/A030/A030979.seq"1011# MEMBERSHIP1213def held(p, q, n):14 m = 2 * q15 a = p % m16 seen = set()17 while a not in seen:18 if a > q:19 return False20 seen.add(a)21 a = n * a % m22 return True2324def held_ifs(x, n):25 seen = set()26 while x not in seen:27 if x < 0 or x > 1:28 return False29 seen.add(x)30 y = n * x31 d = math.floor(y)32 if d % 2:33 if y != d:34 return False35 d -= 136 x = y - d37 return True3839def residues(p, q):40 return [r for r in range(1, 2 * q, 2) if held(p, q, r)]4142def parity_rule(p, q, r):43 s = p % q44 seen = set()45 while s not in seen:46 if s and (s - p) % 2:47 return False48 seen.add(s)49 s = r * s % q50 return True5152def is_prime(n):53 return n > 1 and all(n % d for d in range(2, isqrt(n) + 1))5455def phi(n):56 return sum(1 for a in range(1, n + 1) if gcd(a, n) == 1)5758def order(r, q):59 k, a = 1, r % q60 while a != 1:61 a = a * r % q62 k += 163 return k6465def primitive_root(q):66 return next(g for g in range(2, q) if order(g, q) == q - 1)6768def prime_share(p, q):69 m = q - 170 while m % 2 == 0:71 m //= 272 g = primitive_root(q)73 units = 074 for d in range(1, m + 1):75 if m % d:76 continue77 h = pow(g, (q - 1) // d, q)78 coset = {p * pow(h, i, q) % q for i in range(d)}79 if all((s - p) % 2 == 0 for s in coset):80 units += phi(d)81 return Fraction(1 + units, q)8283def least_period(bits):84 n = len(bits)85 return next(t for t in range(1, n + 1) if n % t == 0 and all(bits[i] == bits[i % t] for i in range(n)))8687def period(top):88 t0 = time.time()89 checks = law_bad = 090 for q in range(2, top + 1):91 for p in range(0, q + 1):92 if gcd(p, q) != 1:93 continue94 res = set(residues(p, q))95 last = 6 * q + 1 if q <= 40 else 2 * q + 196 for n in range(3, last + 1, 2):97 checks += 198 law_bad += held_ifs(Fraction(p, q), n) != (n % (2 * q) in res)99 print(f"law: orbit rule against the digit walk at every p/q in [0, 1], q <= {top}, odd sides 3..6q+1 at q <= 40 and 3..2q+1 above: {checks} checks, {law_bad} failures")100 count = sym_bad = flip_bad = rule_bad = unit_bad = prime_bad = low_bad = high_bad = 0101 best = []102 low_eq = []103 short = []104 primes_big = []105 for q in range(2, top + 1):106 for p in range(1, q):107 if gcd(p, q) != 1:108 continue109 count += 1110 res = residues(p, q)111 rs = set(res)112 share = Fraction(len(res), q)113 if res != residues(q - p, q):114 sym_bad += 1115 floor_ = Fraction(1, 2) if q == 2 else Fraction(2, q)116 ceil_ = Fraction(q + 1, 2 * q) if q % 2 else Fraction(1, 2)117 low_bad += share < floor_118 high_bad += share > ceil_119 if share == floor_:120 low_eq.append(q)121 best.append((share, p, q))122 bits = [held(p, q, 2 * i + 1) for i in range(q)]123 t = least_period(bits)124 if t < q:125 short.append((p, q, t))126 if q % 2 == 0:127 flip_bad += any(((q - r) % (2 * q) in rs) != (r in rs) for r in range(1, 2 * q, 2))128 else:129 rule_bad += any(parity_rule(p, q, r) != (r in rs) for r in range(1, 2 * q, 2))130 units = sum(1 for r in res if gcd(r, q) == 1)131 seen = {}132 for r in range(1, q):133 if gcd(r, q) != 1:134 continue135 g = frozenset(pow(r, i, q) for i in range(order(r, q)))136 seen[g] = all((p * s % q - p) % 2 == 0 for s in g)137 formula = sum(phi(len(g)) for g, ok in seen.items() if ok)138 unit_bad += formula != units139 if is_prime(q):140 ps = prime_share(p, q)141 prime_bad += ps != share142 if p == 1 and ps > Fraction(2, q):143 primes_big.append(q)144 best.sort(key=lambda r: (-r[0], r[2], r[1]))145 print(f"fractions p/q in (0,1), q <= {top}: {count}")146 print(f"symmetry p -> q-p: {sym_bad} failures; even q, r -> q-r: {flip_bad} failures")147 print(f"odd q parity rule mod q: {rule_bad} failures; unit formula over cyclic subgroups: {unit_bad} failures; prime formula: {prime_bad} failures")148 print(f"lower bound 2/q (1/2 at q=2): {low_bad} failures, attained at {len(low_eq)} fractions, least q {sorted(set(low_eq))[:12]}")149 print(f"upper bound (q+1)/(2q) odd q, 1/2 even q: {high_bad} failures")150 print("largest shares:", ", ".join(f"{p}/{q} {s}" for s, p, q in best[:8]))151 above = [(p, q, s) for s, p, q in best if q > 3 and s > Fraction(1, 2)]152 print(f"shares above 1/2 with q > 3: {len(above)}")153 odd = max((s, -q, p) for s, p, q in best if q % 2 and q >= 5)154 half = sorted({q for s, p, q in best if q % 2 == 0 and s == Fraction(1, 2)})155 print(f"largest share at odd q >= 5: {odd[2]}/{-odd[1]} {odd[0]}; even q with some share 1/2: {len(half)}, first {half[:10]}")156 print(f"least period in n = (N-1)/2 below q: {len(short)} fractions, first {short[:6]}")157 print(f"odd primes q <= {top} with share(1/q) > 2/q: {primes_big}")158 print(f"runtime {time.time() - t0:.2f} s")159160# EISENSTEIN161162def eisenstein(top):163 t0 = time.time()164 tests = bad = full_match = 0165 for q in range(3, top + 1):166 if not is_prime(q):167 continue168 for p in range(1, q):169 c = sum(1 for n in range(3, q, 2) if (p * n // q) % 2)170 a = p if p % 2 else q - p171 leg = pow(a, (q - 1) // 2, q)172 tests += 1173 bad += (c % 2 == 1) != (leg == q - 1)174 full = sum(p * n // q for n in range(1, 2 * q, 2))175 assert full == (2 * p - 1) * (q - 1) // 2 + p176 full_match += (full % 2 == 1) == (leg == q - 1)177 print(f"odd primes q <= {top}, every p: parity of #(odd 3 <= N < q, first base-N digit of p/q odd) against (a/q), a = p or q-p odd")178 print(f"tests {tests}, failures {bad}")179 print(f"sum of floor(pN/q) over odd N < 2q equals (2p-1)(q-1)/2 + p at all {tests}; its parity agrees with the symbol at {full_match}")180 print(f"runtime {time.time() - t0:.2f} s")181182# INTEGER COUNT183184def in_k(k, n):185 h = (n - 1) // 2186 while k:187 if k % n > h:188 return False189 k //= n190 return True191192def counts(top):193 t = isqrt(2 * top) + 1194 c = np.zeros(top + 1, dtype=np.int64)195 for n in range(3, t + 1, 2):196 h = (n - 1) // 2197 digits = np.arange(h + 1, dtype=np.int64)198 vals = digits[digits <= top]199 pw = n200 while pw <= top:201 step = digits * pw202 step = step[step <= top]203 vals = (vals[None, :] + step[:, None]).ravel()204 vals = vals[vals <= top]205 pw *= n206 vals = vals[vals >= (n + 1) // 2]207 c += np.bincount(vals, minlength=top + 1)208 diff = np.zeros(top + 2, dtype=np.int64)209 first = t + 1 if (t + 1) % 2 else t + 2210 j = 2211 while j * first <= 2 * top:212 ns = np.arange(first, 2 * top // j + 1, 2, dtype=np.int64)213 lo = j * ns // 2214 hi = np.minimum(((j + 1) * ns - 1) // 2, top)215 diff += np.bincount(lo, minlength=top + 2)216 diff -= np.bincount(hi + 1, minlength=top + 2)217 j += 2218 return c + np.cumsum(diff)[: top + 1]219220def count(top):221 t0 = time.time()222 c = counts(top)223 t1 = time.time()224 small = 3000225 bad = sum(1 for k in range(1, small + 1) if c[k] != sum(1 for n in range(3, 2 * k + 1, 2) if in_k(k, n)))226 print(f"c(k) = #(odd 3 <= N <= 2k : 2k in Z_N) for all k <= {top}: {t1 - t0:.2f} s; direct digit test at k <= {small}: {bad} failures")227 k = np.arange(top + 1, dtype=np.float64)228 err = c - (1 - log(2)) * k229 ratio = np.zeros_like(err)230 ratio[1:] = err[1:] / np.sqrt(k[1:])231 i_max = int(np.argmax(np.abs(ratio[1:]))) + 1232 print(f"max over 1 <= k <= {top} of |c(k) - (1 - log 2) k| / sqrt(k): {math.ceil(abs(ratio[i_max]) * 10**6) / 10**6:.6f} at k = {i_max}")233 for lo in [10**2, 10**3, 10**4, 10**5]:234 hi = min(10 * lo, top)235 seg = ratio[lo:hi + 1]236 print(f" k in [{lo}, {hi}]: err/sqrt(k) min {math.floor(seg.min() * 10**6) / 10**6:.6f} max {math.ceil(seg.max() * 10**6) / 10**6:.6f} mean {seg.mean():.6f}")237 for kk in [10**3, 10**4, 10**5, 10**6]:238 if kk <= top:239 print(f" c({kk}) = {c[kk]}, (1 - log 2) k = {(1 - log(2)) * kk:.3f}, err/sqrt(k) = {ratio[kk]:.6f}")240 bound = np.sqrt(2 * k[1:]) + 1241 print(f"largest |err| / (sqrt(2k) + 1) over 1 <= k <= {top}: {math.ceil(np.max(np.abs(err[1:]) / bound) * 10**6) / 10**6:.6f}")242 print(f"runtime {time.time() - t0:.2f} s")243244def fast(k):245 x = 2 * k246 s = isqrt(x)247 b = sum(1 for n in range(3, s + 1, 2) if in_k(k, n))248 a, j = 0, 2249 while x // j > s:250 lo, hi = max(x // (j + 1), s), x // j251 a += (hi + 1) // 2 - (lo + 1) // 2252 j += 2253 return a + b, b254255def second(top):256 t0 = time.time()257 bad = sum(1 for k in range(1, 3001) if fast(k)[0] != sum(1 for n in range(3, 2 * k + 1, 2) if in_k(k, n)))258 print(f"block count against the direct digit test at k <= 3000: {bad} failures")259 zeta_half = -1.4603545088095868260 kappa = -(2 - math.sqrt(2)) * zeta_half / 4261 beta = (math.sqrt(2) + (2 - math.sqrt(2)) * zeta_half) / 4262 print(f"kappa = -(2 - sqrt 2) zeta(1/2)/4 = {kappa:.6f}; beta = (sqrt 2 + (2 - sqrt 2) zeta(1/2))/4 = {beta:.6f}")263 rng = np.random.default_rng(376)264 for e in range(6, top + 1):265 ks = [int(v) for v in rng.integers(10**e, 2 * 10**e, size=60)]266 r = []267 rb = []268 for k in ks:269 v, b = fast(k)270 r.append((v - (1 - log(2)) * k) / math.sqrt(k))271 rb.append(b / math.sqrt(k))272 print(f" 60 k in [10^{e}, 2 10^{e}): err/sqrt(k) mean {np.mean(r):.6f} min {min(r):.6f} max {max(r):.6f}; small-side part / sqrt(k) mean {np.mean(rb):.6f}")273 print(f"runtime {time.time() - t0:.2f} s")274275# FAMILY276277def carries_ok(k, p, a):278 c = 0279 i = 0280 while k or c:281 s = 2 * (k % p) + c282 c = 1 if s >= p else 0283 if c and i % a == a - 1:284 return False285 k //= p286 i += 1287 return True288289def family():290 t0 = time.time()291 meet = [m for m in range(0, 3000) if all(m % 2 == 0 and in_k(m // 2, n) for n in range(3, m + 3, 2))]292 print(f"integers below 3000 held by every odd side 3 <= N <= m+1: {meet}")293 miss = [(p, q) for q in range(2, 61) for p in range(1, q) if gcd(p, q) == 1 and held(p, q, 2 * q - 1)]294 print(f"p/q in (0,1), q <= 60, held at side 2q-1: {len(miss)}")295 sub_int = all(in_k(k, n ** e) for n in (3, 5, 7) for e in (2, 3) for k in range(0, 10**5) if in_k(k, n))296 sub_rat = all(held(p, q, pow(n, e, 2 * q)) for q in range(2, 61) for p in range(0, q + 1) if gcd(p, q) == 1 for n in range(3, 2 * q + 1, 2) for e in (2, 3, 4) if held(p, q, n))297 print(f"E_N inside E_(N^e): integers k < 10^5 at N = 3, 5, 7, e = 2, 3: {sub_int}; rationals q <= 60, e = 2, 3, 4: {sub_rat}")298 strict = all(not held(n + 1, n ** e, n) and held(n + 1, n ** e, n ** e) and not in_k((n + 1) // 2, n) and in_k((n + 1) // 2, n ** e) for n in range(3, 40, 2) for e in (2, 3, 4))299 print(f"strict at e >= 2: (N+1)/N^e and N+1 lie at side N^e and not at side N, odd N < 40, e = 2, 3, 4: {strict}")300 kum = all((math.comb(2 * k, k) % p != 0) == in_k(k, p) for p in (3, 5, 7, 11, 13, 17, 19, 23) for k in range(0, 1500))301 print(f"Kummer: K_p = (k : p does not divide C(2k,k)) at odd primes p <= 23, k < 1500: {kum}")302 pw = all(in_k(k, p ** a) == carries_ok(k, p, a) for p, a in ((3, 2), (3, 3), (5, 2), (7, 2)) for k in range(0, 10**5))303 print(f"K_(p^a) = no carry of k + k in base p out of a position = a-1 mod a, (p,a) in (3,2),(3,3),(5,2),(7,2), k < 10^5: {pw}")304 v3 = next(k for k in range(10**5) if in_k(k, 9) and (math.comb(2 * k, k) % 9 == 0))305 print(f"least k in K_9 with 9 | C(2k,k): {v3}; least k with 9 not dividing C(2k,k) outside K_9: {next(k for k in range(10**5) if not in_k(k, 9) and math.comb(2 * k, k) % 9)}")306 a = next(k for k in range(10**5) if in_k(k, 3) and in_k(k, 5) and not in_k(k, 15))307 b = next(k for k in range(10**5) if in_k(k, 15) and not in_k(k, 3))308 cc = next(k for k in range(10**5) if in_k(k, 15) and not in_k(k, 5))309 d = next(k for k in range(10**5) if in_k(k, 15) and not (math.comb(2 * k, k) % 15))310 print(f"least k in K_3 cap K_5 outside K_15: {a}; least in K_15 outside K_3: {b}; outside K_5: {cc}; least k in K_15 with 15 | C(2k,k): {d}")311 print(f"runtime {time.time() - t0:.2f} s")312313# FINITE INTERSECTIONS314315def next_in(n, lo):316 h = (n - 1) // 2317 while True:318 x, i, bad = lo, 0, -1319 while x:320 if x % n > h:321 bad = i322 x //= n323 i += 1324 if bad < 0:325 return lo326 p = n ** (bad + 1)327 lo = (lo // p + 1) * p328329def common(sides, top):330 b, rest = sides[0], sides[1:]331 h = (b - 1) // 2332 depth = 0333 while b ** depth <= top:334 depth += 1335 pw = [b ** i for i in range(depth + 1)]336 span = [h * (pw[i] - 1) // (b - 1) for i in range(depth + 1)]337 out, nodes, stack = [], 0, [(0, depth)]338 while stack:339 base, j = stack.pop()340 nodes += 1341 if base > top:342 continue343 hi = min(base + span[j], top)344 if any(next_in(n, base) > hi for n in rest):345 continue346 if j == 0:347 out.append(base)348 continue349 for d in range(h, -1, -1):350 stack.append((base + d * pw[j - 1], j - 1))351 return sorted(out), nodes352353def a030979(path):354 terms = []355 with open(path) as f:356 for line in f:357 if line.startswith("A030979 "):358 return [int(t) for t in line.split(",")[1:] if t.strip()]359 if line[:2] in ("%S", "%T", "%U") and "A030979" in line:360 terms += [int(t) for t in line.split(None, 2)[2].split(",") if t.strip()]361 return terms362363def inter(x, path=None):364 t0 = time.time()365 top = (x - 1) // 2366 small = 10**6367 for sides in ((3, 5), (3, 5, 7), (3, 5, 7, 11), (3, 5, 15)):368 brute = [k for k in range(small) if all(in_k(k, n) for n in sides)]369 got, _ = common(sides, small - 1)370 assert got == brute, sides371 print(f"control: pruned walk equals the direct digit test below k = {small} at four side sets")372 oeis = a030979(path) if path else None373 if not oeis:374 print(f"A030979 comparison skipped: pass a copy of {A030979_URL} as the second argument")375 for sides in ((3, 5), (3, 5, 15), (3, 5, 7), (3, 5, 7, 9), (3, 5, 7, 11), (3, 5, 7, 13), (3, 5, 7, 15), (3, 5, 7, 11, 13), (3, 5, 7, 11, 15)):376 t = time.time()377 got, nodes = common(sides, top)378 ints = [2 * k for k in got]379 tail = f": {ints}" if len(ints) <= 20 else ""380 print(f"sides {sides}: {len(ints)} integers below {x:.0e}, {nodes} nodes, {time.time() - t:.2f} s{tail}")381 if sides == (3, 5, 7) and oeis:382 print(f" equals twice the A030979 terms below {x:.0e}: {got == [k for k in oeis if k <= top]}")383 print(f"runtime {time.time() - t0:.2f} s")384385def deep(e, sides):386 t0 = time.time()387 got, nodes = common(sides, (10**e - 1) // 2)388 print(f"sides {sides}: integers below 10^{e}: {[2 * k for k in got]}, {nodes} nodes, {time.time() - t0:.2f} s")389390# IRRATIONALS391392def weyl(top, levels):393 t0 = time.time()394 scale = 10**90395 import mpmath396 mpmath.mp.dps = 110397 xs = {398 "sqrt(2) - 1": isqrt(2 * scale * scale) - scale,399 "(sqrt(5) - 1)/2": (isqrt(5 * scale * scale) - scale) // 2,400 "2^(1/3) - 1": int(mpmath.floor(mpmath.cbrt(2) * scale)) - scale,401 "pi - 3": int(mpmath.floor(mpmath.pi * scale)) - 3 * scale,402 "e - 2": int(mpmath.floor(mpmath.e * scale)) - 2 * scale,403 }404 sides = range(3, 2 * top + 2, 2)405 print(f"share of odd sides 3 <= N <= {2 * top + 1} holding x to level L (first L digits even), against 2^-L")406 for name, v in xs.items():407 hits = [0] * (levels + 1)408 for n in sides:409 a = v410 for lev in range(1, levels + 1):411 a = a * n % (2 * scale)412 if a >= scale:413 break414 hits[lev] += 1415 row = " ".join(f"{hits[lev] / len(sides) * 2**lev:.4f}" for lev in range(1, levels + 1))416 print(f" {name}: share times 2^L at L = 1..{levels}: {row}")417 print(f"runtime {time.time() - t0:.2f} s")418419if __name__ == "__main__":420 verb = sys.argv[1] if len(sys.argv) > 1 else ""421 words = sys.argv[2:]422 args = [int(float(a)) for a in words if a.replace(".", "").replace("e", "").isdigit()]423 if verb == "period":424 period(*(args or [200]))425 elif verb == "eisenstein":426 eisenstein(*(args or [200]))427 elif verb == "count":428 count(*(args or [10**6]))429 elif verb == "family":430 family()431 elif verb == "inter":432 files = [a for a in words if not a.replace(".", "").replace("e", "").isdigit()]433 inter(args[0] if args else 10**12, files[0] if files else None)434 elif verb == "second":435 second(*(args or [11]))436 elif verb == "deep":437 deep(args[0] if args else 30, (3, 5, 7, 11))438 elif verb == "weyl":439 weyl(*(args or [10**5, 6]))440 else:441 raise SystemExit("verbs: period Q, eisenstein Q, count K, second E, family, inter X [FILE], deep D, weyl M L")