gaussian_franel.py
13.5 kB · python · 440 lines
1from fractions import Fraction2from math import gcd, isqrt34import numpy as np5from sympy import Catalan, Poly, Symbol, cyclotomic_poly6from sympy import zeta as szeta78X = Symbol("x")9ZETA_K2 = float(szeta(2) * Catalan)10ZETA_2 = float(szeta(2))11CHECK_N = 5012CHECK_LAMBDA = 2013RESIDUE_N = 20014SAYOUS_TS = [2, 3, 4, 5, 6]15T2_NS = [20, 50]16T2_CUT = 20000017CLASSICAL_Q = 4018CLASSICAL_QS = [125, 250, 500, 1000, 2000, 4000, 8000]19METER_NS = [100, 250, 500, 1000, 2000, 4000, 8000, 16000, 32000, 64000]2021# GAUSSIAN ARITHMETIC2223def norm(z):24 return z[0] * z[0] + z[1] * z[1]2526def gmul(z, w):27 return (z[0] * w[0] - z[1] * w[1], z[0] * w[1] + z[1] * w[0])2829def gconj(z):30 return (z[0], -z[1])3132def gsub(z, w):33 return (z[0] - w[0], z[1] - w[1])3435def gnearest(z, w):36 p = gmul(z, gconj(w))37 n = norm(w)38 return ((2 * p[0] + n) // (2 * n), (2 * p[1] + n) // (2 * n))3940def ggcd(z, w):41 while w != (0, 0):42 q = gnearest(z, w)43 z, w = w, gsub(z, gmul(q, w))44 return z4546def gdivides(w, z):47 p = gmul(z, gconj(w))48 n = norm(w)49 return p[0] % n == 0 and p[1] % n == 05051def gquo(z, w):52 p = gmul(z, gconj(w))53 n = norm(w)54 return (p[0] // n, p[1] // n)5556def canon(z):57 if z == (0, 0):58 return z59 for c in (z, (-z[1], z[0]), (-z[0], -z[1]), (z[1], -z[0])):60 if c[0] > 0 and c[1] >= 0:61 return c62 return z6364def classes_up_to(n):65 out = []66 for a in range(1, isqrt(n) + 1):67 for b in range(0, isqrt(n - a * a) + 1):68 out.append((a, b))69 return sorted(out, key=lambda z: (norm(z), z))7071def elements_up_to(n):72 r = isqrt(n)73 out = []74 for a in range(-r, r + 1):75 for b in range(-r, r + 1):76 if 0 < a * a + b * b <= n:77 out.append((a, b))78 return sorted(out, key=lambda z: (norm(z), z))7980def is_rational_prime(n):81 if n < 2:82 return False83 i = 284 while i * i <= n:85 if n % i == 0:86 return False87 i += 188 return True8990def is_prime_class(z):91 n = norm(z)92 if is_rational_prime(n):93 return True94 r = isqrt(n)95 return r * r == n and r % 4 == 3 and is_rational_prime(r)9697def factor_class(z):98 out = []99 cur = canon(z)100 while norm(cur) > 1:101 for p in classes_up_to(norm(cur)):102 if norm(p) > 1 and is_prime_class(p) and gdivides(p, cur):103 out.append(p)104 cur = canon(gquo(cur, p))105 break106 else:107 raise ValueError108 return out109110def mobius_class(z):111 f = factor_class(z)112 return 0 if len(set(f)) != len(f) else (-1) ** len(f)113114def totient_class(z):115 v = norm(z)116 if v == 1:117 return 1118 for p in set(factor_class(z)):119 v = v // norm(p) * (norm(p) - 1)120 return v121122# RESIDUES MOD A GAUSSIAN INTEGER123124def residue_box(d):125 g = abs(gcd(d[0], d[1]))126 return g, norm(d) // g127128def residues(d):129 g, h = residue_box(d)130 return [(x, y) for x in range(g) for y in range(h)]131132def check_residues(n_max):133 bad_count = 0134 bad_distinct = 0135 bad_totient = 0136 for d in classes_up_to(n_max):137 n = norm(d)138 r = residues(d)139 if len(r) != n:140 bad_count += 1141 reps = {tuple(v % n for v in gmul(u, gconj(d))) for u in r}142 if len(reps) != n:143 bad_distinct += 1144 if len(coprime_residues(d)) != totient_class(d):145 bad_totient += 1146 return len(classes_up_to(n_max)), bad_count, bad_distinct, bad_totient147148def coprime_residues(d):149 return [u for u in residues(d) if norm(ggcd(u, d)) == 1]150151# EXACT SUMS OF ROOTS OF UNITY152153def roots_of_unity_sum(counts, n, phi):154 deg = len(phi) - 1155 c = list(counts)156 for k in range(n - 1, deg - 1, -1):157 f = c[k]158 if f:159 for j in range(deg + 1):160 c[k - j] -= f * phi[j]161 if any(c[1:deg]):162 return None163 return c[0]164165# THEOREM 1166167def ramanujan_exact(d, lam, phi_cache):168 n = norm(d)169 counts = [0] * n170 for u in coprime_residues(d):171 k = gmul(gmul(lam, u), gconj(d))[0] % n172 counts[k] += 1173 if n not in phi_cache:174 phi_cache[n] = [int(c) for c in Poly(cyclotomic_poly(n, X), X).all_coeffs()]175 return roots_of_unity_sum(counts, n, phi_cache[n])176177def ramanujan_formula(d, lam, mob):178 tot = 0179 for e in classes_up_to(norm(d)):180 if gdivides(e, d) and (lam == (0, 0) or gdivides(e, lam)):181 tot += mob[canon(gquo(d, e))] * norm(e)182 return tot183184def mertens_classes(t, mob):185 return sum(v for z, v in mob.items() if norm(z) <= t)186187def sum_formula(n_max, lam, mob):188 tot = 0189 for e in classes_up_to(n_max if lam == (0, 0) else min(n_max, norm(lam))):190 if lam == (0, 0) or gdivides(e, lam):191 tot += norm(e) * mertens_classes(n_max // norm(e), mob)192 return tot193194def node_set(n_max):195 out = []196 for d in classes_up_to(n_max):197 n = norm(d)198 for u in coprime_residues(d):199 w = gmul(u, gconj(d))200 out.append((Fraction(w[0] % n, n), Fraction(w[1] % n, n)))201 return out202203def torus_point(w, n):204 x, y = w[0] % n, w[1] % n205 g = gcd(gcd(x, y), n)206 return (x // g, y // g, n // g)207208def sayous_set(t):209 out = set()210 for q in elements_up_to(t * t):211 n = norm(q)212 cj = gconj(q)213 for x in range(n):214 for y in range(n):215 out.add(torus_point(gmul((x, y), cj), n))216 return out217218def node_points(n_max):219 out = set()220 for x, y in node_set(n_max):221 n = x.denominator * y.denominator // gcd(x.denominator, y.denominator)222 out.add(torus_point((x.numerator * (n // x.denominator), y.numerator * (n // y.denominator)), n))223 return out224225def check_sayous(ts):226 rows = []227 for t in ts:228 s = sayous_set(t)229 g = node_points(t * t)230 rows.append((t, t * t, len(s), len(g), s == g))231 return rows232233def literal_sum(nodes, lam):234 p, q = lam235 x = np.array([float(a) for a, _ in nodes])236 y = np.array([float(b) for _, b in nodes])237 return complex(np.sum(np.exp(2j * np.pi * (p * x - q * y))))238239def check_theorem_1(n_max, lam_bound):240 mob = {z: mobius_class(z) for z in classes_up_to(n_max)}241 phi_cache = {}242 nodes = node_set(n_max)243 bad_ram = 0244 bad_sum = 0245 bad_lit = 0246 worst = 0.0247 lams = elements_up_to(lam_bound)248 for lam in lams:249 tot = 0250 for d in classes_up_to(n_max):251 ex = ramanujan_exact(d, lam, phi_cache)252 fo = ramanujan_formula(d, lam, mob)253 if ex is None or ex != fo:254 bad_ram += 1255 tot += fo256 if tot != sum_formula(n_max, lam, mob):257 bad_sum += 1258 lit = literal_sum(nodes, lam)259 err = abs(lit - tot)260 worst = max(worst, err)261 if err > 1e-9:262 bad_lit += 1263 return len(nodes), len(lams), bad_ram, bad_sum, bad_lit, worst264265# THEOREM 2266267def gauss_weights(n_max, mob):268 out = []269 for e in classes_up_to(n_max):270 w = norm(e) * mertens_classes(n_max // norm(e), mob)271 if w:272 out.append((e, w))273 return out274275def franel_exact(n_max, mob):276 cls = classes_up_to(n_max)277 mert = {e: mertens_classes(n_max // norm(e), mob) for e in cls}278 tot = Fraction(0)279 for a in cls:280 if mert[a] == 0:281 continue282 for b in cls:283 if mert[b] == 0:284 continue285 g = norm(canon(ggcd(a, b)))286 tot += Fraction(g * g * mert[a] * mert[b], norm(a) * norm(b))287 return tot288289def lambda_side(n_max, cut, mob, a_arr):290 r = isqrt(cut)291 side = 2 * r + 1292 acc = np.zeros((side, side), dtype=np.int64)293 for e, w in gauss_weights(n_max, mob):294 lim = cut // norm(e)295 rr = isqrt(lim)296 s, t = np.meshgrid(np.arange(-rr, rr + 1), np.arange(-rr, rr + 1), indexing="ij")297 keep = (s * s + t * t <= lim) & ((s != 0) | (t != 0))298 s, t = s[keep], t[keep]299 xs = e[0] * s - e[1] * t300 ys = e[0] * t + e[1] * s301 np.add.at(acc, (xs + r, ys + r), w)302 ii, jj = np.meshgrid(np.arange(-r, r + 1), np.arange(-r, r + 1), indexing="ij")303 nn = ii * ii + jj * jj304 keep = (nn > 0) & (nn <= cut)305 vals = acc[keep].astype(np.float64)306 nrm = nn[keep].astype(np.float64)307 lhs = float(np.sum(vals * vals / (nrm * nrm)))308 weight = float(np.sum(1.0 / (nrm * nrm)))309 return lhs, weight310311def check_theorem_2(n_max, cut, mob, a_arr):312 m = sum(totient_class(d) for d in classes_up_to(n_max))313 exact = franel_exact(n_max, mob)314 rhs = 4.0 * ZETA_K2 * float(exact)315 lhs, weight = lambda_side(n_max, cut, mob, a_arr)316 tail = m * m * (4.0 * ZETA_K2 - weight)317 return m, exact, rhs, lhs, rhs - lhs, tail318319# THE SIEVES320321def chi4(n):322 return 0 if n % 2 == 0 else (1 if n % 4 == 1 else -1)323324def sieve(n_max, gaussian):325 a = np.zeros(n_max + 1, dtype=np.int64)326 if gaussian:327 for d in range(1, n_max + 1, 2):328 a[d::d] += chi4(d)329 else:330 a[1:] = 1331 mob = np.zeros(n_max + 1, dtype=np.int64)332 mob[1] = 1333 for n in range(1, n_max // 2 + 1):334 v = mob[n]335 if v:336 mob[2 * n :: n] -= a[2 : n_max // n + 1] * v337 j2 = np.zeros(n_max + 1, dtype=np.int64)338 ph = np.zeros(n_max + 1, dtype=np.int64)339 for d in range(1, n_max + 1):340 if a[d]:341 j2[d::d] += (d * d * a[d]) * mob[1 : n_max // d + 1]342 ph[d::d] += (d * a[d]) * mob[1 : n_max // d + 1]343 mg = np.zeros(n_max + 1, dtype=np.int64)344 mg[1:] = np.cumsum(mob[1:])345 return a, mob, mg, j2, np.cumsum(ph)346347def inner_h(t, a, mg):348 ms = np.arange(1, t + 1)349 return float(np.sum(a[1 : t + 1] * mg[t // ms] / ms))350351def franel_form(n, a, mg, j2):352 ns = np.arange(1, n + 1)353 ts = n // ns354 cache = {int(t): inner_h(int(t), a, mg) for t in np.unique(ts)}355 hv = np.array([cache[int(t)] for t in ts])356 return float(np.sum(j2[1 : n + 1] * hv * hv / (ns.astype(np.float64) ** 2)))357358# THE CLASSICAL CONTROL359360def farey_delta_square(q):361 nodes = sorted({Fraction(a, b) for b in range(1, q + 1) for a in range(1, b + 1)})362 m = len(nodes)363 return m, sum((r - Fraction(j + 1, m)) ** 2 for j, r in enumerate(nodes))364365def classical_exact(q):366 mu = {1: 1}367 for n in range(2, q + 1):368 mu[n] = -sum(mu[d] for d in range(1, n) if n % d == 0)369 mert = {n: sum(mu[k] for k in range(1, q // n + 1)) for n in range(1, q + 1)}370 tot = Fraction(0)371 for x in range(1, q + 1):372 for y in range(1, q + 1):373 g = gcd(x, y)374 tot += Fraction(g * g * mert[x] * mert[y], x * y)375 return tot376377# THE RUN378379def readout(x, a, mg):380 ms = np.arange(1, x + 1)381 return int(np.sum(a[1 : x + 1] * mg[x // ms]))382383def main():384 nodes, lams, bad_ram, bad_sum, bad_lit, worst = check_theorem_1(CHECK_N, CHECK_LAMBDA)385 ncls, bc, bd, bt = check_residues(RESIDUE_N)386 print(f"RESIDUE SYSTEMS to norm bound {RESIDUE_N}: {ncls} classes, {bc} wrong sizes, {bd} collisions mod d, {bt} totient mismatches")387 print("SAYOUS SET")388 for t, n, ns, ng, same in check_sayous(SAYOUS_TS):389 print(f" T = {t}, N = T^2 = {n}: literal G_T has {ns} points, node set has {ng}, equal is {same}")390 print(f"zeta_K(2) = zeta(2) * Catalan = {ZETA_K2:.6f}")391 print("THEOREM 1")392 print(f" nodes at N = {CHECK_N}: {nodes}; lambda with N(lambda) <= {CHECK_LAMBDA}: {lams}")393 print(f" exact Ramanujan mismatches: {bad_ram} of {lams * len(classes_up_to(CHECK_N))}")394 print(f" Mertens-form mismatches: {bad_sum} of {lams}")395 print(f" literal node-sum mismatches at 1e-9: {bad_lit} of {lams}, worst {worst:.3e}")396397 a_big, mob_big, mg_big, j2_big, ph_big = sieve(max(METER_NS), True)398 print("THEOREM 2")399 for n in T2_NS:400 mob = {z: mobius_class(z) for z in classes_up_to(n)}401 m, exact, rhs, lhs, diff, tail = check_theorem_2(n, T2_CUT, mob, a_big)402 ok = 0.0 <= diff <= tail403 print(f" N = {n}: m = {m}, F(N) = {exact} = {float(exact):.6f}")404 print(f" identity {rhs:.6f}, truncated Fourier side {lhs:.6f}, gap {diff:.6f}, tail bound {tail:.6f}, inside {ok}")405406 print("THE ZERO MODE")407 for n in T2_NS + [200]:408 mob = {z: mobius_class(z) for z in classes_up_to(n)}409 m = sum(totient_class(d) for d in classes_up_to(n))410 print(f" N = {n}: sum Phi(d) = {m}, sum N(e) M_G(N/N(e)) = {sum_formula(n, (0, 0), mob)}")411 print(f" sum_(N(a) <= x) M_G(x/N(a)) over x = 1..2000: always 1 is {all(readout(x, a_big, mg_big) == 1 for x in range(1, 2001))}")412413 print("CLASSICAL CONTROL")414 m_cl, s2 = farey_delta_square(CLASSICAL_Q)415 c_q = classical_exact(CLASSICAL_Q)416 print(f" Q = {CLASSICAL_Q}: m = {m_cl}, sum delta^2 = {float(s2):.10f}")417 print(f" C(Q) - 1 = 12 m sum delta^2 exactly: {c_q - 1 == 12 * m_cl * s2}")418 a_cl, mob_cl, mg_cl, j2_cl, ph_cl = sieve(max(METER_NS), False)419 print(f" sieve C(Q) = {franel_form(CLASSICAL_Q, a_cl, mg_cl, j2_cl):.9f} against exact {float(c_q):.9f}")420421 print(" Q | Phi(Q) | S2 * Q from the identity")422 for q in CLASSICAL_QS:423 cq = franel_form(q, a_cl, mg_cl, j2_cl)424 print(f" {q} | {int(ph_cl[q])} | {(cq - 1.0) * q / (12.0 * int(ph_cl[q])):.4f}")425426 print("THE METER")427 print(" N | m | F(N) | F(N)/N | D_2(N)^2 * N^3 | slope | M_G(N) | M_G(N)^2 / (zeta_K(2) F(N)) | classical C(N)/N")428 prev = None429 for n in METER_NS:430 f = franel_form(n, a_big, mg_big, j2_big)431 m = int(ph_big[n])432 d2sq = 4.0 * ZETA_K2 * f / (m * m)433 c = franel_form(n, a_cl, mg_cl, j2_cl)434 slope = "-" if prev is None else f"{np.log(f / prev[1]) / np.log(n / prev[0]):.4f}"435 mg = int(mg_big[n])436 print(f" {n} | {m} | {f:.4f} | {f / n:.6f} | {d2sq * n ** 3:.4f} | {slope} | {mg} | {mg * mg / (ZETA_K2 * f):.6f} | {c / n:.6f}")437 prev = (n, f)438439if __name__ == "__main__":440 main()