#!/usr/bin/env python3
"""Certify the Johnston--Yang global absolute-psi decay 0.2043325."""

from __future__ import annotations

from dataclasses import dataclass
from fractions import Fraction
from hashlib import sha256
from pathlib import Path

from flint import arb, ctx


if not __debug__:
    raise SystemExit("verify.py requires assertions")


HERE = Path(__file__).resolve().parent


SOURCE_VK_DENOMINATOR = Fraction("51.34")
INTERIOR_VK_DENOMINATOR = Fraction("51.3401")
FORD_VK_DENOMINATOR = Fraction("57.54")
PUBLISHED_DECAY = Fraction("0.1853")
PRIOR_DECAY = Fraction("0.2043")
TARGET_DECAY = Fraction("0.2043325")
TARGET_AMPLITUDE = Fraction("0.026")
TARGET_LOG_POWER = Fraction("1.801")
C0067_AMPLITUDE = Fraction("9.2202181")
C0067_LOG_POWER = Fraction(3, 2)
C0067_DECAY = Fraction("0.88178")

H = 3_000_175_332_800
C1 = Fraction("17.418")
C2 = Fraction("2.077")
SIGMA = Fraction("0.9999714")
TAIL_LOG_HANDOFF = 950_000_000
TRIAL_THETA = Fraction("0.08917")
B3 = Fraction("0.21617")
PUBLISHED_TRIAL_THETA = Fraction("0.082")

SOURCE_PINS = {
    "source/johnston-yang-2204.01980v2.pdf":
        "565993a6def48b237a68a92acba604f2c42f99165e0e71e390f8e21a313b74b2",
    "source/johnston-yang-2204.01980v2.tar":
        "572d5739936ad3f5e867a142eccb0193b001dd2ee9b27b1d7183124071ec7edd",
    "source/johnston-yang-2204.01980v2.txt":
        "509a6c574dc459925b533337b5ea7e36a5a544955a3689b08b1fe8c6604c3edb",
    "source/tex/main.tex":
        "2a79d56dbd6da46d46c6ddd8852d9fa763c716110ebf08b7b029e52346f92602",
    "source/yang-thesis-2025.pdf":
        "b67b44e9608e3b73f261c4a4bd1ca70cc6000baee2aa43a4d749e44c5bd14d39",
    "source/yang-thesis-2025.txt":
        "189572ed43e57219c62499e65ef094efcb3034a02e542e1ee8a9529e5bfa51cc",
    "source/dependency/c0067/claim.md":
        "fffc5899666782461d560662b64e807ddfb4b516eae3a84a637c0c3354bf6571",
    "source/dependency/c0067/verify.py":
        "57fe65dbef4958579ba235f97bb3b657fefa7b39f91dfb25345db507d5206568",
    "source/dependency/c0070/claim.md":
        "3517deb12d395af9474c0b4fe226d72e988fea376d6ee851a0a1300afc7b9d23",
    "source/dependency/c0070/verify.py":
        "e8c9156b786b474c49449017ceeee32c75f1d16f2c39785d6625d2e2c2b1ce0d",
    "source/prior-art/tao-2503.20170v4.tar.gz":
        "68bc40d2ebf6d6a5b366ce380e1f27b0888482e88ffb6ff9599df2d15d2f4053",
    "source/prior-art/tao-2503.20170v4-notes.tex":
        "5cb98e1d90bb722abecf4428b38cf68a74abae2cd5294350dfeda3eecb8c8acc",
}

# Each pair is (prime power, underlying prime). Summing log(prime) gives psi.
PRIME_POWERS = (
    (2, 2),
    (3, 3),
    (4, 2),
    (5, 5),
    (7, 7),
    (8, 2),
    (9, 3),
    (11, 11),
    (13, 13),
    (16, 2),
    (17, 17),
    (19, 19),
    (23, 23),
    (25, 5),
    (27, 3),
    (29, 29),
    (31, 31),
    (32, 2),
    (37, 37),
    (41, 41),
    (43, 43),
    (47, 47),
    (49, 7),
    (53, 53),
    (59, 59),
)


