"""
Independent audit of the analytic steps behind the certificate of Lambda <= 0.174 (Paper IX, Sections 3-4 and
Appendix A).  This script is NOT part of the certificate; it checks, with mpmath and plain numpy, that

  1. the inputs of the certificate are the files that were certified (SHA-256 against CERTIFICATE_0.174.txt);
  2. the algebra of Lemma 3 is right: the coefficientwise expansion of M*conj(B), the identity
     |M|(|A|-|gamma||B|) = |MA| - |gamma||M conj(B)|, and the final inequality, at points where s*, kappa, gamma
     are computed exactly from Polymath's definitions (17)-(18) (random weights, moderate x);
  3. the side conditions used in the code hold: D^{|kappa|}-1 <= 1e-12 for x >= X, 3^y+3^{-y} <= 3*y_b+1 on the
     y-range of hypothesis (ii) and <= 4 on [0,1], sigma(x_v) > 0 and e_y(v) < 0 for v >= 2;
  4. the derivative bounds in the proof of Lemma 5 (monotonicity of F) and of Lemma 6 (tail) have the stated sign
     on a logarithmic grid, and F itself decreases on that grid (the proof is analytic; the grid is a sanity check);
  5. a plain floating-point evaluation of the lower bound of Lemma 3 at the first chunk of hypothesis (ii)
     (x = X + sqrt(1-y0^2), y in [0.2, 0.205], N = 690891, D = 240), written independently of cert_region_ii.py
     and without its block bounds, agrees with (and is slightly larger than) the certified value 0.02188.

Run:  python audit_lemmas.py            (about 3 minutes; needs numpy and mpmath)
P. Caceres, October 2026.
"""
import hashlib
import math
import os
import re

import mpmath as mp
import numpy as np

HERE = os.path.dirname(os.path.abspath(__file__))
T0, Y0, X = 0.15399, 0.2, 5998313329980.3135


def check_hashes():
    cert = open(os.path.join(HERE, "CERTIFICATE_0.174.txt"), encoding="utf-8").read()
    pairs = re.findall(r"(\S+\.(?:py|npy)) sha256=([0-9a-f]{64})", cert)
    for name, h in pairs:
        got = hashlib.sha256(open(os.path.join(HERE, name), "rb").read()).hexdigest()
        print("  %-18s %s" % (name, "matches the certificate" if got == h else "DIFFERS: " + got))
        assert got == h
    return len(pairs)


def alpha(s):
    return 1 / (2 * s) + 1 / (s - 1) + mp.log(s / (2 * mp.pi)) / 2


def log_m0(s):
    return (mp.log(s) + mp.log(s - 1) - (s / 2) * mp.log(mp.pi) + mp.log(mp.sqrt(2 * mp.pi) / 16)
            + (s / 2 - mp.mpf(1) / 2) * mp.log(s / 2) - s / 2)


def polymath_parameters(x, y, t):
    s1 = (1 + y - 1j * x) / 2
    s2 = (1 - y + 1j * x) / 2
    s_star = s1 + (t / 2) * alpha(s1)
    kappa = (t / 2) * (alpha(s2) - alpha((1 + y + 1j * x) / 2))
    gamma = mp.exp((t / 4) * (alpha(s2) ** 2 - alpha(s1) ** 2) + log_m0(s2) - log_m0(s1))
    return s_star, kappa, gamma


