#!/usr/bin/env python3
"""Certify the global FKS pi-minus-Li amplitude 9.221056."""

from __future__ import annotations

from dataclasses import dataclass
from fractions import Fraction

from flint import arb, ctx


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


PUBLISHED_AMPLITUDE = Fraction("9.2211")
NEW_AMPLITUDE = Fraction("9.221056")
R = Fraction("5.5666305")
A_PI = Fraction("121.107")
PUBLISHED_DECAY = Fraction("0.84768363")


@dataclass(frozen=True)
class Certificate:
    exact_amplitude: arb
    amplitude_margin: arb
    exact_decay: arb
    decay_margin: arb
    endpoint_rhs: arb
    exact_decrease: Fraction


def as_arb(value: Fraction) -> arb:
    """Convert an exact rational to an Arb enclosure."""
    return arb(value.numerator) / value.denominator


def certify() -> Certificate:
    ctx.prec = 256
    radius = as_arb(R)
    square_root_radius = radius.sqrt()

    exact_amplitude = as_arb(A_PI) / (radius * square_root_radius)
    amplitude_margin = as_arb(NEW_AMPLITUDE) - exact_amplitude
    assert exact_amplitude > arb("9.2210558")
    assert exact_amplitude < as_arb(NEW_AMPLITUDE)
    assert amplitude_margin > arb("1.18e-7")

    exact_decay = 2 / square_root_radius
    decay_margin = exact_decay - as_arb(PUBLISHED_DECAY)
    assert exact_decay > as_arb(PUBLISHED_DECAY)
    assert decay_margin > arb("3.66e-9")

    exact_decrease = PUBLISHED_AMPLITUDE - NEW_AMPLITUDE
    assert exact_decrease == Fraction(11, 250_000)
    assert NEW_AMPLITUDE < PUBLISHED_AMPLITUDE

    # FKS use Li(x)=integral_2^x dt/log(t), so Li(2)=0 and pi(2)=1.
    logarithm_two = arb(2).log()
    square_root_logarithm = logarithm_two.sqrt()
    endpoint_rhs = (
        as_arb(NEW_AMPLITUDE)
        * 2
        * square_root_logarithm
        * (-as_arb(PUBLISHED_DECAY) * square_root_logarithm).exp()
    )
    assert endpoint_rhs > 1

    return Certificate(
        exact_amplitude=exact_amplitude,
        amplitude_margin=amplitude_margin,
        exact_decay=exact_decay,
        decay_margin=decay_margin,
        endpoint_rhs=endpoint_rhs,
        exact_decrease=exact_decrease,
    )


def main() -> None:
    certificate = certify()
    print("published amplitude:", float(PUBLISHED_AMPLITUDE))
    print("improved amplitude:", float(NEW_AMPLITUDE))
    print("exact normalized amplitude:", certificate.exact_amplitude)
    print("amplitude rounding margin:", certificate.amplitude_margin)
    print("exact decay coefficient:", certificate.exact_decay)
    print("decay rounding margin:", certificate.decay_margin)
    print("exact coefficient decrease:", float(certificate.exact_decrease))
    print("x=2 right-hand side:", certificate.endpoint_rhs)
    print("PASS: 9.2211 is replaced globally by 9.221056")


if __name__ == "__main__":
    main()
