"""Verification script for "Robin's Inequality at Colossally Abundant Numbers and
Nicolas' Inequality at Their Largest Prime Factor" (F. Vega).

Written with the help of Claude Opus 5.5 (Anthropic), as was the script of the earlier
version of this work; run by the author.

Every numerical inequality between explicit constants used in the paper is checked
with mpmath at 60 significant digits; the computed left-hand side is printed.

Part 1  constants of Lemma Pboot, Lemma LlogL, Lemma corekernel and Theorem nicSA
        (counterexamples and superabundant counterexamples).
Part 2  constants of Lemmas theta, S2, T1, IN and of the proof of the main theorem
        (colossally abundant numbers with log xi >= 31).
Part 3  constants of the corollaries and of Proposition NicRH.

Self-contained: needs only Python 3.12 and mpmath (no package install).
Run:  python verify_robin_nicolas.py
"""

import sys
from fractions import Fraction as D
from math import gcd
from mpmath import mp, mpf, euler, exp, log, pi, sqrt

mp.dps = 60

results = []


def check(label, ok, value=None):
    results.append(bool(ok))
    detail = f"   ({mp.nstr(value, 15)})" if value is not None else ""
    print(f"[{'PASS' if ok else 'FAIL'}] {label}{detail}")


def eps(P):
    """Axler's epsilon(P) = 1/(20 P^3) + 3/(16 P^4)."""
    return 1 / (20 * P ** 3) + mpf(3) / (16 * P ** 4)


eg = exp(euler)
L31 = mpf(31)

# ---------------------------------------------------------------------------
print("=== Part 1: counterexamples and superabundant counterexamples ===")
LMP = mpf("13.11485") * log(10) + log(log(10))
check("loglog(10^(10^13.11485)) > 31.032", LMP > mpf("31.032"), LMP)
x0 = mpf(46909038)
check("log(46909038) < 17.6638", log(x0) < mpf("17.6638"), log(x0))
v = 1 / (1 - eps(log(x0)))
check("(1 - eps(log 46909038))^{-1} < 1.000011", v < mpf("1.000011"), v)
check("17.6638 * 1.000011 < 31.032", D("17.6638") * D("1.000011") < D("31.032"))
v = mpf("31.032") / mpf("1.000011")
check("31.032/1.000011 > 31.0316", v > mpf("31.0316"), v)
v = 1 / (1 - eps(mpf("31.0316")))
check("(1 - eps(31.0316))^{-1} < 1.000002", v < mpf("1.000002"), v)
v = mpf("31.032") / mpf("1.000002")
check("31.032/1.000002 > 31.0319", v > mpf("31.0319"), v)
P1 = mpf("31.0319")
v = (mpf("0.05") + mpf(3) / (16 * P1)) / (1 - eps(P1))
check("(1/20 + 3/(16*31.0319))/(1 - eps(31.0319)) < 0.0562", v < mpf("0.0562"), v)
Q0 = mpf("31.03")
v = exp(mpf("0.0562") / Q0 ** 2) * (1 + mpf("0.0562") / Q0 ** 3)
check("e^(0.0562/31.03^2) (1 + 0.0562/31.03^3) < 1.0000603", v < mpf("1.0000603"), v)
v = exp(mpf("15.515"))
check("e^15.515 > 5.47e6", v > mpf("5.47e6"), v)
v = 1 - mpf("1.94") / sqrt(mpf("5.47e6"))
check("1 - 1.94/sqrt(5.47e6) > 0.99917", v > mpf("0.99917"), v)
v = mpf("0.99917") / mpf("1.0000603")
check("0.99917/1.0000603 > 0.9991", v > mpf("0.9991"), v)
v = eg * mpf("31.032")
check("e^gamma * 31.032 > 55", v > 55, v)
phi = lambda m: sum(1 for j in range(1, m + 1) if gcd(j, m) == 1)
v = max(mpf(m) / phi(m) for m in range(1, 16))
check("max m/phi(m) over 1 <= m <= 15 equals 3", v == 3, v)