@dataclass(frozen=True)
class Certificate:
    finite_min_margin: arb
    buthe_log_margin: arb
    buthe_derivative_margin: arb
    middle_start_margin: arb
    middle_end_margin: arb
    q_ceiling: arb
    q_formula_gap: arb
    b3_margin: arb
    explicit_lower_margin: arb
    explicit_upper_margin: arb
    native_decay: arb
    tail_s1_derivative_margin: arb
    tail_ratio: arb
    tail_margin: arb
    published_b3_margin: arb
    published_explicit_lower_margin: arb
    published_explicit_upper_margin: arb
    published_native_decay: arb
    published_s1_derivative_margin: arb
    published_repair_ratio: arb
    printed_defect: arb
    exact_improvement: Fraction
    prior_improvement: Fraction


def pin_sources() -> None:
    for relative, expected in SOURCE_PINS.items():
        path = HERE / relative
        assert path.is_file(), path
        assert sha256(path.read_bytes()).hexdigest() == expected, path

    jy_text = (
        HERE / "source/johnston-yang-2204.01980v2.txt"
    ).read_text(encoding="utf-8")
    for anchor in (
        "Theorem 1.4. For all x ≥ 23",
        "|ψ(x) − x| ≤ 0.026x(log x)1.801",
        "B2 = 0.18525 and B3 = 0.20680",
        "Letting σ = 0.9999932 we obtain (1.7)",
    ):
        assert anchor in jy_text

    jy_tex = (HERE / "source/tex/main.tex").read_text(encoding="utf-8")
    for anchor in (
        r"\begin{lemma}[Zero-density estimates]\label{zerodenlem}",
        r"if $\sigma_1<\sigma<\sigma_2$, then one can take $C_1(\sigma_2)$ "
        r"and $C_2(\sigma_1)$",
        r"0.999 & 0.3326 & 0.0545 & 0.3104 & 17.362 & 2.077\\",
        r"1 & 0.3325 & 0.0539 & 0.3105 & 17.418 & 2.069\\",
    ):
        assert anchor in jy_tex

    yang_text = (HERE / "source/yang-thesis-2025.txt").read_text(encoding="utf-8")
    for anchor in (
        "Theorem 1.3. One has ζ(σ + it) ̸= 0 if t ≥ 3",
        "51.34 (log t)2/3 (log log t)1/3",
    ):
        assert anchor in yang_text

    tao_text = (
        HERE / "source/prior-art/tao-2503.20170v4-notes.tex"
    ).read_text(encoding="utf-8")
    for anchor in (
        r"\eps_3(x) &\coloneqq 0.026 (\log^{1.801} x) "
        r"\exp(-0.1853",
        r"D. Johnston, A. Yang, \emph{Some explicit estimates for the "
        r"error term in the prime number theorem}",
    ):
        assert anchor in tao_text

    c0070_claim = (
        HERE / "source/dependency/c0070/claim.md"
    ).read_text(encoding="utf-8")
    c0070_verify = (
        HERE / "source/dependency/c0070/verify.py"
    ).read_text(encoding="utf-8")
    assert "global Vinogradov--Korobov psi decay improves to 0.2043" in (
        c0070_claim
    )
    assert 'TARGET_DECAY = Fraction("0.2043")' in c0070_verify
    assert "PASS: the global absolute-psi VK decay is 0.2043" in c0070_verify


def as_arb(value: Fraction | int) -> arb:
    value = Fraction(value)
    return arb(value.numerator) / value.denominator


def r_of_l(log_x: arb) -> arb:
    return log_x ** as_arb(Fraction(3, 5)) / (
        log_x.log() ** as_arb(Fraction(1, 5))
    )


def q_of_l(log_x: arb) -> arb:
    return as_arb(Fraction(3, 5)) - as_arb(Fraction(1, 5)) / log_x.log()


def certify_u_decreases(log_x_floor: arb) -> arb:
    """Certify u(L)=r(L)q(L)/L decreases for L>=log_x_floor."""
    ell = log_x_floor.log()
    assert ell > 1

    # With ell=log L,
    #
    # d log(u) / d ell
    # = -2/5 - 1/(5 ell) + 1/(ell(3 ell-1)).
    #
    # For ell>1, the last term is below 1/(2 ell), so the
    # derivative is at most -2/5+3/(10 ell). This upper bound
    # decreases with ell and is already negative at the floor.
    exact_at_floor = (
        -as_arb(Fraction(2, 5))
        - 1 / (5 * ell)
        + 1 / (ell * (3 * ell - 1))
    )
    upper_at_floor = -as_arb(Fraction(2, 5)) + 3 / (10 * ell)
    assert 3 * ell - 1 > 2
    assert exact_at_floor < upper_at_floor
    assert upper_at_floor < 0
    return upper_at_floor


