"""Verification script for "Robin's Inequality and Superabundant Counterexamples" (F. Vega). Written with the help of Claude Opus 5.5 (Anthropic). Part 1 checks the constants of the constants lemma and the Mertens-product lemma. Part 2 checks the constants of the theorems with the exponent 1.016. Part 3 checks the constants of the section on the largest prime factor of a counterexample (Axler's Proposition 6.1 and the theta bound). Part 4 enumerates all superabundant numbers up to 10**NEXP (rigorously, as the records of sigma(n)/n among Hardy-Ramanujan integers) and tests on them the superabundant corollary (a)/(b) and the inequality e^0.0562 > 1 + x log^2 p_k, where x = loglog n - log p_k. Part 5 tests on the same superabundant numbers the hypothesis of the sufficient condition for the Riemann Hypothesis, x >= c / log^alpha(p_k), for alpha = 1, ..., 5. Self-contained: needs only Python 3.12 and mpmath (no package install). Run: python car/verify_robin_d.py [N_EXPONENT] (default bound 10**40) """ import sys from mpmath import mp, mpf, euler, exp, log, pi, zeta mp.dps = 60 results = [] def check(label, ok, detail=""): results.append(ok) print(f"[{'PASS' if ok else 'FAIL'}] {label}" + (f" ({detail})" if detail else "")) def primes_upto(m): sieve = bytearray([1]) * (m + 1) sieve[0:2] = b"\x00\x00" for i in range(2, int(m ** 0.5) + 1): if sieve[i]: sieve[i * i :: i] = bytearray(len(sieve[i * i :: i])) return [i for i in range(m + 1) if sieve[i]] def theta(x, primes): return sum(log(mpf(p)) for p in primes if p <= x) def mertens(x, primes): prod = mpf(1) for p in primes: if p > x: break prod *= mpf(p) / (p - 1) return prod # --------------------------------------------------------------------------- # Lemma: numerical constants # --------------------------------------------------------------------------- print("=== Lemma: numerical constants ===") eg = exp(euler) z2 = pi ** 2 / 6 assert abs(z2 - zeta(2)) < mpf(10) ** -50 d = eg / z2 print(f"gamma = {mp.nstr(euler, 15)}") print(f"e^gamma = {mp.nstr(eg, 15)}") print(f"zeta(2) = {mp.nstr(z2, 15)}") print(f"d = {mp.nstr(d, 15)}") check("(i) 1.0827 < d < 1.0828", mpf("1.0827") < d < mpf("1.0828")) c2 = eg * (d - 1) check("(ii) e^gamma (d-1) = e^{2gamma}/zeta(2) - e^gamma > 0.1474", c2 > mpf("0.1474") and abs(c2 - (exp(2 * euler) / z2 - eg)) < mpf(10) ** -50, mp.nstr(c2, 12)) l64sq = log(64) ** 2 check("(iii) log^2(64) >= 17.296", l64sq >= mpf("17.296"), mp.nstr(l64sq, 12)) check("(iii) 0.1474*17.296 > 2.549", mpf("0.1474") * mpf("17.296") > mpf("2.549")) check("(iii) e^gamma (d-1) log^2(64) > 2.51", c2 * l64sq > mpf("2.51"), mp.nstr(c2 * l64sq, 12)) P = primes_upto(2000) th73 = theta(73, P) th71 = theta(71, P) check("(iv) theta(73) > 64", th73 > 64, f"theta(73) = {mp.nstr(th73, 12)}") check("(iv) theta(71) < 64 (so 71 is the largest prime allowed in Case II)", th71 < 64, f"theta(71) = {mp.nstr(th71, 12)}") L285 = log(285) fac285 = 1 / (1 - 1 / (2 * L285 ** 2)) check("(v) (1 - 1/(2 log^2 285))^{-1} < 1.016", fac285 < mpf("1.016"), mp.nstr(fac285, 12)) check("(v) 1.016 < d", mpf("1.016") < d) B285 = eg * L285 * fac285 check("(vi) e^gamma log 285 (1 - 1/(2 log^2 285))^{-1} < 10.23", B285 < mpf("10.23"), mp.nstr(B285, 12)) v313 = eg * log(313) check("(vii) e^gamma log 313 > 10.23", v313 > mpf("10.23"), mp.nstr(v313, 12)) check("(vii) e^313 < 10^136", mpf(313) < 136 * log(10), f"313/ln10 = {mp.nstr(313 / log(10), 10)}") # --------------------------------------------------------------------------- # Mertens-product lemma: true values versus the bounds used # --------------------------------------------------------------------------- print("\n=== Lemma: Mertens products ===") cum = {} prod = mpf(1) for p in P: prod *= mpf(p) / (p - 1) cum[p] = prod m285 = cum[283] m71 = cum[71] check("(a) prod_{p<=285} p/(p-1) < 10.23 (true value)", m285 < mpf("10.23"), mp.nstr(m285, 12)) check(" prod_{p<=71} p/(p-1) < 10.23 (Case II of Theorem 1)", m71 < mpf("10.23"), mp.nstr(m71, 12)) big = [p for p in P if p >= 285] ok_b = all(cum[p] < mpf("1.016") * eg * log(p) for p in big) check("(b) prod_{p<=x} p/(p-1) < 1.016 e^gamma log x for primes 285 <= x <= 2000", ok_b) worst = max((cum[p] / (eg * log(p)), p) for p in big) print(f" max ratio prod/(e^gamma log x) on primes in [285, 2000]: {mp.nstr(worst[0], 10)} at x = {worst[1]}") bad = [p for p in P if cum[p] >= d * eg * log(p)] print(f" primes x <= 2000 with prod_{{p<=x}} p/(p-1) >= d e^gamma log x: largest is {max(bad) if bad else None}") print("\n=== Constants for the exponent 1.016 ===") PB = primes_upto(20000) th14000 = theta(14000, PB) check("theta(14000) > 13000", th14000 > 13000, mp.nstr(th14000, 12)) v = eg * mpf("0.016") * log(13000) ** 2 check("e^gamma * 0.016 * log^2(13000) > 2.51", v > mpf("2.51"), mp.nstr(v, 12)) v = mpf("1.016") * eg * log(14000) check("1.016 * e^gamma * log(14000) < 17.3", v < mpf("17.3"), mp.nstr(v, 12)) m14000 = mertens(14000, PB) check("prod_{p<=14000} p/(p-1) < 17.3 (true value)", m14000 < mpf("17.3"), mp.nstr(m14000, 12)) v = eg * log(16600) check("e^gamma * log(16600) > 17.3", v > mpf("17.3"), mp.nstr(v, 12)) check("e^16600 < 10^7210", mpf(16600) < 7210 * log(10), f"16600/ln10 = {mp.nstr(16600 / log(10), 10)}") # --------------------------------------------------------------------------- # Constants of the section on the largest prime factor of a counterexample # --------------------------------------------------------------------------- print("\n=== Constants for the largest prime factor of a counterexample ===") LMP = mpf("13.11485") * log(10) + log(log(10)) # loglog of 10^(10^13.11485) check("loglog(10^(10^13.11485)) > 31.032", LMP > mpf("31.032"), mp.nstr(LMP, 12)) check("31.032/1.016 > 30.54", mpf("31.032") / mpf("1.016") > mpf("30.54"), mp.nstr(mpf("31.032") / mpf("1.016"), 12)) check("e^30.54 > 7232121212 > 46909038", exp(mpf("30.54")) > 7232121212 > 46909038, mp.nstr(exp(mpf("30.54")), 12)) # Axler, Integers 18 (2018), Proposition 6.1 (Mertens product, x >= 46909038) P0 = mpf("30.54") eps0 = 1 / (20 * P0 ** 3) + mpf(3) / (16 * P0 ** 4) lam = 1 / (1 - eps0) check("(1 - eps(30.54))^{-1} < 1.000002, eps(P) = 1/(20 P^3) + 3/(16 P^4)", lam < mpf("1.000002"), mp.nstr(lam, 15)) c = (mpf("0.05") + mpf(3) / (16 * P0)) / (1 - eps0) check("(1/20 + 3/(16*30.54))/(1 - eps(30.54)) < 0.0562", c < mpf("0.0562"), mp.nstr(c, 12)) check("0.0562/30.54^2 < 0.0000603", mpf("0.0562") / P0 ** 2 < mpf("0.0000603"), mp.nstr(mpf("0.0562") / P0 ** 2, 12)) check("31.032/1.000002 > 31.0319", mpf("31.032") / mpf("1.000002") > mpf("31.0319"), mp.nstr(mpf("31.032") / mpf("1.000002"), 12)) check("0.0562 * 1.000002^2 < 0.056201", mpf("0.0562") * mpf("1.000002") ** 2 < mpf("0.056201"), mp.nstr(mpf("0.0562") * mpf("1.000002") ** 2, 15)) check("0.0101 * 1.000002^3 < 0.010101", mpf("0.0101") * mpf("1.000002") ** 3 < mpf("0.010101"), mp.nstr(mpf("0.0101") * mpf("1.000002") ** 3, 15)) # Theorem Lower: log p_k >= log 7232121212 > 22.7 and y = 0.01/log^3 p_k < 10^-6 check("log(7232121212) > 22.7", log(mpf(7232121212)) > mpf("22.7"), mp.nstr(log(mpf(7232121212)), 12)) check("0.01/22.7^3 < 10^-6", mpf("0.01") / mpf("22.7") ** 3 < mpf("1e-6"), mp.nstr(mpf("0.01") / mpf("22.7") ** 3, 12)) # --------------------------------------------------------------------------- # Superabundant numbers up to 10**NEXP: which satisfy the SA corollary (a)/(b)? # --------------------------------------------------------------------------- NEXP = int(sys.argv[1]) if len(sys.argv) > 1 else 40 N = 10 ** NEXP print(f"\n=== Superabundant numbers up to 10^{NEXP} ===") HP = primes_upto(400) def hr_integers(limit): """Yield (n, sigma(n)) for all Hardy-Ramanujan integers n <= limit.""" out = [(1, 1)] def rec(idx, max_exp, n, sig): p = HP[idx] pk = 1 sp = 1 # sigma(p^a) for a in range(1, max_exp + 1): pk *= p if n * pk > limit: break sp += pk m, s = n * pk, sig * sp out.append((m, s)) rec(idx + 1, a, m, s) rec(0, 200, 1, 1) return out hr = sorted(hr_integers(N)) print(f"Hardy-Ramanujan integers <= 10^{NEXP}: {len(hr)}") sa = [] best_n, best_s = 1, 1 for n, s in hr: if n == 1 or s * best_n > best_s * n: # I(n) > max_{m 2: break return lp cnt_a = cnt_b = 0 robin_fail = [] last_a = last_b = None eA = exp(mpf("0.0562")) cnt_xneg = 0 cnt_A_fail = 0 last_A_fail = None A_hold = [] xP2_vals = [] x_data = [] for n, s in sa: if n <= 5040: continue pk = largest_prime(n) lln = log(log(mpf(n))) In = mpf(s) / n robin = In < eg * lln a = d * log(theta(pk, HP)) <= lln b = mpf(pk) ** d < log(mpf(n)) if not robin: robin_fail.append(n) if a: cnt_a += 1 last_a = (n, pk) if b: cnt_b += 1 last_b = (n, pk) Pl = log(mpf(pk)) xv = lln - Pl # x = loglog n - log p_k xP2 = xv * Pl ** 2 xP2_vals.append((xP2, n, pk)) x_data.append((n, pk, xv, Pl)) if xv < 0: cnt_xneg += 1 if not (eA > 1 + xP2): # Axler inequality e^0.0562 > 1 + x log^2 p_k fails cnt_A_fail += 1 last_A_fail = (n, pk, xv, xP2) else: A_hold.append((n, pk, xP2)) check(f"Robin's inequality holds for every superabundant 5040 < n <= 10^{NEXP}", not robin_fail, f"failures: {robin_fail[:5]}") print(f"SA n > 5040 with d*log(theta(p_k)) <= loglog n (SA corollary (a)): {cnt_a}" + (f", largest such n = {last_a[0]} (p_k = {last_a[1]})" if last_a else "")) print(f"SA n > 5040 with p_k^d < log n (SA corollary (b)): {cnt_b}" + (f", largest such n = {last_b[0]} (p_k = {last_b[1]})" if last_b else "")) total = len(xP2_vals) print(f"\n--- x = loglog n - log p_k on the {total} superabundant numbers 5040 < n <= 10^{NEXP} ---") print(f"SA n with x < 0 (that is, p_k > log n): {cnt_xneg}") print(f"SA n with e^0.0562 <= 1 + x log^2 p_k (Main inequality FAILS): {cnt_A_fail}" + (f", largest such n = {last_A_fail[0]} (p_k = {last_A_fail[1]}, x log^2 p_k = {mp.nstr(last_A_fail[3], 8)})" if last_A_fail else "")) print(f"SA n with e^0.0562 > 1 + x log^2 p_k (Main inequality HOLDS): {total - cnt_A_fail}") for n_, pk_, v_ in A_hold: print(f" holds: n = {n_} p_k = {pk_} x log^2 p_k = {mp.nstr(v_, 10)}") mx = max(xP2_vals) mn = min(xP2_vals) print(f"max of x log^2 p_k: {mp.nstr(mx[0], 10)} at n = {mx[1]} (p_k = {mx[2]})") print(f"min of x log^2 p_k: {mp.nstr(mn[0], 10)} at n = {mn[1]} (p_k = {mn[2]})") print("x log^2 p_k on the last 10 superabundant numbers in range:") for v, n, pk in xP2_vals[-10:]: print(f" n = {n} p_k = {pk} x log^2 p_k = {mp.nstr(v, 10)}") # --------------------------------------------------------------------------- # Part 5: the hypothesis x >= c / log^alpha(p_k) of the sufficient condition for RH # --------------------------------------------------------------------------- print(f"\n=== Hypothesis x >= c / log^alpha(p_k) on superabundant 5040 < n <= 10^{NEXP} ===") nonpos = [(n, pk, xv) for n, pk, xv, Pl in x_data if xv <= 0] print(f"SA n with x <= 0: {len(nonpos)}") for n_, pk_, xv_ in nonpos: print(f" x <= 0: n = {n_} p_k = {pk_} x = {mp.nstr(xv_, 10)}") N0 = max(n for n, _, _ in nonpos) if nonpos else 5040 tail = [(n, pk, xv, Pl) for n, pk, xv, Pl in x_data if n > N0] print(f"So the hypothesis can only hold with N >= {N0}; beyond it there are {len(tail)} superabundant numbers up to 10^{NEXP}, all with x > 0.") last50 = x_data[-50:] for alpha in range(1, 6): vals = [(xv * Pl ** alpha, n, pk) for n, pk, xv, Pl in tail] m_all = min(vals) m_50 = min((xv * Pl ** alpha, n, pk) for n, pk, xv, Pl in last50) print(f"alpha = {alpha}: min x log^{alpha} p_k over n > N = {N0}: {mp.nstr(m_all[0], 10)} " f"(n = {m_all[1]}, p_k = {m_all[2]}); over the last 50 SA numbers: {mp.nstr(m_50[0], 10)} (p_k = {m_50[2]})") print("x log^alpha p_k on the last 10 superabundant numbers in range (alpha = 1, 2, 3):") for n, pk, xv, Pl in x_data[-10:]: print(f" p_k = {pk} x = {mp.nstr(xv, 8)} x log p_k = {mp.nstr(xv * Pl, 8)} " f"x log^2 p_k = {mp.nstr(xv * Pl ** 2, 8)} x log^3 p_k = {mp.nstr(xv * Pl ** 3, 8)}") print(f"\n{sum(results)}/{len(results)} checks passed") sys.exit(0 if all(results) else 1)