# VERIF_FERMAT_69.py
# Ancillary file for the paper "The exponent of p in x**(p-1) - 1: Yamada's constant 283 lowered to 69".
# David Niedbala Giraudin, ORCID 0009-0009-1526-1178.
# Technical assistance: Claude (Anthropic).
# Plain Python + numpy. No argument, no external file. Run in plain Python (not under Sage).
# Running time: a few seconds. Memory: a few MB.
#
# Contents:
#   [0] arithmetic self-test
#   [1] Lemma 3 (rearrangement): exact bound on the spread, checked by brute force, and attained
#   [2] Lower bound for the product of factorials (Lemma 5), checked up to K = 20000
#   [3] Closed bound for beta = (log b)/(log p) on the whole non-trivial range, and control grid
#   [4] Final condition of Theorem 1 (C = 69, L = 5): margin, residual terms and bound for eta(K)
#   [5] Cross-check at explicit primes p, with the integer parameters built exactly
#   [6] Calibration: the same chain with Yamada's coefficient L gives back 282.79 <= 283; table C(L)
#   [7] Degenerate case (Section 4.2): (R2 - 1) log x / log p <= L sqrt(k) g a1 a2 < 69 g a1 a2

import decimal
import math
import random
import numpy as np

LOG2 = math.log(2.0)


def say(*a):
    print(*a, flush=True)


# ---------------------------------------------------------------- [0]
def selftest():
    assert 7 * 6 == 42
    assert int(np.sum(np.arange(10, dtype=np.int64))) == 45
    assert abs(math.log(float(math.e)) - 1.0) < float(1e-15)
    assert math.isqrt(10 ** 20) == 10 ** 10
    assert pow(3, 5) == 243
    say("[0] arithmetic self-test: OK")