def certify_rq_increases(log_x_floor: arb) -> arb:
    """Certify r(L)q(L) increases for L>=log_x_floor."""
    ell = log_x_floor.log()
    q = q_of_l(log_x_floor)
    assert ell > 1
    assert q > 0

    # d log(rq)/d ell = q + 1/(ell(3 ell-1)). Both summands
    # are positive, and q itself increases with ell.
    derivative_lower = q
    assert derivative_lower > 0
    return derivative_lower


def certify_middle_h_increases(log_x_floor: arb) -> arb:
    """Certify H=L^(1/10)(log L)^(-1/5)q(L) increases."""
    ell = log_x_floor.log()
    assert ell > 2

    # d log(H)/d ell
    # = 1/10 - 1/(5 ell) + 1/(ell(3 ell-1)).
    # The first two terms increase with ell and are positive at the floor.
    derivative_lower = as_arb(Fraction(1, 10)) - 1 / (5 * ell)
    assert derivative_lower > 0
    return derivative_lower


def target_value(x: arb) -> arb:
    log_x = x.log()
    return (
        as_arb(TARGET_AMPLITUDE)
        * x
        * log_x ** as_arb(TARGET_LOG_POWER)
        * (-as_arb(TARGET_DECAY) * r_of_l(log_x)).exp()
    )


def vk_ceiling(denominator: Fraction) -> tuple[arb, arb]:
    c = as_arb(denominator)
    theta_star = (
        (2 / (3 * c)) ** as_arb(Fraction(3, 5))
        * as_arb(Fraction(5, 3)) ** as_arb(Fraction(1, 5))
    )
    q_ceiling = as_arb(Fraction(5, 2)) * theta_star
    source_formula = (
        (5 / (3 * c**3)) ** as_arb(Fraction(1, 5))
        * (
            as_arb(Fraction(3, 2)) ** as_arb(Fraction(2, 5))
            + as_arb(Fraction(2, 3)) ** as_arb(Fraction(3, 5))
        )
    )
    return q_ceiling, source_formula


def finite_range_certificate() -> arb:
    log_23 = arb(23).log()
    u_23 = r_of_l(log_23) * q_of_l(log_23) / log_23

    certify_u_decreases(log_23)
    assert as_arb(TARGET_DECAY) * u_23 < 1
    target_log_derivative_lower = 1 - as_arb(TARGET_DECAY) * u_23
    assert target_log_derivative_lower > 0

    psi_value = arb(0)
    relevant: list[tuple[int, arb]] = []
    for prime_power, prime in PRIME_POWERS:
        psi_value += arb(prime).log()
        if prime_power >= 23:
            relevant.append((prime_power, psi_value))

    margins: list[arb] = []
    for index, (left, psi_left) in enumerate(relevant):
        right = relevant[index + 1][0] if index + 1 < len(relevant) else left
        lower_target = target_value(arb(left))
        margins.append(lower_target - abs(psi_left - left))
        margins.append(lower_target - abs(psi_left - right))

    for margin in margins:
        assert margin > 0
    finite_min_margin = min(margins)
    assert finite_min_margin > arb("0.12")
    return finite_min_margin


def buthe_bridge_certificate() -> tuple[arb, arb]:
    log_59 = arb(59).log()
    ell_59 = log_59.log()
    r_59 = r_of_l(log_59)
    u_59 = r_59 * q_of_l(log_59) / log_59
    certify_u_decreases(log_59)

    # Logarithm of target divided by sqrt(x) log(x)^2/(8*pi).
    log_margin = (
        (as_arb(TARGET_AMPLITUDE) * 8 * arb.pi()).log()
        - as_arb(Fraction("0.199")) * ell_59
        + log_59 / 2
        - as_arb(TARGET_DECAY) * r_59
    )
    derivative_margin = (
        as_arb(Fraction(1, 2))
        - as_arb(Fraction("0.199")) / log_59
        - as_arb(TARGET_DECAY) * u_59
    )

    assert log_margin > arb("0.8899")
    assert derivative_margin > arb("0.39")
    # Both positive subtrahends in the derivative decrease for L>=log 59,
    # so this endpoint is the global derivative minimum.
    assert as_arb(Fraction("0.199")) > 0
    assert as_arb(TARGET_DECAY) > 0
    assert arb("2.169e25").log() > 58
    return log_margin, derivative_margin