# ---------------------------------------------------------------------------
print("\n=== Part 2: colossally abundant numbers with log xi >= 31 ===")
h = log(2) / 2
check("e^31 > 10^9 + 7 and e^31 > 89967803", exp(L31) > 10 ** 9 + 7 and exp(L31) > 89967803, exp(L31))
v = 1 - log(2) / (2 * L31)
check("1 - log 2/(2*31) > 0.98882", v > mpf("0.98882"), v)
v = sqrt(2) * mpf("0.98882")
check("sqrt(2) * 0.98882 > 1.39838", v > mpf("1.39838"), v)
v = sqrt(mpf("1.39838"))
check("sqrt(1.39838) > 1.18253", v > mpf("1.18253"), v)
v = mpf("1.94") * L31 * exp(-L31 / 4) / mpf("1.18253")
check("1.94 * 31 * e^(-31/4) / 1.18253 < 0.022", v < mpf("0.022"), v)
v = 1 / log(mpf(10) ** 19) ** 3
check("1/log^3(10^19) < 0.000012", v < mpf("0.000012"), v)
v = 2 / log(mpf(10) ** 19) ** 2
check("2 log(10^19)/log^3(10^19) < 0.022", v < mpf("0.022"), v)
v = mpf(8) / L31 ** 2
check("8/31^2 < 0.022", v < mpf("0.022"), v)
v = -log(1 - log(2) / (2 * L31))
check("-log(1 - log 2/62) < 0.011306", v < mpf("0.011306"), v)
v = 2 * (h - mpf("0.011306"))
check("2*(log 2/2 - 0.011306) >= 0.670535", v >= mpf("0.670535"), v)
v = log(2) - mpf("0.323")
check("log 2 - 0.323 < 0.370148", v < mpf("0.370148"), v)
v = mpf("0.670535") - h - h * mpf("0.670535") / L31
check("0.670535 - log2/2 - (log2/2)(0.670535)/31 > 0.316464", v > mpf("0.316464"), v)
v = mpf("0.316464") - mpf("0.316464") ** 2 / L31
check("0.316464 - 0.316464^2/31 > 0.313233", v > mpf("0.313233"), v)
v = 2 * (1 - (mpf("0.370148") + log(2)) / L31)
check("2*(1 - (0.370148 + log 2)/31) > 1.9314", v > mpf("1.9314"), v)
v = D("0.066") + D("0.370148") + 2
check("0.066 + 0.370148 + 2 <= 2.436148", v <= D("2.436148"), mpf(v.numerator) / v.denominator)
v = mpf("1.0022") * L31 * exp(-L31 / 2) / sqrt(2)
check("1.0022 * 31 * e^(-31/2)/sqrt(2) < 4.08e-6", v < mpf("4.08e-6"), v)
check("1 + 3*0.022/31 < 1.0022", 1 + 3 * mpf("0.022") / L31 < mpf("1.0022"))
check("2.436148 + 4.08e-6 < 2.43616", D("2.436148") + D("4.08e-6") < D("2.43616"))
v = 8 * (1 + mpf("0.022") / L31)
check("8*(1 + 0.022/31) < 8.006", v < mpf("8.006"), v)
v = mpf("8.006") / L31
check("8.006/31 < 0.258259", v < mpf("0.258259"), v)
v = D("0.313233") - D("0.066") + D("1.9314") - D("0.258259")
check("0.313233 - 0.066 + 1.9314 - 0.258259 >= 1.920374", v >= D("1.920374"), mpf(v.numerator) / v.denominator)
check("0.022 + log 2/2 <= 0.368574", mpf("0.022") + h <= mpf("0.368574"), mpf("0.022") + h)
check("0.323 - 0.022 >= 0.301", D("0.323") - D("0.022") >= D("0.301"))
check("e^(31/3) > 853 and > 1418", exp(L31 / 3) > 1418, exp(L31 / 3))
check("0.945 * e^(31/3) > 31", mpf("0.945") * exp(L31 / 3) > L31, mpf("0.945") * exp(L31 / 3))
v = mpf("1.000012") * mpf("1.9769") / sqrt(2)
check("1.000012 * 1.9769/sqrt(2) < 1.397905", v < mpf("1.397905"), v)
v = mpf("1.397905") * L31 * exp(-L31 / 6)
check("1.397905 * 31 * e^(-31/6) < 0.24717", v < mpf("0.24717"), v)
check("0.301 - 0.24717 >= 0.0538", D("0.301") - D("0.24717") >= D("0.0538"))
v = 1 / L31 ** 2 + sqrt(2) * L31 * exp(-L31 / 2)
check("1/31^2 + sqrt(2)*31*e^(-31/2) < 0.0010488", v < mpf("0.0010488"), v)
v = mpf("0.0010488") / L31
check("0.0010488/31 < 0.0000339", v < mpf("0.0000339"), v)
v = mpf("0.0010488") + mpf("0.0000339") + mpf("0.0010488") * mpf("0.0000339") / L31
check("0.0010488 + 0.0000339 + (0.0010488*0.0000339)/31 < 0.001085", v < mpf("0.001085"), v)
v = 1 / L31 ** 2 + mpf("1.0001") / L31 ** 3
check("1/31^2 + 1.0001/31^3 < 0.00108", v < mpf("0.00108"), v)
v = 1 / (1 - L31 ** -3)
check("1/(1 - 31^-3) < 1.0001 (so -log(1 - x) <= 1.0001 x)", v < mpf("1.0001"), v)
check("0.368574 + 0.001085 <= 0.369659", D("0.368574") + D("0.001085") <= D("0.369659"))
check("0.0538 - 0.00108 >= 0.05272", D("0.0538") - D("0.00108") >= D("0.05272"))
v = mpf("0.086") / (1 - mpf(1418) ** -3)
check("0.086/(1 - 1418^-3) < 0.0861", v < mpf("0.0861"), v)
v = 3 * mpf("1.26") * mpf(3) ** (mpf(1) / 3)
check("3 * 1.26 * 3^(1/3) < 5.4518", v < mpf("5.4518"), v)
v = mpf("5.4518") / L31 * (1 + mpf(10) ** -12)
check("(5.4518/31)(1 + 1e-12) < 0.17587", v < mpf("0.17587"), v)
v = mpf("0.0861") + mpf("0.17587") + exp(-L31 / 3)
check("0.0861 + 0.17587 + e^(-31/3) < 0.26201", v < mpf("0.26201"), v)
v = mpf("0.26201") * L31 ** 2 * exp(-L31 / 6) / sqrt(2)
check("0.26201 * 31^2 * e^(-31/6)/sqrt(2) < 1.0155", v < mpf("1.0155"), v)
v = 2 * mpf("1.39838") ** -2 * exp(-L31) * L31
check("2 xi_2^-2 * 31 < 1e-4 with xi_2 >= 1.39838 e^(31/2)", v < mpf("1e-4"), v)
check("1.0155 + 1e-4 <= 1.0156", D("1.0155") + D("1e-4") <= D("1.0156"))
check("0.369659 + 2.43616 <= 2.805819", D("0.369659") + D("2.43616") <= D("2.805819"))
check("0.05272 + 1.920374 - 1.0156 >= 0.957494",
      D("0.05272") + D("1.920374") - D("1.0156") >= D("0.957494"))
