#!/usr/bin/env python3
"""Reproduce the zero-free-region sector-cutoff arithmetic.

The conservative cutoff is justified by the elementary inequalities displayed
in the manuscript, not by rounding the high-precision candidate to an integer.
The candidate endpoint is a numerical diagnostic only until a directed-rounding
interval enclosure of log(H) and atan has been supplied.
"""
from __future__ import annotations

from fractions import Fraction
import mpmath as mp

H = 3_000_175_332_800
C = Fraction(4_896, 1_000)  # Theorem 1 of Bellotti--Trudgian--Yang (2026)
C_UPPER = Fraction(49, 10)  # 4.9
R_SAFE = 9_555_000_000_000
R_STRIP = 9_425_328_785_004


def main() -> None:
    # Exact rational checks for the coarse arithmetic after the cited theorem.
    assert H > 3_000_000_000_000
    assert C < C_UPPER
    # e > 27/10 follows from e > 1 + 1 + 1/2 + 1/6 + 1/24;
    # hence H < (27/10)^29 < e^29 proves log(H) < 29 by exact arithmetic.
    e_partial = sum((Fraction(1, __import__('math').factorial(k)) for k in range(5)), Fraction(0, 1))
    assert e_partial > Fraction(27, 10)
    assert Fraction(27, 10) ** 29 > H
    assert Fraction(1, 1) / (C_UPPER * 29) > Fraction(7, 1000)
    # theta_H < 0.986 / (3*10^12) < 3.287e-13.
    assert Fraction(986, 1000) / 3_000_000_000_000 < Fraction(3287, 10**16)
    # pi > 3.1415 implies pi/theta_H > this exact rational lower bound.
    lower_bound_A = Fraction(31415, 10000) / Fraction(3287, 10**16)
    assert lower_bound_A > 9_557_347_125_038
    assert lower_bound_A > R_SAFE + 1
    assert R_SAFE - R_STRIP == 129_671_214_996

    mp.mp.dps = 80
    Hm = mp.mpf(H)
    Cm = mp.mpf('4.896')
    logH = mp.log(Hm)
    delta = 1 / (Cm * logH)
    theta_strip = 2 * mp.atan(mp.mpf('0.5') / Hm)
    theta_zf = 2 * mp.atan((mp.mpf('0.5') - delta) / Hm)
    A_strip = mp.pi / theta_strip
    A_zf = mp.pi / theta_zf
    R_star = int(mp.floor(A_zf)) - 1

    print(f'H = {H:,}')
    print(f'C = {mp.nstr(Cm, 10)}')
    print(f'log(H) = {mp.nstr(logH, 40)}')
    print(f'delta(H) = {mp.nstr(delta, 40)}')
    print(f'theta (critical strip only) = {mp.nstr(theta_strip, 40)}')
    print(f'theta (zero-free region) = {mp.nstr(theta_zf, 40)}')
    print(f'pi/theta strip = {mp.nstr(A_strip, 50)}')
    print(f'pi/theta zero-free = {mp.nstr(A_zf, 50)}')
    print(f'R_strip = {R_STRIP:,}')
    print(f'R_safe = {R_SAFE:,}')
    print(f'conservative gain = {R_SAFE - R_STRIP:,}')
    print(f'R_star (high-precision candidate; not interval-certified) = {R_star:,}')
    print(f'candidate gain (not certified) = {R_star - R_STRIP:,}')
    print('Exact-rational checks of displayed coarse inequalities: PASS')


if __name__ == '__main__':
    main()