def middle_bridge_value(log_x: arb) -> arb:
    return (
        (as_arb(TARGET_AMPLITUDE) / as_arb(C0067_AMPLITUDE)).log()
        + (as_arb(TARGET_LOG_POWER) - as_arb(C0067_LOG_POWER))
        * log_x.log()
        + as_arb(C0067_DECAY) * log_x.sqrt()
        - as_arb(TARGET_DECAY) * r_of_l(log_x)
    )


def middle_bridge_certificate() -> tuple[arb, arb]:
    start = arb(58)
    end = arb(TAIL_LOG_HANDOFF)
    h_derivative_lower = certify_middle_h_increases(start)

    def h_value(log_x: arb) -> arb:
        return (
            log_x ** as_arb(Fraction(1, 10))
            * log_x.log() ** -as_arb(Fraction(1, 5))
            * q_of_l(log_x)
        )

    def b_value(log_x: arb) -> arb:
        return (
            as_arb(C0067_DECAY) / 2
            - as_arb(TARGET_DECAY) * h_value(log_x)
        )

    def g_value(log_x: arb) -> arb:
        return (
            as_arb(TARGET_LOG_POWER - C0067_LOG_POWER)
            + log_x.sqrt() * b_value(log_x)
        )

    h_start = h_value(start)
    h_end = h_value(end)
    b_start = b_value(start)
    b_end = b_value(end)
    g_start = g_value(start)
    g_end = g_value(end)

    assert h_derivative_lower > 0
    assert h_end > h_start
    assert b_start > 0
    assert b_end < 0
    assert g_start > 0
    assert g_end < 0
    assert TARGET_LOG_POWER - C0067_LOG_POWER > 0

    # F'(L)=G(L)/L with G=0.301+sqrt(L)B(L). H increases, so B
    # strictly decreases. While B>=0, G>0. Once B<0,
    #
    # G' = B/(2 sqrt(L)) + sqrt(L) B' < 0.
    #
    # Therefore G has one positive-to-negative crossing at most, and F
    # can have an interior maximum but no interior minimum.
    assert abs(
        as_arb(TARGET_LOG_POWER - C0067_LOG_POWER)
        + start.sqrt() * b_start
        - g_start
    ) < arb("1e-60")

    start_margin = middle_bridge_value(start)
    end_margin = middle_bridge_value(end)
    assert start_margin > arb("0.30")
    assert end_margin > arb("20")
    return start_margin, end_margin


def tail_terms(
    log_x: arb,
    q_ceiling: arb,
    sigma: arb,
    b3: arb,
) -> tuple[arb, arb, arb, arb, arb]:
    r = r_of_l(log_x)
    p = 5 - 2 * sigma
    a = (5 - 8 * sigma) / 3
    h_term = (arb(H) / (2 * arb.pi())).log() ** 2 / (2 * arb.pi())

    s1_low = (-log_x / 2).exp() * h_term
    s1_high = ((sigma - 1) * log_x).exp() * (
        (b3 * r) ** 2 / (2 * arb.pi()) - h_term + as_arb(Fraction("1.8642"))
    )
    s2_density = (
        2
        * as_arb(C1)
        * (a * q_ceiling * r).exp()
        * (q_ceiling * r) ** p
    )
    s2_log = (
        2
        * as_arb(C2)
        * (-q_ceiling * r).exp()
        * (q_ceiling * r) ** 2
    )
    s3 = (
        as_arb(Fraction("4.3128"))
        * log_x ** as_arb(Fraction(3, 5))
        * (-q_ceiling * r).exp()
    )
    return s1_low, s1_high, s2_density, s2_log, s3


def normalized_tail_target(log_x: arb, decay: Fraction) -> arb:
    return (
        as_arb(TARGET_AMPLITUDE)
        * log_x ** as_arb(TARGET_LOG_POWER)
        * (-as_arb(decay) * r_of_l(log_x)).exp()
    )