v = mpf("0.5") / L31 ** 2 + 1 / L31 ** 3
check("1/(2*31^2) + 1/31^3 < 0.000554", v < mpf("0.000554"), v)
v = mpf("0.957494") * (1 - L31 ** -4)
check("0.957494 (1 - 31^-4) > 0.9574", v > mpf("0.9574"), v)
check("2.805819 + 2*0.000554 < 2.8070", D("2.805819") + 2 * D("0.000554") < D("2.8070"))
v = mpf("0.9574") * sqrt(2)
check("0.9574 sqrt(2) > 1.353", v > mpf("1.353"), v)
v = mpf("2.8070") * sqrt(2)
check("2.8070 sqrt(2) < 3.970", v < mpf("3.970"), v)
check("0.9574/2 >= 0.4787 and 2.8070/2 <= 1.4035",
      D("0.9574") / 2 >= D("0.4787") and D("2.8070") / 2 <= D("1.4035"))

# ---------------------------------------------------------------------------
print("\n=== Part 3: corollaries ===")
v = 2 * sqrt(2) * (1 - mpf("1.4035") / mpf("31.0319"))
check("2 sqrt(2) (1 - 1.4035/31.0319) > 2.70", v > mpf("2.70"), v)
check("e^31.0319 > 10^9", exp(mpf("31.0319")) > mpf(10) ** 9, exp(mpf("31.0319")))
tau = 2 + euler - log(4 * pi)
check("tau = 2 + gamma - log(4 pi) < 0.0462", tau < mpf("0.0462"), tau)
check("2 + tau < 2.0462 < 2.70", 2 + tau < mpf("2.0462") < mpf("2.70"), 2 + tau)
v = 2 + tau - 2 * sqrt(2) + 2 * sqrt(2) * mpf("1.4035") / L31
check("2 + tau - 2 sqrt(2) + 2 sqrt(2) 1.4035/31 < -0.654", v < mpf("-0.654"), v)
v = 2 * sqrt(2) - 2 - tau
print(f"      2 sqrt(2) - 2 - tau = {mp.nstr(v, 15)}")

print(f"\n{sum(results)}/{len(results)} checks passed")
sys.exit(0 if all(results) else 1)