# ---------------------------------------------------------------- [1]
def lemma_R():
    # l_j = floor((j-1)/K), j = 1..N, N = K L ; r_j integers in [0, R-1].
    # max_sigma | sum l_i r_sigma(i) - (L-1)/2 sum r | <= (R-1) K floor(L*L/4) / 2.
    rng = random.Random(20261022)
    worst = 0.0
    n_tests = 0
    for K in range(1, 7):
        for L in range(1, 8):
            N = K * L
            l = np.repeat(np.arange(L, dtype=np.float64), K)
            for R in (1, 2, 3, 5, 9, 17):
                G = (R - 1) * K * ((L * L) // 4) / 2.0
                for t in range(60):
                    r = np.array([rng.randrange(R) for _ in range(N)], dtype=np.float64)
                    M = (L - 1) / 2.0 * float(np.sum(r))
                    rs = np.sort(r)
                    hi = float(np.dot(l, rs))          # rearrangement: maximum
                    lo = float(np.dot(l, rs[::-1]))    # minimum
                    dev = max(hi - M, M - lo)
                    assert dev <= G + float(1e-9), (K, L, R, dev, G)
                    if G > 0:
                        worst = max(worst, dev / G)
                    n_tests += 1
                # extremal configuration: half at 0, half at R-1 -> equality (L even, K even)
                if L % 2 == 0 and K % 2 == 0 and R > 1:
                    r = np.array([0.0] * (N // 2) + [float(R - 1)] * (N // 2))
                    M = (L - 1) / 2.0 * float(np.sum(r))
                    hi = float(np.dot(l, np.sort(r)))
                    assert abs((hi - M) - G) < float(1e-9)
    say("[1] Lemma 3: %d draws, maximal ratio spread/bound %.6f ; equality attained (L, K even): OK"
        % (n_tests, worst))


# ---------------------------------------------------------------- [2]
def lemma_F():
    # (2/(K*K-K)) sum for k from 1 to K-1 of log k! >= log(K-1) - 3/2 - log(K-1)/K   (K >= 2)
    s = 0.0
    worst = float("inf")
    for K in range(2, 20001):
        s += math.lgamma(K)          # adds log (K-1)!
        lhs = 2.0 * s / (K * K - K)
        rhs = math.log(K - 1) - 1.5 - math.log(K - 1) / K
        worst = min(worst, lhs - rhs)
        assert lhs >= rhs - float(1e-12), K
    say("[2] Lemma 5: checked for 2 <= K <= 20000, minimal margin %.3e" % worst)


# ---------------------------------------------------------------- [3]
C = 69.0
L = 5
k = C / (L * L)
CF = 2.0 * ((L * L) // 4) / L       # 12/5
EPS = float(1e-8)                   # bounds eps(K) of BL96 Lemma 13 and log(K-1)/K (K > 4e18)


def beta_closed():
    c1 = 1.5 - LOG2 - 0.5 * math.log(k) + EPS
    A = math.log(C) + c1
    assert A > 2.0
    return 1.0 + A / (C * LOG2), A


def beta_grid():
    c1 = 1.5 - LOG2 - 0.5 * math.log(k) + EPS
    worst = 0.0
    arg = None
    for y in range(2, 300):
        ly = math.log(y)
        for x in list(range(2, 3000)) + [10 ** e + 1 for e in range(4, 60)]:
            lx = math.log(x)
            lxy = lx + ly
            P0 = C * ly * lxy / lx
            E = math.log(1.0 / lxy + 1.0 / ly) + c1
            # sup over P >= P0 of (log P + E)/P
            if math.log(P0) + E >= 1.0:
                v = (math.log(P0) + E) / P0
            else:
                v = math.exp(E - 1.0)
            if v > worst:
                worst, arg = v, (x, y)
    return 1.0 + worst, arg


# ---------------------------------------------------------------- [4]
def final_condition(beta):
    Gam0 = (math.exp(C * LOG2) - 1.0) * LOG2 * math.log(6.0) / (C * LOG2) ** 2
    # residual terms, normalised by Gamma * log p (see Section 4)
    P = C * LOG2
    Nmax_log = math.log(k * L * L * Gam0 + L)            # log N, increasing in Gamma: a crude bound is taken
    d1 = 3.0 * (Nmax_log + 50.0) / (Gam0 * P)            # 3 log N / (Gamma log p), generous overestimate
    d2 = 2.0 * CF * math.sqrt(L / Gam0)                   # terms R1 - 1, S1 - 1
    delta = d1 + d2
    lhs = k * L * (L - 1 - beta)
    rhs = 2.0 * CF * L * math.sqrt(k)
    say("[4] k = %.4f, L = %d, coefficient = %.4f, Gamma >= %.3e" % (k, L, CF, Gam0))
    say("    k L (L-1-beta) = %.6f   2 c L sqrt(k) = %.6f   residual <= %.2e" % (lhs, rhs, delta))
    m = lhs - rhs - delta
    say("    margin = %.6f  -> %s" % (m, "OK" if m > 0 else "FAILED"))
    assert m > 0
    # bound for eta(K) = log(sqrt(K)(1+sqrt(K-1))/(K-1)) + log(K-1)/K for K >= k L Gamma0 (eta decreasing)
    Kmin = k * L * Gam0
    eta = math.log(math.sqrt(Kmin) * (1.0 + math.sqrt(Kmin - 1.0)) / (Kmin - 1.0)) + math.log(Kmin - 1.0) / Kmin
    for Kt in (Kmin, 10.0 * Kmin, 1000.0 * Kmin):
        et = math.log(math.sqrt(Kt) * (1.0 + math.sqrt(Kt - 1.0)) / (Kt - 1.0)) + math.log(Kt - 1.0) / Kt
        assert et <= eta + float(1e-15)
    say("    eta(K) <= %.3e < 1e-8 for K >= %.3e : %s" % (eta, Kmin, "OK" if eta < EPS else "FAILED"))
    assert eta < EPS
    # final integer: N - 1 = K L - 1 <= k L L Gamma + L - 1 = 69 Gamma + 4
    assert abs(k * L * L - 69.0) < float(1e-12)


# ---------------------------------------------------------------- [5]
def is_prime(n):
    if n < 2:
        return False
    for q in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37):
        if n % q == 0:
            return n == q
    d, s = n - 1, 0
    while d % 2 == 0:
        d //= 2
        s += 1
    for a in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41):
        x = pow(a, d, n)
        if x == 1 or x == n - 1:
            continue
        for _ in range(s - 1):
            x = pow(x, 2, n)
            if x == n - 1:
                break
        else:
            return False
    return True


def next_prime(n):
    n |= 1
    while not is_prime(n):
        n += 2
    return n


def direct_instance(p, x, y):
    # exact construction of the integer parameters in decimal arithmetic (3 times the number of digits of p, plus 80)
    D = decimal.Decimal
    decimal.getcontext().prec = 3 * len(str(p)) + 80
    P = D(p).ln()
    g = p - 1
    lA1, lA2 = D(x * y).ln(), D(y).ln()
    a1, a2 = lA1 / P, lA2 / P
    Gam = D(g) * a1 * a2
    K = int((D(69) / D(25) * L * Gam).to_integral_value(rounding=decimal.ROUND_FLOOR)) + 1
    fl = lambda t: int(t.to_integral_value(rounding=decimal.ROUND_FLOOR))
    R1 = fl((D(g) * L * a2 / a1).sqrt()) + 1
    S1 = fl((D(g) * L * a1 / a2).sqrt()) + 1
    R2 = fl((D(g) * (K - 1) * L * a2 / a1).sqrt()) + 1
    S2 = fl((D(g) * (K - 1) * L * a1 / a2).sqrt()) + 1
    assert R1 * S1 > g * L and R2 * S2 > g * (K - 1) * L, (p, x, y, R1 * S1 - g * L, R2 * S2 - g * (K - 1) * L)     # pigeonhole over the classes
    R, S, N = R1 + R2 - 1, S1 + S2 - 1, K * L
    # log (product for k from 1 to K-1 of k!) = log G(K+1), Barnes G; asymptotic expansion, error O(1/K/K)
    z = D(K)
    zeta1 = D("-0.1654211437004509292139196602427806")
    two_pi = D(2) * D("3.14159265358979323846264338327950288419716939937510582097494459")
    lsf = z * z / 2 * z.ln() - D(3) / 4 * z * z + z / 2 * two_pi.ln() - z.ln() / 12 + zeta1
    logb = (D(g) * (R + S - 2) / 2).ln() - 2 * lsf / (z * z - z)
    lhs = K * (L - 1) * P
    rhs = 3 * D(N).ln() + (K - 1) * logb + D(12) / D(5) * ((R - 1) * lA1 + (S - 1) * lA2)
    marg = float((lhs - rhs) / (K * P))
    okb = (N - 1) <= D(69) * Gam + 4
    okd = (R2 - 1) * D(x).ln() / P <= D(69) * Gam
    return marg, okb, okd


def direct_checks():
    cases = []
    for (x, y) in ((10 ** 9 + 7, 2), (3, 2), (5, 2), (2, 3), (7, 3), (10 ** 12 + 39, 2), (6, 5), (101, 7)):
        P0 = C * math.log(y) * math.log(x * y) / math.log(x)
        for f in (1.0005, 1.01, 1.2, 2.0, 5.0):
            Pt = P0 * f
            if Pt > 2000:
                continue
            p = next_prime(int(decimal.Decimal(Pt).exp()))
            if (x * y) % p == 0:
                continue
            cases.append((p, x, y, f))
    worst = None
    for (p, x, y, f) in cases:
        m, okb, okd = direct_instance(p, x, y)
        assert m > 0 and okb and okd, (p, x, y, m, okb, okd)
        if worst is None or m < worst[0]:
            worst = (m, x, y, f, p.bit_length())
    say("[5] %d explicit primes (up to %d bits): condition (3) of the paper holds everywhere; minimal normalised margin %.5f (x=%d, y=%d, P/P0=%.4f, p of %d bits)"
        % (len(cases), max(c[0].bit_length() for c in cases), worst[0], worst[1], worst[2], worst[3], worst[4]))


# ---------------------------------------------------------------- [6]
def solve_C(Lv, coef_fn):
    # least constant C = k L L with k L (L-1-beta(C,k)) >= 2 coef L sqrt(k), beta at worst (y = 2, x -> infinity)
    lo, hi = 10.0, 2000.0
    for _ in range(200):
        Cm = 0.5 * (lo + hi)
        kk = Cm / (Lv * Lv)
        A = math.log(Cm) + 1.5 - LOG2 - 0.5 * math.log(kk)
        bet = 1.0 + A / (Cm * LOG2)
        ok = kk * Lv * (Lv - 1 - bet) >= 2.0 * coef_fn(Lv) * Lv * math.sqrt(kk) and Lv - 1 - bet > 0
        if ok:
            hi = Cm
        else:
            lo = Cm
    return hi


def calibration():
    cy = solve_C(5, lambda Lv: float(Lv))
    say("[6] Calibration, Yamada's coefficient L (his Theorem 1), L = 5: minimal C = %.4f (Yamada states 283)" % cy)
    assert 282.0 < cy < 283.0
    line = []
    for Lv in range(3, 10):
        c_new = solve_C(Lv, lambda t: 2.0 * ((t * t) // 4) / t)
        c_old = solve_C(Lv, lambda t: float(t))
        line.append("L=%d : %.2f (Yamada %.2f)" % (Lv, c_new, c_old))
    for s in line:
        say("    " + s)


# ---------------------------------------------------------------- [7]
def degenerate():
    v = L * math.sqrt(k)
    say("[7] degenerate case: (R2-1) log x / log p <= L sqrt(k) g a1 a2 = %.3f g a1 a2 < 69 g a1 a2 : %s"
        % (v, "OK" if v < 69 else "FAILED"))
    assert v < 69


if __name__ == "__main__":
    selftest()
    lemma_R()
    lemma_F()
    b_closed, A = beta_closed()
    b_grid, arg = beta_grid()
    say("[3] beta (closed form, y = 2, x -> infinity) = %.6f ; maximum over the grid = %.6f at (x, y) = %s : %s"
        % (b_closed, b_grid, arg, "OK" if b_grid <= b_closed + float(1e-12) else "FAILED"))
    assert b_grid <= b_closed + float(1e-12)
    final_condition(b_closed)
    direct_checks()
    calibration()
    degenerate()
    say("ALL CHECKS PASSED.")