def b3_upper_certificate(
    log_x0: arb,
    denominator: Fraction,
    b3: Fraction,
    trial_theta_value: Fraction,
) -> arb:
    """Certify log(T)<=B3*r(L) for every L>=log_x0."""
    c = as_arb(denominator)
    ell0 = log_x0.log()
    r0 = r_of_l(log_x0)
    trial_theta = as_arb(trial_theta_value)
    trial_log_t = trial_theta * r0
    trial_coefficient = (
        trial_theta
        + (ell0 / trial_log_t.log()) ** as_arb(Fraction(1, 3))
        / (c * trial_theta ** as_arb(Fraction(2, 3)))
    )
    margin = as_arb(b3) - trial_coefficient

    assert trial_log_t > arb(H).log()
    assert margin > 0

    # For fixed theta, the nonconstant part is a positive power of
    # ell/d(ell), where
    # d(ell)=log(theta)+(3/5)ell-(1/5)log(ell).
    # Its derivative has the sign of the expression below, which is
    # negative at ell0 and strictly decreases thereafter.
    derivative_numerator = (
        trial_theta.log()
        - as_arb(Fraction(1, 5)) * ell0.log()
        + as_arb(Fraction(1, 5))
    )
    assert derivative_numerator < 0
    return margin


def tail_regime_certificate(
    log_x0: arb,
    q_ceiling: arb,
    sigma: arb,
    b3: arb,
    decay: Fraction,
) -> tuple[arb, arb, arb, arb]:
    """Certify all range and termwise monotonicity claims beyond log_x0."""
    ell0 = log_x0.log()
    r0 = r_of_l(log_x0)
    q0 = q_of_l(log_x0)
    u0 = r0 * q0 / log_x0
    p = 5 - 2 * sigma
    a = (5 - 8 * sigma) / 3
    decay_arb = as_arb(decay)

    assert ell0 > 1
    assert q0 > 0
    certify_u_decreases(log_x0)
    certify_rq_increases(log_x0)

    # The lower explicit-formula margin has derivative
    # (Q*r*q-1)/L. Since rq increases, its endpoint sign is uniform.
    explicit_lower_margin = (
        q_ceiling * r0 - (as_arb(Fraction("1.8")) * log_x0).log()
    )
    assert explicit_lower_margin > 0
    assert q_ceiling * r0 * q0 > 1

    # The upper explicit-formula margin has derivative 1/35-B3*u(L).
    # The function u decreases, so its endpoint sign is uniform.
    explicit_upper_margin = (
        log_x0 / 35
        - b3 * r0
        - as_arb(Fraction(80, 18)).log()
    )
    assert log_x0 / 35 > arb(4).log()
    assert explicit_upper_margin > 0
    assert as_arb(Fraction(1, 35)) - b3 * u0 > 0

    # Both density majorants decrease with log(T) once Qr is beyond
    # their stationary points. Since rq>0, Qr increases thereafter.
    assert a < 0
    assert q_ceiling * r0 > p / (-a)
    assert q_ceiling * r0 > 2

    native_decay = -a * q_ceiling
    assert native_decay > decay_arb
    assert p * as_arb(Fraction(3, 5)) < as_arb(TARGET_LOG_POWER)

    # The normalized x^(-1/2) term decreases because u decreases and
    # decay*u(L0)<1/2.
    assert decay_arb * u0 < as_arb(Fraction(1, 2))

    # For the second s1 term, write its positive polynomial factor as
    # alpha*r^2-k. Its logarithmic derivative is bounded by the
    # expression below. The 1/L factor and u decrease, while r and the
    # denominator increase, so the endpoint is the global maximum.
    h_term = (arb(H) / (2 * arb.pi())).log() ** 2 / (2 * arb.pi())
    k = h_term - as_arb(Fraction("1.8642"))
    alpha = b3**2 / (2 * arb.pi())
    assert k > 0
    assert alpha * r0**2 > k
    s1_high_log_derivative_upper = (
        decay_arb * u0
        + as_arb(Fraction(6, 5))
        / log_x0
        / (1 - k / (alpha * r0**2))
    )
    s1_high_derivative_margin = (
        1 - sigma - s1_high_log_derivative_upper
    )
    assert s1_high_derivative_margin > 0

    # These exact power comparisons certify the remaining normalized
    # density, logarithmic-density, and explicit-formula terms.
    assert 2 * as_arb(Fraction(3, 5)) < as_arb(TARGET_LOG_POWER)
    assert as_arb(Fraction(3, 5)) < as_arb(TARGET_LOG_POWER)
    assert q_ceiling > decay_arb

    return (
        explicit_lower_margin,
        explicit_upper_margin,
        native_decay,
        s1_high_derivative_margin,
    )


