stack_algebra.py
14.1 kB · python · 374 lines
1import math2from fractions import Fraction34import numpy as np56N_CONV = 607N_SELECT = [30, 61, 200, 501]8HARMONIC_NS = [1000, 4000, 16000]9S_VALUES = [3, 4, 5]10PRIME_LS = [5, 10, 100, 1000]11ODD_LS = [10, 100, 1000, 4000]12SQUAREFREE_LS = [10, 50, 200, 1000]13DAVENPORT_NS = [10 ** 3, 10 ** 4, 10 ** 5]14DAVENPORT_X = [Fraction(1, 7), Fraction(1, 4), Fraction(1, 3), Fraction(2, 5), Fraction(3, 11)]1516# SIEVES1718def mobius_table(n):19 mu = np.ones(n + 1, dtype=np.int64)20 composite = np.zeros(n + 1, dtype=bool)21 for p in range(2, n + 1):22 if composite[p]:23 continue24 composite[p::p] = True25 mu[p::p] *= -126 square = p * p27 if square <= n:28 mu[square::square] = 029 mu[0] = 030 return mu3132def primes_upto(n):33 flag = np.ones(n + 1, dtype=bool)34 flag[:2] = False35 for p in range(2, int(n ** 0.5) + 1):36 if flag[p]:37 flag[p * p :: p] = False38 return np.nonzero(flag)[0]3940def first_primes(count):41 limit = 6442 while True:43 found = primes_upto(limit)44 if found.size >= count:45 return [int(p) for p in found[:count]]46 limit *= 24748def prime_powers_upto(n):49 out = {}50 for p in primes_upto(n):51 p = int(p)52 power = p53 exponent = 154 while power <= n:55 out[power] = (p, exponent)56 power *= p57 exponent += 158 return out5960def totients(n):61 phi = np.arange(n + 1, dtype=np.int64)62 for p in range(2, n + 1):63 if phi[p] == p:64 phi[p::p] -= phi[p::p] // p65 return phi6667# DIRICHLET CONVOLUTION6869def convolve(u, v, n):70 out = [Fraction(0)] * (n + 1)71 for k in range(1, n + 1):72 if u[k] == 0:73 continue74 for m in range(1, n // k + 1):75 out[k * m] += u[k] * v[m]76 return out7778def literal_stack(w, n):79 nodes = {}80 for scale in range(1, n + 1):81 if w[scale] == 0:82 continue83 for j in range(1, scale + 1):84 g = math.gcd(j, scale)85 key = (j // g, scale // g)86 nodes[key] = nodes.get(key, Fraction(0)) + w[scale]87 return nodes8889def literal_double_stack(u, v, n):90 nodes = {}91 for k in range(1, n + 1):92 if u[k] == 0:93 continue94 for m in range(1, n // k + 1):95 weight = u[k] * v[m]96 if weight == 0:97 continue98 scale = k * m99 for j in range(1, scale + 1):100 g = math.gcd(j, scale)101 key = (j // g, scale // g)102 nodes[key] = nodes.get(key, Fraction(0)) + weight103 return nodes104105def nested_stack(u, v, n):106 nodes = {}107 for k in range(1, n + 1):108 if u[k] == 0:109 continue110 for m in range(1, n // k + 1):111 weight = u[k] * v[m]112 if weight == 0:113 continue114 for i in range(1, k + 1):115 for j in range(1, m + 1):116 num = (i - 1) * m + j117 den = k * m118 g = math.gcd(num, den)119 key = (num // g, den // g)120 nodes[key] = nodes.get(key, Fraction(0)) + weight121 return nodes122123def rectangular_weights(u, v, outer, inner):124 out = {}125 for k in range(1, outer + 1):126 for m in range(1, inner + 1):127 out[k * m] = out.get(k * m, Fraction(0)) + u[k] * v[m]128 return out129130def node_formula(w, n, b):131 return sum((w[k * b] for k in range(1, n // b + 1)), Fraction(0))132133def convolution_check():134 n = N_CONV135 mu = mobius_table(n)136 one = [Fraction(0)] + [Fraction(1)] * n137 mob = [Fraction(0)] + [Fraction(int(mu[k])) for k in range(1, n + 1)]138 inv = [Fraction(0)] + [Fraction(1, k) for k in range(1, n + 1)]139 pairs = [("1 * 1", one, one), ("1 * mu", one, mob), ("mu * mu", mob, mob), ("1 * n^-1", one, inv)]140 print("CONVOLUTION hyperbolic cut k*m <= N, N =", n)141 for name, u, v in pairs:142 w = convolve(u, v, n)143 left = literal_double_stack(u, v, n)144 right = literal_stack(w, n)145 nested = nested_stack(u, v, n)146 keys = set(left) | set(right) | set(nested)147 bad = sum(1 for key in keys if left.get(key, Fraction(0)) != right.get(key, Fraction(0)))148 geo = sum(1 for key in keys if nested.get(key, Fraction(0)) != right.get(key, Fraction(0)))149 form = sum(1 for (a, b) in keys if right.get((a, b), Fraction(0)) != node_formula(w, n, b))150 rect = rectangular_weights(u, v, 24, 40)151 cut = sum(1 for m in range(1, 25) if rect.get(m, Fraction(0)) != w[m])152 above = sum(1 for m in range(25, 41) if rect.get(m, Fraction(0)) != w[m])153 print(" ", name, "nodes", len(keys), "scaled-copy mismatches", geo, "hyperbolic mismatches", bad, "formula mismatches", form)154 print(" ", name, "rectangular cut 24 x 40: mismatches at m <= 24:", cut, " at 25 <= m <= 40:", above)155 e = convolve(one, mob, n)156 print(" 1 * mu identity e(1) =", e[1], " nonzero beyond 1:", sum(1 for k in range(2, n + 1) if e[k] != 0))157 d = convolve(one, one, n)158 print(" 1 * 1 is d(n), first ten", [int(d[k]) for k in range(1, 11)])159 print(" scale weight W(M) = (u*v)(M) for every M <= N under the hyperbolic cut")160161# SELECTIONS162163def selection_literal(weight, n, b):164 return sum(1 for k in range(1, n // b + 1) if weight(k * b))165166def even_form(n, b):167 return n // (b if b % 2 == 0 else 2 * b)168169def odd_form(n, b):170 return 0 if b % 2 == 0 else (n // b + 1) // 2171172def prime_form(n, b, prime_set, pi_n):173 if b == 1:174 return pi_n175 return 1 if b in prime_set else 0176177def coprime_count(y, b_divisors_mu):178 return sum(m * (y // e) for e, m in b_divisors_mu)179180def squarefree_form(n, b, mu):181 if mu[b] == 0:182 return 0183 divisors = [(e, int(mu[e])) for e in range(1, b + 1) if b % e == 0 and mu[e] != 0]184 x = n // b185 total = 0186 d = 1187 while d * d <= x:188 if mu[d] != 0 and math.gcd(d, b) == 1:189 total += int(mu[d]) * coprime_count(x // (d * d), divisors)190 d += 1191 return total192193def prime_power_form(n, b, powers):194 if b == 1:195 return sum(int(math.log(n, p) + 1e-9) for p in primes_upto(n))196 if b not in powers:197 return 0198 p, i = powers[b]199 top = int(math.log(n, p) + 1e-9)200 return max(0, top - i + 1)201202def selection_closed_forms():203 print("SELECTIONS closed form against literal stacking, all b <= N")204 for n in N_SELECT:205 mu = mobius_table(n)206 prime_set = set(int(p) for p in primes_upto(n))207 powers = prime_powers_upto(n)208 pi_n = len(prime_set)209 rows = []210 checks = [211 ("evens", lambda m: m % 2 == 0, lambda b: even_form(n, b)),212 ("odds", lambda m: m % 2 == 1, lambda b: odd_form(n, b)),213 ("primes", lambda m: m in prime_set, lambda b: prime_form(n, b, prime_set, pi_n)),214 ("squarefree", lambda m: mu[m] != 0, lambda b: squarefree_form(n, b, mu)),215 ("prime powers", lambda m: m in powers, lambda b: prime_power_form(n, b, powers)),216 ]217 for name, weight, form in checks:218 bad = sum(1 for b in range(1, n + 1) if selection_literal(weight, n, b) != form(b))219 rows.append("%s %d" % (name, bad))220 print(" N =", n, " mismatches:", ", ".join(rows))221 n = N_SELECT[-1]222 prime_set = set(int(p) for p in primes_upto(n))223 lit = sum(1 for b in range(1, n + 1) if selection_literal(lambda m: m in prime_set, n, b) > 0)224 print(" primes-only line stack at N =", n, "lights", lit, "denominators: b = 1 and the", len(prime_set), "primes")225226# HARMONIC STACK227228def zeta_em(s, terms=200000):229 k = np.arange(1, terms + 1, dtype=np.float64)230 total = float(np.sum(k ** (-float(s))))231 m = float(terms)232 return total + m ** (1 - s) / (s - 1) - 0.5 * m ** (-s) + s * m ** (-s - 1) / 12.0233234def harmonic_stack():235 print("HARMONIC STACK weights n^-s, node a/b reads b^-s H_s(floor(N/b))")236 for s in S_VALUES:237 target = zeta_em(s - 1)238 row = []239 for n in HARMONIC_NS:240 k = np.arange(1, n + 1, dtype=np.float64)241 partial = np.concatenate((np.zeros(1), np.cumsum(k ** (-float(s)))))242 phi = totients(n)243 b = np.arange(1, n + 1, dtype=np.float64)244 mass = float(np.sum(phi[1:] * b ** (-float(s)) * partial[n // np.arange(1, n + 1)]))245 row.append("N=%d %.6f" % (n, mass))246 print(" s =", s, " zeta(s-1) = %.6f" % target, " total node mass", ", ".join(row))247 n = HARMONIC_NS[-1]248 s = S_VALUES[0]249 k = np.arange(1, n + 1, dtype=np.float64)250 partial = np.concatenate((np.zeros(1), np.cumsum(k ** (-float(s)))))251 z = zeta_em(s)252 for b in (1, 2, 5, 17):253 node = float(partial[n // b]) * b ** (-float(s))254 literal = float(np.sum((b * np.arange(1, n // b + 1, dtype=np.float64)) ** (-float(s))))255 print(" s = %d b = %2d node %.9f literal %.9f zeta(s)/b^s %.9f" % (s, b, node, literal, z / b ** s))256257# CARPET LAYERS258259def sign_at(n, i, length):260 return 1 if ((n * i) // length) % 2 == 0 else -1261262def mean_sign(m):263 return Fraction(sum(sign_at(m, i, m) for i in range(m)), m)264265def pair_integral(m, n):266 length = m * n // math.gcd(m, n)267 return Fraction(sum(sign_at(m, i, length) * sign_at(n, i, length) for i in range(length)), length)268269def carpet_exact(m, n):270 am, an = mean_sign(m), mean_sign(n)271 a = pair_integral(m, n)272 mu_m, mu_n = (1 - am) / 2, (1 - an) / 2273 nu = (1 - am - an + a) / 4274 return nu * nu - (mu_m * mu_n) ** 2275276def carpet_var(n):277 q = Fraction(n - 1, 2 * n)278 return q ** 2 - q ** 4279280def carpet_cov(m, n):281 d = math.gcd(m, n)282 return Fraction((d * d - 1) * (2 * (m - 1) * (n - 1) + d * d - 1), 16 * m * m * n * n)283284def pair_sum(layers):285 a = np.array(layers, dtype=np.int64)286 m = a.astype(np.float64)287 var = float(np.sum((m - 1) ** 2 / (4 * m ** 2) - (m - 1) ** 4 / (16 * m ** 4)))288 off = 0.0289 hits = 0290 step = 512291 for start in range(0, a.size, step):292 block = a[start : start + step]293 d = np.gcd.outer(block, a).astype(np.float64)294 x = block.astype(np.float64)[:, None]295 y = m[None, :]296 cov = (d * d - 1) * (2 * (x - 1) * (y - 1) + d * d - 1) / (16 * x * x * y * y)297 cov[np.arange(block.size), np.arange(start, start + block.size)] = 0.0298 off += float(np.sum(cov))299 hits += int(np.count_nonzero(d > 1))300 return var, off, (hits - a.size) // 2301302def prime_carpet_variance():303 print("CARPET LAYER LAW closed form against exact rational integration")304 for m, n in ((3, 5), (5, 7), (3, 9), (15, 21), (5, 15)):305 print(" (%d,%d) exact %s closed %s equal %s" % (m, n, carpet_exact(m, n), carpet_cov(m, n), carpet_exact(m, n) == carpet_cov(m, n)))306 twelve = first_primes(13)[1:]307 bad = sum(1 for i, p in enumerate(twelve) for q in twelve[i + 1 :] if carpet_exact(p, q) != 0)308 print(" odd prime pairs among", twelve, "with nonzero exact covariance:", bad)309 print("PRIME CARPET STACK L*Var of the L-layer mean, odd primes")310 for count in PRIME_LS:311 ps = first_primes(count + 1)[1:]312 total = sum((carpet_var(p) for p in ps), Fraction(0))313 var, off, sharing = pair_sum(ps)314 exact = str(total / count) if count <= 5 else "-"315 print(" L = %4d pairs sharing a factor %d covariance sum %.1f L*Var = %.10f independent %.10f ratio %.10f" % (count, sharing, off, float(total / count), var / count, (var + off) / var))316 if exact != "-":317 print(" exact L*Var = %s" % exact)318 print(" limit 3/16 = %.10f c = sqrt(3)/4 = %.7f c factor over independent exactly 1 at every L" % (3 / 16, math.sqrt(3) / 4))319 bad = sum(1 for p in first_primes(1001)[1:] if 16 * carpet_var(p) * p ** 4 != 3 * p ** 4 - 4 * p ** 3 - 2 * p ** 2 + 4 * p - 1)320 print(" 16 p^4 Var_p = 3p^4 - 4p^3 - 2p^2 + 4p - 1 breaches over the first 1000 odd primes:", bad)321 print("ODD CARPET STACK the tree's full odd stack for comparison")322 for count in ODD_LS:323 layers = list(range(3, 2 * count + 3, 2))324 var, off, _ = pair_sum(layers)325 actual = (var + off) / count326 indep = var / count327 print(" L = %4d L*Var = %.7f independent %.7f ratio %.6f c factor %.6f" % (count, actual, indep, actual / indep, math.sqrt(actual / indep)))328329def squarefree_carpet_variance():330 print("SQUAREFREE ODD STACK the rival selection")331 mu = mobius_table(20000)332 pool = [n for n in range(3, 20001, 2) if mu[n] != 0]333 print(" Cov(15,21) = %s = %.10f Pearson r = %.10f" % (carpet_cov(15, 21), float(carpet_cov(15, 21)), float(carpet_cov(15, 21) / (carpet_var(15) * carpet_var(21)) ** Fraction(1, 2))))334 for count in SQUAREFREE_LS:335 layers = pool[:count]336 var, off, sharing = pair_sum(layers)337 actual = (var + off) / count338 indep = var / count339 print(" L = %4d pairs sharing a factor %6d L*Var = %.7f independent %.7f ratio %.6f" % (count, sharing, actual, indep, actual / indep))340341# DAVENPORT342343def davenport():344 print("DAVENPORT sum mu(n)/n ((nx)) against -sin(2 pi x)/pi")345 top = DAVENPORT_NS[-1]346 mu = mobius_table(top)347 n = np.arange(1, top + 1, dtype=np.int64)348 weight = mu[1:].astype(np.float64) / n.astype(np.float64)349 errors = {}350 for x in DAVENPORT_X:351 a, q = x.numerator, x.denominator352 r = (n * a) % q353 frac = r.astype(np.float64) / q354 saw = np.where(r == 0, 0.0, frac - 0.5)355 target = -math.sin(2 * math.pi * float(x)) / math.pi356 row = []357 for cut in DAVENPORT_NS:358 row.append("%d %.6f" % (cut, float(np.sum(weight[:cut] * saw[:cut]))))359 tail = float(np.sum(weight * frac))360 errors[x] = [abs(float(np.sum(weight[:cut] * saw[:cut])) - target) for cut in DAVENPORT_NS]361 print(" x = %s target %.6f partial sums %s {nx} form %.6f" % (x, target, ", ".join(row), tail))362 for i, cut in enumerate(DAVENPORT_NS):363 print(" cut n <= %6d max error over the five points %.2e" % (cut, max(e[i] for e in errors.values())))364 print(" sum mu(n)/n to %d = %.8f the {nx} form needs this zero, which is PNT-equivalent" % (top, float(np.sum(weight))))365366def main():367 convolution_check()368 selection_closed_forms()369 harmonic_stack()370 prime_carpet_variance()371 squarefree_carpet_variance()372 davenport()373374main()