def check_lemma3(trials=6):
    """Lemma 3 at exact Polymath parameters, random weights; returns the worst slack."""
    mp.mp.dps = 40
    rng = np.random.default_rng(1)
    worst = mp.inf
    for k in range(trials):
        x = mp.mpf(800 + 900 * k)
        y = mp.mpf(0.2 + 0.1 * k) if k < 3 else mp.mpf("0.95")
        t = mp.mpf("0.15") if k < 3 else mp.mpf("0.5")
        s, kap, gam = polymath_parameters(x, y, t)
        n_max = int(mp.floor(mp.sqrt(x / (4 * mp.pi)) + t / 16))
        d_max = 6
        if k < 3:
            lam = [mp.mpf(1)] + [mp.mpf(float(v)) for v in rng.normal(0, 0.6, d_max - 1)]
        else:   # Euler-type weights -b_p on 2, 3, 5 plus a random perturbation: here the bound L is positive
            lam = [mp.mpf(1), -mp.e ** ((t / 4) * mp.log(2) ** 2), -mp.e ** ((t / 4) * mp.log(3) ** 2), mp.mpf(0),
                   -mp.e ** ((t / 4) * mp.log(5) ** 2), mp.mpf(0)]
            lam = [lam[0]] + [v + mp.mpf(float(r)) for v, r in zip(lam[1:], rng.normal(0, 0.05, d_max - 1))]
        b = [None] + [mp.e ** ((t / 4) * mp.log(n) ** 2) for n in range(1, n_max + 1)]
        big_a = mp.fsum(b[n] * mp.power(n, -s) for n in range(1, n_max + 1))
        big_b = mp.fsum(mp.power(n, y) * b[n] * mp.power(n, -mp.conj(s) - kap) for n in range(1, n_max + 1))
        mm = mp.fsum(lam[d - 1] * mp.power(d, -s) for d in range(1, d_max + 1))
        # coefficientwise expansion of M conj(B)
        expansion = mp.mpc(0)
        for n in range(1, d_max * n_max + 1):
            inner = mp.fsum(lam[d - 1] * mp.power(d, -y + mp.conj(kap)) * b[n // d]
                            for d in range(1, d_max + 1) if n % d == 0 and n // d <= n_max)
            expansion += mp.power(n, y - s - mp.conj(kap)) * inner
        assert abs(expansion - mm * mp.conj(big_b)) < mp.mpf(10) ** -25, "expansion of M conj(B) fails"
        lhs = abs(mm) * (abs(big_a) - abs(gam) * abs(big_b))
        rhs = abs(mm * big_a) - abs(gam) * abs(mm * mp.conj(big_b))
        assert abs(lhs - rhs) < mp.mpf(10) ** -25
        sigma = mp.re(s)
        c = {}
        cp = {}
        eps = {}
        dk = mp.power(d_max, abs(kap)) - 1
        for d in range(1, d_max + 1):
            for m in range(1, n_max + 1):
                n = d * m
                c[n] = c.get(n, 0) + lam[d - 1] * b[m]
                cp[n] = cp.get(n, 0) + lam[d - 1] * mp.power(d, -y) * b[m]
                eps[n] = eps.get(n, 0) + abs(lam[d - 1]) * mp.power(d, -y) * b[m] * dk
        bound = (1 - mp.fsum(abs(c[n]) * mp.power(n, -sigma) for n in c if n >= 2)
                 - abs(gam) * mp.fsum((abs(cp[n]) + eps[n]) * mp.power(n, y - sigma + abs(kap)) for n in cp))
        m_plus = mp.fsum(abs(lam[d - 1]) * mp.power(d, -sigma) for d in range(1, d_max + 1))
        f_t = big_a + gam * big_b
        slack = abs(f_t) - bound / m_plus
        print("  x=%5d y=%.2f t=%.2f N=%3d: |f_t|=%.4f  L/M+=%+.4f  |M|(|A|-|g||B|) - L = %.3e" % (
            int(x), float(y), float(t), n_max, float(abs(f_t)), float(bound / m_plus), float(lhs - bound)))
        assert lhs >= bound and slack >= 0
        worst = min(worst, slack)
    return worst


def check_side_conditions():
    iv = mp.iv
    iv.dps = 30
    t = iv.mpf("0.15399")
    x = iv.mpf("5998313329980.3135")
    kap = t * 1 / (2 * (x - 6))                       # y <= 1
    dk = iv.exp(kap * iv.log(240)) - 1
    print("  D^|kappa| - 1 <= %.3e  (code uses 1e-12)" % float(dk.b))
    assert dk.b < 1e-12
    y_top = float(mp.sqrt(1 - 2 * mp.mpf("0.15399")))
    ys = np.linspace(0, y_top, 2001)
    assert np.all(3 ** ys + 3.0 ** -ys <= 3 * y_top + 1 - 0.5)
    assert np.all(3 ** np.linspace(0, 1, 2001) + 3.0 ** -np.linspace(0, 1, 2001) <= 4)
    print("  3^y+3^-y <= %.4f < 3*y_b+1 = %.4f on [0, y_b];  <= 10/3 < 4 on [0,1]" % (3 ** y_top + 3 ** -y_top, 3 * y_top + 1))


def sigma_lb(x, y, t):
    return (1 + y) / 2 + (t / 4) * mp.log(x / (4 * mp.pi)) - t / x ** 2


def gamma_ub(x, y):
    return mp.exp(mp.mpf("0.02") * y) * mp.exp(-(y / 2) * mp.log(x / (4 * mp.pi)))


def kappa_ub(x, y, t):
    return t * y / (2 * (x - 6))


def big_f(v, ya, yb, t, lam_abs, logd):
    """F(v) of Lemma 5 (including the factor D^{kappa_+})."""
    xv = 4 * mp.pi * (v - t / 16) ** 2
    lv = mp.log(v)
    s = sigma_lb(xv, ya, t)
    e = max(ya - sigma_lb(xv, ya, t) + kappa_ub(xv, ya, t), yb - sigma_lb(xv, yb, t) + kappa_ub(xv, yb, t))
    bv = mp.exp((t / 4) * lv ** 2)
    s1 = mp.fsum(a * mp.exp(-s * l) for a, l in zip(lam_abs, logd))
    s2 = mp.fsum(a * mp.exp((e - ya) * l) for a, l in zip(lam_abs, logd))
    dk = mp.exp(kappa_ub(xv, yb, t) * logd[-1])
    return bv * mp.exp(-s * lv) * s1 + gamma_ub(xv, ya) * dk * bv * mp.exp(e * lv) * s2


def check_lemma5(lam):
    mp.mp.dps = 30
    t = mp.mpf("0.15399")
    lam_abs = [mp.mpf(abs(float(v))) for v in lam]
    logd = [mp.log(d) for d in range(1, len(lam) + 1)]
    for ya, yb in ((0.2, 0.205), (0.515, 0.831877)):
        ya, yb = mp.mpf(ya), mp.mpf(yb)
        grid = [mp.mpf(2) * mp.mpf(1.07) ** k for k in range(0, 2100, 7)]      # v from 2 to about e^{140}
        vals = [big_f(v, ya, yb, t, lam_abs, logd) for v in grid]
        assert all(b < a for a, b in zip(vals, vals[1:])), "F not decreasing on the grid"
        worst1 = worst2 = -mp.inf
        for v in grid:
            xv = 4 * mp.pi * (v - t / 16) ** 2
            lv, ell = mp.log(v), mp.log(v - t / 16)
            sv = sigma_lb(xv, ya, t)
            # Lemma 5: phi_1'(v) <= (t/2) L (1/v - 1/(v-t/16)) - sigma(x_v)/v
            d1 = (t / 2) * lv * (1 / v - 1 / (v - t / 16)) - sv / v
            # Lemma 5: psi_y'(v) <= (t/2) L (1/v - 1/(v-t/16)) + e_y(v)/v
            ey = max(y - sigma_lb(xv, y, t) + kappa_ub(xv, y, t) for y in (ya, yb))
            d2 = (t / 2) * lv * (1 / v - 1 / (v - t / 16)) + ey / v
            assert sv > 0 and ey < 0 and d1 < 0 and d2 < 0
            worst1, worst2 = max(worst1, d1 * v), max(worst2, d2 * v)
            if lv >= 60:   # Lemma 6: v d/dv log(T_i v^2) <= 2 - (t/2) log(v - t/16) + t/x_v^2 + 1/10^6
                assert 2 - (t / 2) * ell + t / xv ** 2 + mp.mpf(10) ** -6 < 0
        print("  y in [%.3f,%.3f]: F decreasing on %d grid points in [2, e^140];  max v*phi1' = %.3f, max v*psi' = %.3f"
              % (float(ya), float(yb), len(grid), float(worst1), float(worst2)))
        vmax = mp.e ** 60
        chk = big_f(vmax, ya, yb, t, lam_abs, logd) * vmax ** 2
        print("      F(e^60) e^120 = %.3e < 1  (Lemma 6)" % float(chk))
        assert chk < 1


def float_lower_bound(lam, ya=0.2, yb=0.205):
    """Plain float64 evaluation of the bound of Lemma 3 at the first chunk (no block bounds, no outward rounding)."""
    t = T0
    x = X + math.sqrt(1 - Y0 ** 2)
    n_trunc = int(math.floor(math.sqrt(x / (4 * math.pi)) + t / 16))
    sig = (1 + ya) / 2 + (t / 4) * math.log(x / (4 * math.pi)) - t / x ** 2
    e_b = max(y - ((1 + y) / 2 + (t / 4) * math.log(x / (4 * math.pi)) - t / x ** 2) + t * y / (2 * (x - 6))
              for y in (ya, yb))
    g = math.exp(0.02 * ya) * (x / (4 * math.pi)) ** (-ya / 2)
    d_count = len(lam)
    m = np.arange(1, n_trunc + 1, dtype=np.float64)
    b = np.exp(t / 4 * np.log(m) ** 2)
    dd = np.arange(1, d_count + 1, dtype=np.float64)
    dya, dyb = dd ** -ya, dd ** -yb
    total_a = total_b = 0.0
    chunk = 20_000_000
    n_max = d_count * n_trunc
    for lo_n in range(1, n_max + 1, chunk):
        hi_n = min(lo_n + chunk, n_max + 1)
        ca = np.zeros(hi_n - lo_n)
        cb = np.zeros(hi_n - lo_n)
        ce = np.zeros(hi_n - lo_n)
        for j in range(d_count):
            d = j + 1
            m0 = max(1, -(-lo_n // d))
            m1 = min(n_trunc, (hi_n - 1) // d)
            if m1 < m0:
                continue
            idx = d * np.arange(m0, m1 + 1) - lo_n
            ca[idx] += lam[j] * b[m0 - 1:m1]
            cb[idx] += lam[j] * dya[j] * b[m0 - 1:m1]
            ce[idx] += abs(lam[j]) * (dya[j] - dyb[j] + 1e-12 * dya[j]) * b[m0 - 1:m1]
        n = np.arange(lo_n, hi_n, dtype=np.float64)
        if lo_n == 1:
            ca[0] = 0.0
        total_a += float(np.sum(np.abs(ca) * n ** -sig))
        total_b += float(np.sum((np.abs(cb) + ce) * n ** e_b))
    lower = 1 - total_a - g * total_b
    m_plus = float(np.sum(np.abs(lam) * dd ** -sig))
    return lower, m_plus, n_trunc


if __name__ == "__main__":
    lam = np.load(os.path.join(HERE, "lam_L0174.npy"))
    print("1. Inputs")
    print("  %d files checked; weights: %d float64 numbers, lambda_1 = %r" % (check_hashes(), len(lam), lam[0]))
    print("2. Lemma 3 at exact Polymath parameters (D = 6; random weights, and perturbed Euler weights where L > 0)")
    print("  smallest slack |f_t| - L/M+ = %.4f" % float(check_lemma3()))
    print("3. Side conditions")
    check_side_conditions()
    print("4. Lemmas 5 and 6 (monotonicity of F, tail)")
    check_lemma5(lam)
    print("5. Float evaluation of L at the first chunk of (ii)")
    lower, m_plus, n_trunc = float_lower_bound(lam)
    print("  N = %d, L = %.5f (certified, with block bounds: 0.02188), |M| <= %.4f (certified 1.9646)" % (n_trunc, lower, m_plus))
    assert lower > 0.02188
    print("All audit checks passed.")