def published_repair_certificate() -> tuple[
    arb,
    arb,
    arb,
    arb,
    arb,
    arb,
    arb,
]:
    old_q, old_formula = vk_ceiling(FORD_VK_DENOMINATOR)
    assert abs(old_q - old_formula) < arb("1e-60")
    assert old_q > arb("0.19084")
    old_sigma = as_arb(Fraction("0.9999932"))

    printed_native = (
        as_arb(Fraction("0.18525")) * (8 * old_sigma - 5) / 3
    )
    printed_defect = as_arb(PUBLISHED_DECAY) - printed_native
    assert printed_defect > arb("0.000053")

    old_log_handoff = arb("2.8e10")
    old_b3_margin = b3_upper_certificate(
        old_log_handoff,
        FORD_VK_DENOMINATOR,
        Fraction("0.20680"),
        PUBLISHED_TRIAL_THETA,
    )
    (
        old_explicit_lower_margin,
        old_explicit_upper_margin,
        old_native,
        old_s1_derivative_margin,
    ) = tail_regime_certificate(
        old_log_handoff,
        old_q,
        old_sigma,
        as_arb(Fraction("0.20680")),
        PUBLISHED_DECAY,
    )
    old_terms = tail_terms(
        old_log_handoff,
        old_q,
        old_sigma,
        as_arb(Fraction("0.20680")),
    )
    old_ratio = sum(old_terms, arb(0)) / normalized_tail_target(
        old_log_handoff,
        PUBLISHED_DECAY,
    )
    assert old_ratio < arb("1e-2000")
    return (
        old_ratio,
        printed_defect,
        old_b3_margin,
        old_explicit_lower_margin,
        old_explicit_upper_margin,
        old_native,
        old_s1_derivative_margin,
    )


def tail_certificate() -> tuple[
    arb,
    arb,
    arb,
    arb,
    arb,
    arb,
    arb,
    arb,
    arb,
]:
    q_ceiling, source_formula = vk_ceiling(INTERIOR_VK_DENOMINATOR)
    q_formula_gap = abs(q_ceiling - source_formula)
    assert q_formula_gap < arb("1e-60")
    assert q_ceiling > arb("0.2043529")
    assert q_ceiling < arb("0.2043531")
    assert q_ceiling < 1

    log_x0 = arb(TAIL_LOG_HANDOFF)
    b3 = as_arb(B3)
    sigma = as_arb(SIGMA)

    assert INTERIOR_VK_DENOMINATOR > SOURCE_VK_DENOMINATOR
    assert (
        INTERIOR_VK_DENOMINATOR - SOURCE_VK_DENOMINATOR
        == Fraction(1, 10_000)
    )

    b3_margin = b3_upper_certificate(
        log_x0,
        INTERIOR_VK_DENOMINATOR,
        B3,
        TRIAL_THETA,
    )
    assert b3_margin > arb("0.000004")
    (
        explicit_lower_margin,
        explicit_upper_margin,
        native_decay,
        tail_s1_derivative_margin,
    ) = tail_regime_certificate(
        log_x0,
        q_ceiling,
        sigma,
        b3,
        TARGET_DECAY,
    )
    assert explicit_lower_margin > arb("27000")
    assert explicit_upper_margin > arb("2.71e7")
    assert native_decay > arb("0.2043373")
    assert tail_s1_derivative_margin > arb("0.000011")

    terms = tail_terms(log_x0, q_ceiling, sigma, b3)
    target = normalized_tail_target(log_x0, TARGET_DECAY)
    tail_ratio = sum(terms, arb(0)) / target
    tail_margin = 1 - tail_ratio
    assert tail_ratio < arb("0.951")
    assert tail_margin > arb("0.049")

    return (
        q_ceiling,
        q_formula_gap,
        b3_margin,
        explicit_lower_margin,
        explicit_upper_margin,
        native_decay,
        tail_s1_derivative_margin,
        tail_ratio,
        tail_margin,
    )


