#!/usr/bin/env python3 """ experiments.py -- numerical verification for "Explicit bounds for odd perfect and quasiperfect numbers" (article_revised.tex). Libraries: numpy (sieves), sympy (primes, factorization, divisors), fractions.Fraction (exact rational identities), mpmath (high precision). Usage: python experiments.py # full run (sieve to 10^7) python experiments.py --sieve-limit 1000000 --brute-limit 200000 # quicker python experiments.py --skip-large # skip the tau = 3906250 row of Experiment E Every check prints PASS or FAIL; the script exits with status 1 if any check fails. Checks whose expected value is tied to the default limits are only asserted when the limits are at least the defaults; otherwise the value is just printed. """ import argparse import itertools import math import sys import time from fractions import Fraction import numpy as np import sympy from mpmath import mp, mpf FAILURES = [] NCHECKS = 0 def check(name, cond): global NCHECKS NCHECKS += 1 ok = bool(cond) print((" PASS " if ok else " FAIL ") + name) if not ok: FAILURES.append(name) return ok def header(title): print() print("=" * 78) print(title) print("=" * 78) # --------------------------------------------------------------------------- # Exact arithmetic helpers # --------------------------------------------------------------------------- def I_from_fact(fact): """Exact abundancy index from a factorization {p: a}.""" r = Fraction(1) for p, a in fact.items(): r *= Fraction(p ** (a + 1) - 1, p ** a * (p - 1)) return r def divisors_from_fact(fact): divs = [1] for p, a in fact.items(): powers = [p ** e for e in range(1, a + 1)] divs = divs + [d * q for q in powers for d in divs] return divs def shifted_exact(n): """Return (I, I+, I-) of n as exact Fractions.""" divs = sympy.divisors(n) I = Fraction(sum(divs), n) Ip = Fraction(1) Im = Fraction(1) for d in divs: if d > 1: Ip += Fraction(1, d + 1) Im += Fraction(1, d - 1) return I, Ip, Im def ratios_exact(n): I, Ip, Im = shifted_exact(n) return I / Ip, Im / I, Im / Ip def n_odd_exponents(fact): return sum(1 for a in fact.values() if a % 2 == 1) def r_fun(x): return math.log(x / (x - 1.0)) / math.log((x + 1.0) / x) # Extremal constants of the shifted-index sections R_STAR = Fraction(3254203392, 3001892677) # R(315) R_E = Fraction(29332068105860657280, 27090788267918857153) # R(2205) S_STAR = Fraction(255, 224) # S(15) S_E = Fraction(18195, 16016) # S(45) T_STAR = Fraction(612, 497) # T(15) T_E = Fraction(836970, 680449) # T(45) # --------------------------------------------------------------------------- # Sieves # --------------------------------------------------------------------------- def odd_sieve(N): """I, I+, I- for odd n <= N (entries at even indices are meaningless).""" I = np.ones(N + 1) Ip = np.ones(N + 1) Im = np.ones(N + 1) for d in range(3, N + 1, 2): step = 2 * d I[d::step] += 1.0 / d Ip[d::step] += 1.0 / (d + 1) Im[d::step] += 1.0 / (d - 1) return I, Ip, Im def full_sieve(N): """I, I+, I- for all n <= N.""" I = np.ones(N + 1) Ip = np.ones(N + 1) Im = np.ones(N + 1) for d in range(2, N + 1): I[d::d] += 1.0 / d Ip[d::d] += 1.0 / (d + 1) Im[d::d] += 1.0 / (d - 1) return I, Ip, Im def spf_sieve(N): """Smallest prime factor of every n <= N.""" spf = np.zeros(N + 1, dtype=np.int64) for i in range(2, int(math.isqrt(N)) + 1): if spf[i] == 0: v = spf[i::i] v[v == 0] = i idx = np.arange(N + 1, dtype=np.int64) mask = spf == 0 spf[mask] = idx[mask] spf[0] = 0 spf[1] = 1 return spf # --------------------------------------------------------------------------- # Experiment A: exhaustive divisor sieve # --------------------------------------------------------------------------- def experiment_A(args): header("Experiment A: exhaustive divisor sieve over odd n <= %d" % args.sieve_limit) N = args.sieve_limit t0 = time.time() I, Ip, Im = odd_sieve(N) print(" sieve time: %.1f s" % (time.time() - t0)) odd = np.arange(3, N + 1, 2) Io, Ipo, Imo = I[odd], Ip[odd], Im[odd] R = Io / Ipo S = Imo / Io T = Imo / Ipo for name, arr, star, where in (("R = I/I+", R, R_STAR, 315), ("S = I-/I", S, S_STAR, 15), ("T = I-/I+", T, T_STAR, 15)): k = int(np.argmax(arr)) print(" max %-10s = %.14f at n = %d (sharp constant %.14f)" % (name, arr[k], odd[k], float(star))) check("max %s over odd n <= N is attained at n = %d" % (name, where), odd[k] == where) check("max %s agrees with the sharp constant to 1e-13" % name, abs(arr[k] - float(star)) < 1e-13) check("no odd n <= N violates the sharp bound for %s" % name, np.all(arr <= float(star) + 1e-13)) check("I+(n) < I(n) < I-(n) for every odd 1 < n <= N", np.all(Ipo < Io) and np.all(Io < Imo)) d1 = Imo - Io d2 = Io - Ipo gap = Imo - Ipo print(" max (I- - I) = %.14f (limit 1 - log 2 = %.14f)" % (d1.max(), 1 - math.log(2))) print(" max (I - I+) = %.14f (limit log 2 - 1/2 = %.14f)" % (d2.max(), math.log(2) - 0.5)) check("max (I- - I) < 1 - log 2 over odd n <= N", d1.max() < 1 - math.log(2)) check("max (I - I+) < log 2 - 1/2 over odd n <= N", d2.max() < math.log(2) - 0.5) k = int(np.argmax(gap)) print(" max (I- - I+) = %.14f at n = %d = %s" % (gap[k], odd[k], sympy.factorint(int(odd[k])))) check("I-(n) - I+(n) < 1/2 for every odd n <= N", gap.max() < 0.5) if N >= 10 ** 7: check("paper value: max (I- - I) = 0.28790584906385", abs(d1.max() - 0.28790584906385) < 1e-12) check("paper value: max (I - I+) = 0.17517080090658", abs(d2.max() - 0.17517080090658) < 1e-12) check("paper value: max (I- - I+) = 0.46307664997043 at n = 3828825", abs(gap[k] - 0.46307664997043) < 1e-12 and odd[k] == 3828825) return I, Ip, Im # --------------------------------------------------------------------------- # Experiments B and C: the reduction lists (exact) # --------------------------------------------------------------------------- def build_lists(): def build(primes, maxexp, euler): out = [] for exps in itertools.product(*[range(m + 1) for m in maxexp]): if euler and sum(1 for e in exps if e % 2 == 1) > 1: continue n = 1 for p, e in zip(primes, exps): n *= p ** e out.append(n) return out L = build([3, 5, 7, 11], [2, 1, 1, 1], False) LE = build([3, 5, 7, 11], [3, 2, 2, 2], True) M = build([3, 5, 7], [1, 1, 1], False) ME = build([3, 5, 7], [2, 2, 2], True) Nl = build([3, 5, 7], [2, 1, 1], False) NE = build([3, 5, 7], [3, 2, 2], True) return L, LE, M, ME, Nl, NE def experiment_BC(args): header("Experiments B and C: the six reduction lists in exact rational arithmetic") L, LE, M, ME, Nl, NE = build_lists() spec = [ ("L (R, odd)", L, 24, 0, R_STAR, 315, Fraction(14, 13)), ("L_E (R, Euler)", LE, 56, 0, R_E, 2205, Fraction(14, 13)), ("M (S, odd)", M, 8, 1, S_STAR, 15, Fraction(9, 8)), ("M_E (S, Euler)", ME, 20, 1, S_E, 45, Fraction(9, 8)), ("N (T, odd)", Nl, 12, 2, T_STAR, 15, Fraction(6, 5)), ("N_E (T, Euler)", NE, 24, 2, T_E, 45, Fraction(6, 5)), ] gaps = [] for name, lst, card, idx, star, where, thr in spec: vals = sorted(((ratios_exact(n)[idx], n) for n in lst), reverse=True) print(" list %-15s cardinality %d; largest values:" % (name, len(lst))) for v, n in vals[:6]: print(" n = %-7d %-28s %.15f" % (n, str(sympy.factorint(n)), float(v))) check("%s has cardinality %d" % (name, card), len(lst) == card) check("%s: maximum equals the stated exact fraction, at n = %d" % (name, where), vals[0][0] == star and vals[0][1] == where) check("%s: the maximum is attained only once" % name, vals[1][0] < vals[0][0]) check("%s: maximum exceeds the threshold %s (exact)" % (name, thr), star > thr) gaps.append(float(vals[0][0] - vals[1][0])) print(" gaps to the runners-up:", ", ".join("%.4e" % g for g in gaps)) expected = [2.4159e-5, 1.1795e-5, 2.4666e-3, 1.0058e-3, 1.4242e-5, 4.1143e-3] check("gaps to the runners-up agree with the paper to 4 significant digits", all(abs(g - e) < 5e-5 * e for g, e in zip(gaps, expected))) # individual values quoted in the text check("R(105) = 5372928/4956455", ratios_exact(105)[0] == Fraction(5372928, 4956455)) check("R(45) = 9568/8837", ratios_exact(45)[0] == Fraction(9568, 8837)) check("R(15) = 384/355", ratios_exact(15)[0] == Fraction(384, 355)) check("S(105) = 385597/339456", ratios_exact(105)[1] == Fraction(385597, 339456)) check("S(75) = 291625/256928", ratios_exact(75)[1] == Fraction(291625, 256928)) check("S(3) = 9/8, S(21) = 721/640, S(35) = 10325/9792", ratios_exact(3)[1] == Fraction(9, 8) and ratios_exact(21)[1] == Fraction(721, 640) and ratios_exact(35)[1] == Fraction(10325, 9792)) check("T(105) = 1348818306/1095376555", ratios_exact(105)[2] == Fraction(1348818306, 1095376555)) check("T(75) = 5762510/4700591", ratios_exact(75)[2] == Fraction(5762510, 4700591)) I15, Ip15, Im15 = shifted_exact(15) check("I+(15) = 71/48 and I-(15) = 51/28", Ip15 == Fraction(71, 48) and Im15 == Fraction(51, 28)) I315, Ip315, _ = shifted_exact(315) check("I(315) = 208/105 and I+(315) = 5574943543/3050815680", I315 == Fraction(208, 105) and Ip315 == Fraction(5574943543, 3050815680)) mp.dps = 40 cp = mp.log(mpf(R_STAR.numerator) / R_STAR.denominator, 2) cpp = mp.log(mpf(S_STAR.numerator) / S_STAR.denominator, 2) t = mp.log(mpf(T_STAR.numerator) / T_STAR.denominator, 2) print(" c' = %s\n c'' = %s\n t = %s" % (mp.nstr(cp, 22), mp.nstr(cpp, 22), mp.nstr(t, 22))) check("t < c' + c'' (T* < R* S*, exact)", T_STAR < R_STAR * S_STAR) check("c' + c'' - t = 0.0031447388723478...", abs(cp + cpp - t - mpf("0.0031447388723478")) < mpf("1e-15")) # remarks on even integers check("R(12) = 12740/10937 > R*", ratios_exact(12)[0] == Fraction(12740, 10937) and Fraction(12740, 10937) > R_STAR) check("S(6) = 27/20 > S*", ratios_exact(6)[1] == Fraction(27, 20)) # --------------------------------------------------------------------------- # Experiment F: certificates (exact integer inequalities) # --------------------------------------------------------------------------- def experiment_F(args): header("Experiment F: threshold certificates as exact integer inequalities") check("13 * 3254203392 = 42304644096 > 42026497478 = 14 * 3001892677", 13 * 3254203392 == 42304644096 and 14 * 3001892677 == 42026497478 and 42304644096 > 42026497478) check("13 * num(R_E) > 14 * den(R_E)", 13 * 29332068105860657280 == 381316885376188544640 and 14 * 27090788267918857153 == 379271035750864000142 and 381316885376188544640 > 379271035750864000142) check("8 * 255 = 2040 > 2016 = 9 * 224", 8 * 255 == 2040 and 9 * 224 == 2016) check("8 * 18195 = 145560 > 144144 = 9 * 16016", 8 * 18195 == 145560 and 9 * 16016 == 144144) check("5 * 612 = 3060 > 2982 = 6 * 497", 5 * 612 == 3060 and 6 * 497 == 2982) check("5 * 836970 = 4184850 > 4082694 = 6 * 680449", 5 * 836970 == 4184850 and 6 * 680449 == 4082694) # --------------------------------------------------------------------------- # Experiment D and the general inequalities, by brute force over odd n # --------------------------------------------------------------------------- def experiment_D_and_general(args, I, Ip, Im, delta): L2 = min(args.sieve_limit, args.brute_limit) header("Experiment D and general inequalities: brute force over odd n <= %d" % L2) t0 = time.time() spf = spf_sieve(L2).tolist() Il, Ipl, Iml = I[:L2 + 1].tolist(), Ip[:L2 + 1].tolist(), Im[:L2 + 1].tolist() c_lin = 8.0 / math.pi ** 2 tol = 1e-12 best = {"R": (0, 0), "S": (0, 0), "T": (0, 0), "RE": (0, 0), "SE": (0, 0), "TE": (0, 0)} viol = {k: 0 for k in ("sieve-vs-multiplicative", "Theorem B > Phi^(1/r(p0))", "Theorem B > (8/pi^2) Phi", "Lemma B <= I / prod(1+1/(p(p+1)))", "identity for I(n rad n)", "I(n rad n) <= I prod(1+1/(p(p+1)))", "I(n rad n) < n/phi(n)", "0 < sum eps(p) < delta", "exp(T - delta) < B < exp(T)", "I- - I+ >= 2 sum 1/(p^2-1)", "sum 1/(p^2-1) > log C", "B > Phi exp(-(I- - I+)/2)", "I < n/phi(n)")} for n in range(3, L2 + 1, 2): m = n B = Phi = C = Imul = Irad = up = radfac = 1.0 T = eps = s2 = logC = 0.0 nodd = 0 omega = 0 p0 = spf[m] while m > 1: p = spf[m] a = 0 while m % p == 0: m //= p a += 1 omega += 1 nodd += a & 1 ip = 1.0 / p B *= 1.0 + ip Phi *= p / (p - 1.0) C *= p * p / (p * p - 1.0) sig = (p ** (a + 1) - 1) // (p - 1) Imul *= sig / float(p ** a) Irad *= ((p ** (a + 2) - 1) // (p - 1)) / float(p ** (a + 1)) radfac *= 1.0 + 1.0 / (p * sig) if a >= 2: up *= 1.0 + 1.0 / (p * (p + 1.0)) T += ip eps += ip - math.log1p(ip) s2 += 1.0 / (p * p - 1.0) logC += math.log1p(1.0 / (p * p - 1.0)) In, Ipn, Imn = Il[n], Ipl[n], Iml[n] if abs(In - Imul) > 1e-9: viol["sieve-vs-multiplicative"] += 1 if omega >= 2 and not B > Phi ** (1.0 / r_fun(p0)): viol["Theorem B > Phi^(1/r(p0))"] += 1 if not B > c_lin * Phi: viol["Theorem B > (8/pi^2) Phi"] += 1 if not B <= Imul / up + tol: viol["Lemma B <= I / prod(1+1/(p(p+1)))"] += 1 if abs(Irad - Imul * radfac) > 1e-12 * Irad: viol["identity for I(n rad n)"] += 1 bound = 1.0 mm = n while mm > 1: p = spf[mm] while mm % p == 0: mm //= p bound *= 1.0 + 1.0 / (p * (p + 1.0)) if not Irad <= Imul * bound + tol: viol["I(n rad n) <= I prod(1+1/(p(p+1)))"] += 1 if not Irad < Phi + tol: viol["I(n rad n) < n/phi(n)"] += 1 if not Imul < Phi: viol["I < n/phi(n)"] += 1 if not 0.0 < eps < delta: viol["0 < sum eps(p) < delta"] += 1 if not (math.exp(T - delta) < B < math.exp(T)): viol["exp(T - delta) < B < exp(T)"] += 1 gap = Imn - Ipn if not gap >= 2.0 * s2 - 1e-11: viol["I- - I+ >= 2 sum 1/(p^2-1)"] += 1 if not s2 > logC: viol["sum 1/(p^2-1) > log C"] += 1 if not B > Phi * math.exp(-gap / 2.0) - tol: viol["B > Phi exp(-(I- - I+)/2)"] += 1 R, S, Tt = In / Ipn, Imn / In, Imn / Ipn if R > best["R"][0]: best["R"] = (R, n) if S > best["S"][0]: best["S"] = (S, n) if Tt > best["T"][0]: best["T"] = (Tt, n) if nodd <= 1: if R > best["RE"][0]: best["RE"] = (R, n) if S > best["SE"][0]: best["SE"] = (S, n) if Tt > best["TE"][0]: best["TE"] = (Tt, n) print(" loop time: %.1f s" % (time.time() - t0)) exp = {"R": (R_STAR, 315), "RE": (R_E, 2205), "S": (S_STAR, 15), "SE": (S_E, 45), "T": (T_STAR, 15), "TE": (T_E, 45)} for key in ("R", "RE", "S", "SE", "T", "TE"): val, n = best[key] star, where = exp[key] cls = "Euler form" if key.endswith("E") else "odd" print(" max %s over %-10s = %.14f at n = %d" % (key[0], cls, val, n)) check("brute force: max %s over %s integers is at n = %d and matches to 1e-13" % (key[0], cls, where), n == where and abs(val - float(star)) < 1e-13) for name, cnt in viol.items(): check("%s: violations = %d" % (name, cnt), cnt == 0) def experiment_even(args): header("Remarks on even integers: full sieve over n <= 200000") N = 200000 I, Ip, Im = full_sieve(N) ev = np.arange(2, N + 1, 2) R = I[ev] / Ip[ev] S = Im[ev] / I[ev] kR, kS = int(np.argmax(R)), int(np.argmax(S)) print(" max R over even n = %.12f at n = %d; max S over even n = %.12f at n = %d" % (R[kR], ev[kR], S[kS], ev[kS])) check("max R(n) over even n <= 2*10^5 is R(12) = 1.164853250434...", ev[kR] == 12 and abs(R[kR] - 12740 / 10937) < 1e-12) check("max S(n) over even n <= 2*10^5 is S(6) = 27/20", ev[kS] == 6 and abs(S[kS] - 1.35) < 1e-12) gap = Im[2:] - Ip[2:] check("I-(n) - I+(n) < 3/2 for every n <= 2*10^5", gap.max() < 1.5) # --------------------------------------------------------------------------- # Experiment E: Eulerian integers with at least ten prime divisors (60 digits) # --------------------------------------------------------------------------- def experiment_E(args): header("Experiment E: members of E_10 at 60 digits, full divisor enumeration") mp.dps = 60 P = [3, 5, 7, 11, 13, 17, 19, 23, 29, 31] def fact(odd_part, squares, power=2): f = dict(odd_part) for p in squares: f[p] = power return f rows = [ ("3 (5..31)^2", fact({3: 1}, P[1:]), (2.871082, 2.699378, 3.153366, 1.063609, 1.098320, 1.168182)), ("(3..29)^2 31", fact({31: 1}, P[:-1]), (3.107207, 2.923302, 3.404504, 1.062910, 1.095680, 1.164609)), ("3^5 (5..31)^2", fact({3: 5}, P[1:]), (3.225537, 3.040003, 3.524570, 1.061031, 1.092708, 1.159397)), ("3 (5..31)^4", fact({3: 1}, P[1:], 4), (2.906542, 2.734750, 3.188915, 1.062818, 1.097151, 1.166072)), ("3 (5..37)^2", fact({3: 1}, P[1:] + [37]), (2.950776, 2.778207, 3.233966, 1.062115, 1.095971, 1.164048)), ("3 (5..41)^2", fact({3: 1}, P[1:] + [37, 41]), (3.024502, 2.851226, 3.308429, 1.060772, 1.093876, 1.160353)), ("101 (103..149)^2", fact({101: 1}, [103, 107, 109, 113, 127, 131, 137, 139, 149]), (1.087571, 1.086865, 1.088290, 1.000650, 1.000661, 1.001312)), ("2205 (101..131)^2", fact({3: 2, 5: 1, 7: 2}, [101, 103, 107, 109, 113, 127, 131]), (2.146937, 1.992203, 2.411320, 1.077670, 1.123144, 1.210379)), ] RE, SE, TE = (mpf(x.numerator) / x.denominator for x in (R_E, S_E, T_E)) print(" %-20s %9s %10s %10s %10s %10s %10s %10s" % ("n", "tau", "I", "I+", "I-", "R", "S", "T")) for name, f, expected in rows: tau = 1 for a in f.values(): tau *= a + 1 if args.skip_large and tau > 10 ** 6: print(" %-20s %9d (skipped: --skip-large)" % (name, tau)) continue t0 = time.time() Ifr = I_from_fact(f) Iv = mpf(Ifr.numerator) / Ifr.denominator Ipv = mpf(1) Imv = mpf(1) one = mpf(1) divs = divisors_from_fact(f) for d in divs: if d > 1: Ipv += one / (d + 1) Imv += one / (d - 1) R, S, T = Iv / Ipv, Imv / Iv, Imv / Ipv got = (Iv, Ipv, Imv, R, S, T) print(" %-20s %9d %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f (%.1f s)" % ((name, tau) + tuple(float(x) for x in got) + (time.time() - t0,))) check("%s: tau = %d, in E_10 (omega >= 10, exactly one odd exponent)" % (name, len(divs)), len(divs) == tau and len(f) >= 10 and n_odd_exponents(f) == 1) check("%s: values agree with Table (6 decimals)" % name, all(abs(float(g) - e) < 1.0e-6 for g, e in zip(got, expected))) check("%s: R < R_E, S < S_E, T < T_E" % name, R < RE and S < SE and T < TE) f = fact({101: 1}, [103, 107, 109, 113, 127, 131, 137, 139, 149]) Ifr = I_from_fact(f) check("I(101 (103..149)^2) = 1.0875712439...", abs(float(Ifr) - 1.0875712439) < 1e-10) # --------------------------------------------------------------------------- # Experiment G: separation, gap identity, witnesses, odd perfect bound # --------------------------------------------------------------------------- def experiment_G(args): header("Experiment G: separation and gap, exact arithmetic for 1 < n <= %d" % args.exact_limit) sep = ident = True bestT = (Fraction(0), 0) for n in range(2, args.exact_limit + 1): divs = sympy.divisors(n) Ip = Fraction(1) Im = Fraction(1) g = Fraction(0) for d in divs: if d > 1: Ip += Fraction(1, d + 1) Im += Fraction(1, d - 1) g += Fraction(2, d * d - 1) if not Im > Ip: sep = False if Im - Ip != g: ident = False if n % 2 == 1 and n_odd_exponents(sympy.factorint(n)) <= 1: Tn = Im / Ip if Tn > bestT[0]: bestT = (Tn, n) check("I-(n) > I+(n) for every 1 < n <= limit (odd and even, exact)", sep) check("I-(n) - I+(n) = sum_{d|n, d>1} 2/(d^2-1) for every n <= limit (exact)", ident) check("max I-/I+ over Euler-form integers <= limit is T_E = 836970/680449 at n = 45", bestT == (T_E, 45)) # witnesses n_k = m_k^2 q for k, lower, omega in ((29, Fraction(7, 15), 10), (31, Fraction(15, 32), 11), (101, Fraction(25, 51), 26)): m = math.lcm(*range(1, k + 1)) while m % 2 == 0: m //= 2 q = int(sympy.nextprime(k)) f = {int(p): 2 * int(a) for p, a in sympy.factorint(m).items()} f[q] = 1 partial = sum(Fraction(2, d * d - 1) for d in range(3, k + 1, 2)) check("witness k = %d: every odd d <= k divides m_k" % k, all(m % d == 0 for d in range(3, k + 1, 2))) check("witness k = %d: partial sum = 1/2 - 1/(k+1) = %s" % (k, lower), partial == Fraction(1, 2) - Fraction(1, k + 1) == lower) check("witness k = %d: n_k in E_10 with omega = %d and one odd exponent" % (k, omega), len(f) == omega and n_odd_exponents(f) == 1 and omega >= 10) # odd perfect bound val = 2 * S_E - 2 / R_E print(" 2 S_E - 2/R_E = %s = %.16f" % (val, float(val))) check("2 S_E - 2/R_E = 43623651349277204129/102662238370512300480", val == Fraction(43623651349277204129, 102662238370512300480)) check("2 S_E - 2/R_E = 0.4249240231041... < 0.4249240232 < 1/2", abs(float(val) - 0.4249240231041) < 1e-12 and val < Fraction(4249240232, 10 ** 10)) check("improvement over 1/2 is 15.02 percent", abs((0.5 - float(val)) / 0.5 * 100 - 15.02) < 0.005) check("2/R_E = 1.847178873998728783... > 1.8471", abs(float(2 / R_E) - 1.847178873998728783) < 1e-14 and 2 / R_E > Fraction(18471, 10000)) check("2 S_E = 18195/8008 = 2.2721028971... < 2.2722", 2 * S_E == Fraction(18195, 8008) and 2 * S_E < Fraction(22722, 10000)) # crude constants check("(log 2 - 1/2) + (1 - log 2) = 1/2", abs((math.log(2) - 0.5) + (1 - math.log(2)) - 0.5) < 1e-15) check("2 - 2^(6/7) = 0.18855267... < log 2 - 1/2", abs(2 - 2 ** (6 / 7) - 0.1885526714721866) < 1e-12 and 2 - 2 ** (6 / 7) < math.log(2) - 0.5) check("2 (S* - 1) = 31/112 < 1 - log 2", 2 * (S_STAR - 1) == Fraction(31, 112) and 31 / 112 < 1 - math.log(2)) # --------------------------------------------------------------------------- # Experiment H: constants of Sections 3-6 and 8 # --------------------------------------------------------------------------- def experiment_H(args): header("Experiment H: constants for B(n), n/phi(n), I(n rad n) (mpmath, 40 digits)") mp.dps = 40 pi, z3 = mp.pi, mp.zeta(3) m = mp.log(mpf(3) / 2) / mp.log(mpf(4) / 3) print(" m = r(3) = %s, 2^(1/m) = %s" % (mp.nstr(m, 20), mp.nstr(mpf(2) ** (1 / m), 20))) check("2^(1/m) = 1.635250087...", abs(mpf(2) ** (1 / m) - mpf("1.635250087185844460878")) < mpf("1e-20")) thr = (pi ** 2 / 8) ** (m / (m - 1)) check("(pi^2/8)^(m/(m-1)) = 2.060569542...", abs(thr - mpf("2.060569542")) < mpf("1e-9")) check("r is decreasing on a grid of x in (1, 200]", all(r_fun(1.0 + 0.01 * (i + 1)) < r_fun(1.0 + 0.01 * i) for i in range(1, 19900))) check("(1+1/p)^r(p) = p/(p-1) for primes p < 1000", all(abs((1 + 1.0 / p) ** r_fun(p) - p / (p - 1.0)) < 1e-12 for p in sympy.primerange(2, 1000))) check("Z(2) = pi^2/6, Z(3) = pi^2/8, Z(5) = pi^2/9, Z(7) = 8 pi^2/75", abs(mp.zeta(2) - pi ** 2 / 6) < mpf("1e-35") and abs(mp.zeta(2) * mpf(3) / 4 - pi ** 2 / 8) < mpf("1e-35") and abs(mp.zeta(2) * mpf(3) / 4 * mpf(8) / 9 - pi ** 2 / 9) < mpf("1e-35") and abs(mp.zeta(2) * mpf(3) / 4 * mpf(8) / 9 * mpf(24) / 25 - 8 * pi ** 2 / 75) < mpf("1e-35")) P3 = 8 / (7 * z3) print(" P3 = 8/(7 zeta(3)) = %s" % mp.nstr(P3, 20)) eps35 = mpf(10) ** -35 check("P3 = 0.950751282949...", abs(P3 - mpf("0.950751282949379")) < mpf("1e-14")) check("zeta(3) < 1.2020569032", z3 < mpf("1.2020569032")) P3low = 8 / (7 * mpf("1.2020569032")) check("zeta(3) < 1.2020569032 gives P3 > 0.9507512829", P3low > mpf("0.9507512829")) check("(2+10^-35)/0.9507512829 < 2.1035995807 < 2.10360", (2 + eps35) / mpf("0.9507512829") < mpf("2.1035995807") < mpf("2.10360")) check("31 (2+10^-35)/(30 * 0.9507512829) < 2.1737195667 < 2.17372", 31 * (2 + eps35) / (30 * mpf("0.9507512829")) < mpf("2.1737195667") < mpf("2.17372")) check("(30/31) P3 = 0.920081886...", abs(mpf(30) / 31 * P3 - mpf("0.920081886")) < mpf("1e-9")) check("(2+10^-35)/P3 = 2.1035995... and 31(2+10^-35)/(30 P3) = 2.1737195...", abs((2 + eps35) / P3 - mpf("2.1035995")) < mpf("1e-7") and abs(31 * (2 + eps35) / (30 * P3) - mpf("2.1737195")) < mpf("1e-7")) check("1 - 1/(q^2+q+1) = (1-q^-2)/(1-q^-3) and is >= 30/31 for primes 5 <= q < 1000", all(Fraction(q * q - 1, q * q) / Fraction(q ** 3 - 1, q ** 3) == 1 - Fraction(1, q * q + q + 1) >= Fraction(30, 31) for q in sympy.primerange(5, 1000))) check("Phi(105) = (3/2)(5/4)(7/6) = 2.1875 > 2.17372", Fraction(3, 2) * Fraction(5, 4) * Fraction(7, 6) == Fraction(35, 16) > Fraction(217372, 100000)) for name, val, dec in (("49/(3 pi^2)", 49 / (3 * pi ** 2), "1.654912666"), ("18/pi^2", 18 / pi ** 2, "1.823781305"), ("6 pi^2/49", 6 * pi ** 2 / 49, "1.208522987"), ("pi^2/9", pi ** 2 / 9, "1.096622711"), ("24/13", mpf(24) / 13, "1.846153846"), ("log(49/(3 pi^2))", mp.log(49 / (3 * pi ** 2)), "0.503748237"), ("log(18/pi^2)", mp.log(18 / pi ** 2), "0.600911986")): check("%s = %s..." % (name, dec), 0 <= val - mpf(dec) < mpf("1e-9")) check("49/(3 pi^2) > 2^(1/m)", 49 / (3 * pi ** 2) > mpf(2) ** (1 / m)) check("(1+1/p)(1+1/(p(p+1))) = 1+1/p+1/p^2 for primes p < 1000", all((1 + Fraction(1, p)) * (1 + Fraction(1, p * (p + 1))) == 1 + Fraction(1, p) + Fraction(1, p * p) for p in sympy.primerange(2, 1000))) # kappa kappa = 84 * z3 / pi ** 4 print(" kappa = 84 zeta(3)/pi^4 = %s" % mp.nstr(kappa, 20)) check("kappa = 1.036584766...", 0 <= kappa - mpf("1.036584766") < mpf("1e-9")) part = mpf(1) for p in sympy.primerange(3, 10 ** 6): part *= 1 + mpf(1) / (p * (p * p + p + 1)) check("partial product of (1+1/(p(p^2+p+1))) over odd primes < 10^6 is below kappa and within 1e-12", part < kappa and kappa - part < mpf("1e-12")) check("1+1/(p(p^2+p+1)) = I(p^3)/I(p^2) = (1-p^-4)/(1-p^-3) for primes p < 1000", all(1 + Fraction(1, p * (p * p + p + 1)) == Fraction(p ** 4 - 1, p ** 3 * (p - 1)) / Fraction(p ** 3 - 1, p ** 2 * (p - 1)) == (1 - Fraction(1, p ** 4)) / (1 - Fraction(1, p ** 3)) for p in sympy.primerange(2, 1000))) check("kappa (2+10^-35) < 2.07317", kappa * (2 + eps35) < mpf("2.07317")) check("(31/30) kappa * 2 = 2.1422751... < 2.14228", abs(mpf(31) / 30 * kappa * 2 - mpf("2.1422751")) < mpf("1e-7") and mpf(31) / 30 * kappa * 2 < mpf("2.14228")) check("sigma(p^(a+1)) = p sigma(p^a) + 1 for p < 50, a <= 8", all((p ** (a + 2) - 1) // (p - 1) == p * ((p ** (a + 1) - 1) // (p - 1)) + 1 for p in sympy.primerange(2, 50) for a in range(1, 9))) # exact check of the rad identity ok = True for n in range(2, 3001): f = sympy.factorint(n) f2 = {p: a + 1 for p, a in f.items()} prod = Fraction(1) for p, a in f.items(): prod *= 1 + Fraction(1, p * ((p ** (a + 1) - 1) // (p - 1))) if I_from_fact(f2) != I_from_fact(f) * prod: ok = False check("I(n rad n) = I(n) prod (1 + 1/(p sigma(p^a))) exactly for 2 <= n <= 3000", ok) # log(1+x+x^2) > x check("log(1+x+x^2) > x on a grid of x in (0,1)", all(math.log(1 + x + x * x) > x for x in (i / 10000.0 for i in range(1, 10000)))) return kappa # --------------------------------------------------------------------------- # Experiment I: the least prime factor and the number of prime factors # --------------------------------------------------------------------------- def experiment_I(args): header("Experiment I: least prime factor versus number of prime factors (Rosser-Schoenfeld)") mp.dps = 30 g = float(mp.euler) # the Rosser-Schoenfeld product bounds on a sample of x prod = 1.0 ok = True primes = list(sympy.primerange(2, 200001)) samples = set([2, 3, 5, 10, 100, 286, 1000, 10 ** 4, 10 ** 5, 2 * 10 ** 5]) j = 0 for x in range(2, 200001): while j < len(primes) and primes[j] <= x: prod *= primes[j] / (primes[j] - 1.0) j += 1 if x in samples or x % 997 == 0: lx = math.log(x) lo = math.exp(g) * lx * (1 - 1 / (2 * lx * lx)) hi = math.exp(g) * lx * (1 + 1 / (lx * lx)) if not lo < prod < hi: ok = False check("e^gamma log x (1 - 1/(2 log^2 x)) < prod_{p<=x} p/(p-1) < e^gamma log x (1 + 1/log^2 x) on samples x <= 2*10^5", ok) check("pi(x) > x/log x for 17 <= x <= 10^5 and pi(x) < 1.25506 x/log x for 2 <= x <= 10^5 (sampled)", all(int(sympy.primepi(x)) > x / math.log(x) for x in range(17, 100001, 37)) and all(int(sympy.primepi(x)) < 1.25506 * x / math.log(x) for x in range(2, 100001, 37))) def X_of(p1): L = math.log(p1 - 1) return (p1 - 1) ** 2 * math.exp(-2.0 / L) def k_bound(p1): return int(sympy.primepi(int(math.floor(X_of(p1))))) - int(sympy.primepi(p1 - 1)) + 1 vals = [k_bound(p) for p in (3, 5, 7, 11, 13, 17, 19, 23)] print(" bound (pi(X) - pi(p1-1) + 1) for p1 = 3,5,7,11,13,17,19,23:", vals) check("these equal 0, 1, 3, 10, 14, 25, 31, 47", vals == [0, 1, 3, 10, 14, 25, 31, 47]) X11 = X_of(11) check("X(11) = 100 exp(-2/log 10) = 41.954... >= 17", abs(X11 - 41.954) < 1e-3 and X11 >= 17) X = X_of(101) expl = X / math.log(X) - 1.25506 * 100 / math.log(100) + 1 print(" p1 = 101: X = %.4f, explicit bound = %.4f, pi(6477) = %d, pi(100) = %d" % (X, expl, int(sympy.primepi(6477)), int(sympy.primepi(100)))) check("p1 = 101: X = 6477.21..., explicit bound 711.80... (k >= 712)", abs(X - 6477.21) < 0.01 and abs(expl - 711.80) < 0.01 and math.ceil(expl) == 712) check("p1 = 101: pi(6477) = 840, pi(100) = 25, so k >= 816", int(sympy.primepi(6477)) == 840 and int(sympy.primepi(100)) == 25 and k_bound(101) == 816) # direct verification: least k with prod of the first k primes >= p1 exceeding 2 ok_y = ok_k = ok_expl = True for p1 in sympy.primerange(3, 200): prodq = mpf(1) q = p1 k = 0 while True: prodq *= mpf(q) / (q - 1) k += 1 if prodq > 2: break q = int(sympy.nextprime(q)) y = q if not y > X_of(p1): ok_y = False if not k >= k_bound(p1): ok_k = False if p1 >= 11: Xp = X_of(p1) if not k > Xp / math.log(Xp) - 1.25506 * (p1 - 1) / math.log(p1 - 1) + 1: ok_expl = False check("for every prime 3 <= p1 < 200: the least admissible k has k-th prime y > X", ok_y) check("for every prime 3 <= p1 < 200: the least admissible k satisfies k >= pi(X) - pi(p1-1) + 1", ok_k) check("for every prime 11 <= p1 < 200: the least admissible k exceeds the explicit bound", ok_expl) # --------------------------------------------------------------------------- # Experiment J: Mertens constant, delta, and the gcd(15, N) bounds # --------------------------------------------------------------------------- def experiment_J(args): header("Experiment J: sum of eps(p), delta, and the bounds according to gcd(15, N)") mp.dps = 40 pi = mp.pi def eps(p): return mpf(1) / p - mp.log(1 + mpf(1) / p) M, gam = mp.mertens, mp.euler D = mp.log(mp.zeta(2)) - gam + M delta = D - eps(2) print(" M = %s\n gamma = %s\n sum_p eps(p) = %s\n delta = %s" % (mp.nstr(M, 20), mp.nstr(gam, 20), mp.nstr(D, 20), mp.nstr(delta, 20))) check("M = 0.2614972128...", 0 <= M - mpf("0.2614972128") < mpf("1e-10")) check("log zeta(2) - gamma + M = 0.1819818504...", 0 <= D - mpf("0.1819818504") < mpf("1e-10")) check("delta = 0.0874469585...", 0 <= delta - mpf("0.0874469585") < mpf("1e-10")) check("eps(2) = 1/2 - log(3/2)", abs(eps(2) - (mpf(1) / 2 - mp.log(mpf(3) / 2))) < mpf("1e-35")) part = mpf(0) partM = mpf(0) limit = 2 * 10 ** 6 for p in sympy.primerange(2, limit): part += eps(p) partM += mp.log(1 - mpf(1) / p) + mpf(1) / p print(" partial sum of eps(p) over p < %d: %s" % (limit, mp.nstr(part, 15))) check("partial sum of eps(p) over primes < 2*10^6 is below the closed form and within 1e-7", part < D and D - part < mpf("1e-7")) check("gamma + partial sum of (log(1-1/p) + 1/p) over primes < 2*10^6 is within 1e-7 of M", abs(gam + partM - M) < mpf("1e-7")) check("log(1+1/p) = log(1-1/p^2) - log(1-1/p) for primes p < 1000", all(abs(mp.log(1 + mpf(1) / p) - (mp.log(1 - mpf(1) / p ** 2) - mp.log(1 - mpf(1) / p))) < mpf("1e-35") for p in sympy.primerange(2, 1000))) # (a) lower bounds valid for odd perfect and quasiperfect numbers lam = {1: 75 / (8 * pi ** 2), 3: 25 / (3 * pi ** 2), 5: 9 / pi ** 2, 15: 49 / (6 * pi ** 2)} dec = {1: "1.899772193", 3: "1.688686394", 5: "1.823781305", 15: "1.654912666"} for g in (1, 3, 5, 15): check("2 lambda_%d = %s..." % (g, dec[g]), 0 <= 2 * lam[g] - mpf(dec[g]) < mpf("1e-9")) check("1/lambda_g = 8 pi^2/75, 3 pi^2/25, pi^2/9, 6 pi^2/49 are the stated bounds for C(N)", abs(1 / lam[1] - mp.zeta(2) * mpf(3) / 4 * mpf(8) / 9 * mpf(24) / 25) < mpf("1e-35") and abs(1 / lam[3] - pi ** 2 / 8 * (1 - mpf(1) / 25)) < mpf("1e-35") and abs(1 / lam[5] - pi ** 2 / 9) < mpf("1e-35") and abs(1 / lam[15] - pi ** 2 / 8 * (1 - mpf(1) / 49)) < mpf("1e-35")) # quoted bounds for T(N) (Cohen, Hagis, Suryanarayana) Tlo = {1: mpf("0.667450"), 3: mpf("0.603831"), 5: mpf("0.647387"), 15: mpf("0.596063")} Thi = {1: mpf("0.693148"), 3: mpf("0.657304"), 5: mpf("0.678036"), 15: mpf("0.673770")} # (b) upper bounds up = {3: mpf(4) / 3 * mp.exp(Thi[3] - mpf(1) / 3), 5: mpf(6) / 5 * mp.exp(Thi[5] - mpf(1) / 5), 15: mpf(8) / 5 * mp.exp(Thi[15] - mpf(8) / 15)} updec = {3: "1.84348", 5: "1.93549", 15: "1.84125"} for g in (3, 5, 15): print(" g = %-2d upper bound for B(N): %s" % (g, mp.nstr(up[g], 12))) check("g = %d: upper bound < %s (and within 1e-5)" % (g, updec[g]), up[g] < mpf(updec[g]) and mpf(updec[g]) - up[g] < mpf("1e-5")) check("upper bounds equal exp(T_hi - sum_{p|g} eps(p))", abs(up[3] - mp.exp(Thi[3] - eps(3))) < mpf("1e-30") and abs(up[5] - mp.exp(Thi[5] - eps(5))) < mpf("1e-30") and abs(up[15] - mp.exp(Thi[15] - eps(3) - eps(5))) < mpf("1e-30")) # (c) lower bounds for odd perfect numbers dg = {1: delta - eps(3) - eps(5), 3: delta - eps(5), 5: delta - eps(3), 15: delta - eps(7)} dgdec = {1: "0.0241172", 3: "0.0697685", 5: "0.0417956", 15: "0.0781212"} lodec = {1: "1.90281", 3: "1.70584", 5: "1.83233", 15: "1.67856"} lo = {} for g in (1, 3, 5, 15): lo[g] = mp.exp(Tlo[g] - dg[g]) print(" g = %-2d delta_g = %s, lower bound for B(N): %s" % (g, mp.nstr(dg[g], 12), mp.nstr(lo[g], 12))) check("g = %d: delta_g = %s..." % (g, dgdec[g]), 0 <= dg[g] - mpf(dgdec[g]) < mpf("1e-7")) check("g = %d: lower bound > %s (and within 1e-5)" % (g, lodec[g]), lo[g] > mpf(lodec[g]) and lo[g] - mpf(lodec[g]) < mpf("1e-5")) check("g = %d: lower bound of (c) exceeds 2 lambda_g of (a)" % g, lo[g] > 2 * lam[g]) # bounds for C(N) A(N) = 2/B(N) ca = {1: ("1", "1.05108", None, "1.90281"), 3: ("1.08490", "1.17245", "1.84348", "1.70584"), 5: ("1.03333", "1.09151", "1.93549", "1.83233"), 15: ("1.08621", "1.19150", "1.84125", "1.67856")} for g, (clo, chi, bup, blo) in ca.items(): ok = 2 / mpf(blo) < mpf(chi) if bup is not None: ok = ok and 2 / mpf(bup) > mpf(clo) check("g = %d: %s < C(N)A(N) < %s follows from the bounds for B(N)" % (g, clo, chi), ok) check("13/12 < 6 pi^2/49 and 1 < pi^2/9", mpf(13) / 12 < 6 * pi ** 2 / 49 and 1 < pi ** 2 / 9) return float(delta) # --------------------------------------------------------------------------- # Experiment L: the functions sigma^+-, phi^+, psi^+ and the series for I^+- # --------------------------------------------------------------------------- def experiment_L(args): header("Experiment L: sigma^+, sigma^-, psi^+, phi^+ (Sandor-Atanassov) and the series for I^+, I^-") mp.dps = 40 pi, z3 = mp.pi, mp.zeta(3) rho = pi ** 2 / (7 * z3) theta = 105 * z3 / pi ** 4 print(" rho = pi^2/(7 zeta(3)) = %s\n theta = 105 zeta(3)/pi^4 = %s" % (mp.nstr(rho, 20), mp.nstr(theta, 20))) check("rho = 1.172942380...", 0 <= rho - mpf("1.172942380") < mpf("1e-9")) check("theta = 1.295730957...", 0 <= theta - mpf("1.295730957") < mpf("1e-9")) check("rho = P3 * pi^2/8 and theta = (2/3) zeta(2) zeta(3)/zeta(6)", abs(rho - 8 / (7 * z3) * pi ** 2 / 8) < mpf("1e-35") and abs(theta - mpf(2) / 3 * mp.zeta(2) * z3 / mp.zeta(6)) < mpf("1e-35")) check("zeta(6) = pi^6/945", abs(mp.zeta(6) - pi ** 6 / 945) < mpf("1e-35")) pr = mpf(1) pt = mpf(1) for p in sympy.primerange(3, 10 ** 6): pr *= 1 + mpf(1) / (p * (p + 1)) pt *= 1 + mpf(1) / (p * (p - 1)) check("partial product of 1+1/(p(p+1)) over odd primes < 10^6 is below rho and within 1e-6", pr < rho and rho - pr < mpf("1e-6")) check("partial product of 1+1/(p(p-1)) over odd primes < 10^6 is below theta and within 1e-6", pt < theta and theta - pt < mpf("1e-6")) check("1+1/(p(p+1)) = (1-p^-3)/(1-p^-2) and 1+1/(p(p-1)) = (1-p^-6)/((1-p^-3)(1-p^-2)) for primes p < 1000", all(1 + Fraction(1, p * (p + 1)) == (1 - Fraction(1, p ** 3)) / (1 - Fraction(1, p ** 2)) and 1 + Fraction(1, p * (p - 1)) == (1 - Fraction(1, p ** 6)) / ((1 - Fraction(1, p ** 3)) * (1 - Fraction(1, p ** 2))) for p in sympy.primerange(2, 1000))) for name, val, dec in (("5/(6 rho)", 5 / (6 * rho), "0.710463998"), ("7 rho/6", 7 * rho / 6, "1.368432777"), ("5 theta/4", 5 * theta / 4, "1.619663697"), ("14 zeta(3)/pi^2 = 2/rho", 2 / rho, "1.705113595"), ("5/(3 rho)", 5 / (3 * rho), "1.420927996")): check("%s = %s..." % (name, dec), 0 <= val - mpf(dec) < mpf("1e-9")) check("14 zeta(3)/pi^2 = 2/rho", abs(14 * z3 / pi ** 2 - 2 / rho) < mpf("1e-35")) def local(n): fact = sympy.factorint(n) sig = sigp = sigm = psi = psip = phi = phip = 1 for p, a in fact.items(): s_ = (p ** (a + 1) - 1) // (p - 1) ps = p ** (a - 1) * (p + 1) ph = p ** (a - 1) * (p - 1) sig *= s_ sigp *= s_ + 1 sigm *= s_ - 1 psi *= ps psip *= ps + 1 phi *= ph phip *= ph + 1 return fact, sig, sigp, sigm, psi, psip, phi, phip # exact identities ok_minus = ok_ratio = True for n in range(2, 3001): fact, sig, sigp, sigm, psi, psip, phi, phip = local(n) rad = 1 for p in fact: rad *= p if sigm != rad * int(sympy.divisor_sigma(n // rad)): ok_minus = False r1 = r2 = r3 = r4 = Fraction(1) for p, a in fact.items(): r1 *= 1 + Fraction(1, (p ** (a + 1) - 1) // (p - 1)) r2 *= 1 + Fraction(1, p * ((p ** a - 1) // (p - 1))) r3 *= 1 + Fraction(1, p ** (a - 1) * (p + 1)) r4 *= 1 + Fraction(1, p ** (a - 1) * (p - 1)) if not (Fraction(sigp, sig) == r1 and Fraction(sig, sigm) == r2 and Fraction(psip, psi) == r3 and Fraction(phip, phi) == r4): ok_ratio = False check("sigma^-(n) = rad(n) sigma(n/rad(n)) for 2 <= n <= 3000 (exact)", ok_minus) check("the four product formulas for sigma^+/sigma, sigma/sigma^-, psi^+/psi, phi^+/phi for 2 <= n <= 3000 (exact)", ok_ratio) # inequalities of the theorem rho_f, theta_f = float(rho), float(theta) limit = min(args.brute_limit, 2 * 10 ** 5) va = vb = 0 na = nb = 0 for n in range(3, limit + 1, 2): fact, sig, sigp, sigm, psi, psip, phi, phip = local(n) ones = [p for p, a in fact.items() if a == 1] if not ones: na += 1 if not (sig / rho_f < sigm < sig < sigp < rho_f * sig and psi < psip < rho_f * psi and phi < phip < theta_f * phi): va += 1 elif len(ones) == 1 and ones[0] >= 5: nb += 1 if not (5 / (6 * rho_f) * sig < sigm < sig < sigp < 7 * rho_f / 6 * sig and psip < 7 * rho_f / 6 * psi and phip < 5 * theta_f / 4 * phi): vb += 1 check("part (a): %d odd n <= %d with all exponents >= 2, violations = %d" % (na, limit, va), va == 0) check("part (b) hypotheses (one exponent 1 at a prime >= 5): %d odd n <= %d, violations = %d" % (nb, limit, vb), vb == 0) # the series for I^- and I^+ ok_series = True for n in list(range(2, 400)) + [945, 2205, 3465, 45045]: divs = [d for d in sympy.divisors(n) if d > 1] Im = 1 + sum(mpf(1) / (d - 1) for d in divs) Ip = 1 + sum(mpf(1) / (d + 1) for d in divs) sm = sp = se = mpf(0) for k in range(1, 140): t = sum(mpf(d) ** (-k) for d in divs) sm += t sp += (-1) ** (k + 1) * t if k % 2 == 0: se += t if abs(1 + sm - Im) > mpf("1e-30") or abs(1 + sp - Ip) > mpf("1e-30") or abs(2 * se - (Im - Ip)) > mpf("1e-30"): ok_series = False check("I^- = 1 + sum_k (sigma_{-k} - 1), I^+ = 1 + sum_k (-1)^(k+1) (sigma_{-k} - 1), I^- - I^+ = 2 sum_j (sigma_{-2j} - 1) on a sample of n", ok_series) check("2 S_E - 2 = 2179/8008 = 0.2721... and 0.4249240232/2 = 0.2124620116", 2 * S_E - 2 == Fraction(2179, 8008) and abs(float(2 * S_E - 2) - 0.2721) < 1e-4 and Fraction(4249240232, 2 * 10 ** 10) == Fraction(2124620116, 10 ** 10)) # --------------------------------------------------------------------------- def main(): ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) ap.add_argument("--sieve-limit", type=int, default=10 ** 7, help="limit of the divisor sieve (Experiment A)") ap.add_argument("--brute-limit", type=int, default=2 * 10 ** 6, help="limit of the brute-force loop (Experiment D)") ap.add_argument("--exact-limit", type=int, default=2 * 10 ** 4, help="limit of the exact loop (Experiment G)") ap.add_argument("--skip-large", action="store_true", help="skip the tau = 3906250 row of Experiment E") args = ap.parse_args() if args.sieve_limit < 3000 or args.brute_limit < 3000 or args.exact_limit < 50: ap.error("limits too small: need sieve-limit, brute-limit >= 3000 and exact-limit >= 50") t0 = time.time() print("numpy %s, sympy %s, mpmath %s" % (np.__version__, sympy.__version__, __import__("mpmath").__version__)) delta = experiment_J(args) experiment_H(args) experiment_I(args) experiment_L(args) experiment_F(args) experiment_BC(args) experiment_G(args) experiment_E(args) I, Ip, Im = experiment_A(args) experiment_D_and_general(args, I, Ip, Im, delta) experiment_even(args) header("Summary") print(" checks run: %d, failed: %d, total time: %.1f s" % (NCHECKS, len(FAILURES), time.time() - t0)) for name in FAILURES: print(" FAILED: " + name) sys.exit(1 if FAILURES else 0) if __name__ == "__main__": main()