def certify() -> Certificate:
    ctx.prec = 256
    pin_sources()

    finite_min_margin = finite_range_certificate()
    buthe_log_margin, buthe_derivative_margin = buthe_bridge_certificate()
    middle_start_margin, middle_end_margin = middle_bridge_certificate()
    (
        q_ceiling,
        q_formula_gap,
        b3_margin,
        explicit_lower_margin,
        explicit_upper_margin,
        native_decay,
        tail_s1_derivative_margin,
        tail_ratio,
        tail_margin,
    ) = tail_certificate()
    (
        published_repair_ratio,
        printed_defect,
        published_b3_margin,
        published_explicit_lower_margin,
        published_explicit_upper_margin,
        published_native_decay,
        published_s1_derivative_margin,
    ) = published_repair_certificate()

    exact_improvement = TARGET_DECAY - PUBLISHED_DECAY
    prior_improvement = TARGET_DECAY - PRIOR_DECAY
    assert exact_improvement == Fraction(7613, 400_000)
    assert prior_improvement == Fraction(13, 400_000)

    return Certificate(
        finite_min_margin=finite_min_margin,
        buthe_log_margin=buthe_log_margin,
        buthe_derivative_margin=buthe_derivative_margin,
        middle_start_margin=middle_start_margin,
        middle_end_margin=middle_end_margin,
        q_ceiling=q_ceiling,
        q_formula_gap=q_formula_gap,
        b3_margin=b3_margin,
        explicit_lower_margin=explicit_lower_margin,
        explicit_upper_margin=explicit_upper_margin,
        native_decay=native_decay,
        tail_s1_derivative_margin=tail_s1_derivative_margin,
        tail_ratio=tail_ratio,
        tail_margin=tail_margin,
        published_b3_margin=published_b3_margin,
        published_explicit_lower_margin=published_explicit_lower_margin,
        published_explicit_upper_margin=published_explicit_upper_margin,
        published_native_decay=published_native_decay,
        published_s1_derivative_margin=published_s1_derivative_margin,
        published_repair_ratio=published_repair_ratio,
        printed_defect=printed_defect,
        exact_improvement=exact_improvement,
        prior_improvement=prior_improvement,
    )


def main() -> None:
    certificate = certify()
    print("interiorized VK denominator:", float(INTERIOR_VK_DENOMINATOR))
    print("finite 23--59 minimum margin:", certificate.finite_min_margin)
    print("Buthe log margin at x=59:", certificate.buthe_log_margin)
    print("Buthe derivative margin:", certificate.buthe_derivative_margin)
    print("C-0067 bridge margin at L=58:", certificate.middle_start_margin)
    print(
        "C-0067 bridge margin at L=9.5e8:",
        certificate.middle_end_margin,
    )
    print("uniform VK coefficient Q(c):", certificate.q_ceiling)
    print("Q formula gap:", certificate.q_formula_gap)
    print("B3 directed margin:", certificate.b3_margin)
    print("explicit-formula lower margin:", certificate.explicit_lower_margin)
    print("explicit-formula upper margin:", certificate.explicit_upper_margin)
    print("native tail decay:", certificate.native_decay)
    print(
        "tail s1 derivative margin:",
        certificate.tail_s1_derivative_margin,
    )
    print("tail coefficient ratio:", certificate.tail_ratio)
    print("tail margin:", certificate.tail_margin)
    print("published B3 margin:", certificate.published_b3_margin)
    print(
        "published explicit-formula lower margin:",
        certificate.published_explicit_lower_margin,
    )
    print(
        "published explicit-formula upper margin:",
        certificate.published_explicit_upper_margin,
    )
    print("published native decay:", certificate.published_native_decay)
    print(
        "published s1 derivative margin:",
        certificate.published_s1_derivative_margin,
    )
    print("published repaired tail ratio:", certificate.published_repair_ratio)
    print("printed B2 decay defect:", certificate.printed_defect)
    print("published decay:", float(PUBLISHED_DECAY))
    print("prior certified decay:", float(PRIOR_DECAY))
    print("improved decay:", float(TARGET_DECAY))
    print("exact decay increase:", certificate.exact_improvement)
    print("exact increase over C-0070:", certificate.prior_improvement)
    print("PASS: the global absolute-psi VK decay is 0.2043325")


if __name__ == "__main__":
    main()
