#!/usr/bin/env python3
"""Independent high-precision corroboration of the VK decay 0.2043325."""

from __future__ import annotations

import hashlib
import math
from pathlib import Path

import mpmath as mp


mp.mp.dps = 80

HERE = Path(__file__).resolve().parent
SOURCE_HASHES = {
    "source/tex/main.tex":
        "2a79d56dbd6da46d46c6ddd8852d9fa763c716110ebf08b7b029e52346f92602",
    "source/yang-thesis-2025.txt":
        "189572ed43e57219c62499e65ef094efcb3034a02e542e1ee8a9529e5bfa51cc",
    "source/dependency/c0067/claim.md":
        "fffc5899666782461d560662b64e807ddfb4b516eae3a84a637c0c3354bf6571",
    "source/dependency/c0070/claim.md":
        "3517deb12d395af9474c0b4fe226d72e988fea376d6ee851a0a1300afc7b9d23",
    "source/dependency/c0070/verify.py":
        "e8c9156b786b474c49449017ceeee32c75f1d16f2c39785d6625d2e2c2b1ce0d",
}

C = mp.mpf("51.3401")
DECAY = mp.mpf("0.2043325")
AMPLITUDE = mp.mpf("0.026")
LOG_POWER = mp.mpf("1.801")
SIGMA = mp.mpf("0.9999714")
L0 = mp.mpf("9.5e8")
TRIAL_THETA = mp.mpf("0.08917")
B3 = mp.mpf("0.21617")
H = mp.mpf("3000175332800")

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),
)


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


def r_value(log_x: mp.mpf) -> mp.mpf:
    return log_x ** mp.mpf("0.6") / mp.log(log_x) ** mp.mpf("0.2")


def q_value(log_x: mp.mpf) -> mp.mpf:
    return mp.mpf("0.6") - mp.mpf("0.2") / mp.log(log_x)


def target_absolute(x: mp.mpf) -> mp.mpf:
    log_x = mp.log(x)
    return (
        AMPLITUDE
        * x
        * log_x**LOG_POWER
        * mp.exp(-DECAY * r_value(log_x))
    )


def q_ceiling(c: mp.mpf) -> mp.mpf:
    return (
        (mp.mpf(5) / (3 * c**3)) ** mp.mpf("0.2")
        * (
            (mp.mpf(3) / 2) ** mp.mpf("0.4")
            + (mp.mpf(2) / 3) ** mp.mpf("0.6")
        )
    )


def finite_check() -> mp.mpf:
    psi = mp.mpf(0)
    intervals: list[tuple[int, mp.mpf]] = []
    for prime_power, prime in PRIME_POWERS:
        psi += mp.log(prime)
        if prime_power >= 23:
            intervals.append((prime_power, psi))

    margins = []
    for index, (left, psi_left) in enumerate(intervals):
        right = intervals[index + 1][0] if index + 1 < len(intervals) else left
        lower = target_absolute(mp.mpf(left))
        margins.append(lower - abs(psi_left - left))
        margins.append(lower - abs(psi_left - right))
    answer = min(margins)
    assert answer > mp.mpf("0.126")
    return answer


def buthe_check() -> tuple[mp.mpf, mp.mpf]:
    log_59 = mp.log(59)
    log_ratio = (
        mp.log(AMPLITUDE * 8 * mp.pi)
        - mp.mpf("0.199") * mp.log(log_59)
        + log_59 / 2
        - DECAY * r_value(log_59)
    )
    derivative = (
        mp.mpf("0.5")
        - mp.mpf("0.199") / log_59
        - DECAY * r_value(log_59) * q_value(log_59) / log_59
    )
    assert log_ratio > mp.mpf("0.8899")
    assert derivative > mp.mpf("0.40")
    assert mp.mpf(58) < mp.log(mp.mpf("2.169e25"))
    return log_ratio, derivative


def middle_log_ratio(log_x: mp.mpf) -> mp.mpf:
    return (
        mp.log(AMPLITUDE / mp.mpf("9.2202181"))
        + mp.mpf("0.301") * mp.log(log_x)
        + mp.mpf("0.88178") * mp.sqrt(log_x)
        - DECAY * r_value(log_x)
    )


def middle_check() -> tuple[mp.mpf, mp.mpf, mp.mpf]:
    start = mp.mpf(58)
    stop = L0
    values = []
    for index in range(20_001):
        weight = mp.mpf(index) / 20_000
        log_x = mp.exp(mp.log(start) * (1 - weight) + mp.log(stop) * weight)
        values.append(middle_log_ratio(log_x))
    minimum = min(values)
    assert minimum > mp.mpf("0.301")

    derivative_numerator = lambda log_x: (
        mp.mpf("0.301")
        + mp.mpf("0.88178") * mp.sqrt(log_x) / 2
        - DECAY * r_value(log_x) * q_value(log_x)
    )
    turning_point = mp.findroot(
        derivative_numerator,
        (mp.mpf("1e8"), mp.mpf("5e8")),
    )
    assert start < turning_point < stop
    assert middle_log_ratio(turning_point) > 1000
    return minimum, turning_point, middle_log_ratio(stop)


def exact_t_coefficient(log_x: mp.mpf, c: mp.mpf) -> tuple[mp.mpf, mp.mpf]:
    r = r_value(log_x)
    equation = lambda y: (
        (2 * mp.log(y) + 1)
        / (y ** (mp.mpf(5) / 3) * mp.log(y) ** (mp.mpf(4) / 3))
        - 3 * c / log_x
    )
    y = mp.findroot(equation, (mp.mpf("0.08") * r, mp.mpf("0.09") * r))
    coefficient = (
        log_x / (c * y ** (mp.mpf(2) / 3) * mp.log(y) ** (mp.mpf(1) / 3))
        + y
    ) / r
    return y / r, coefficient


def tail_components(
    log_x: mp.mpf,
    q_bound: mp.mpf,
    sigma: mp.mpf,
    b3: mp.mpf,
) -> tuple[mp.mpf, ...]:
    r = r_value(log_x)
    p = 5 - 2 * sigma
    a = (5 - 8 * sigma) / 3
    h_term = mp.log(H / (2 * mp.pi)) ** 2 / (2 * mp.pi)
    return (
        mp.exp(-log_x / 2) * h_term,
        mp.exp((sigma - 1) * log_x)
        * ((b3 * r) ** 2 / (2 * mp.pi) - h_term + mp.mpf("1.8642")),
        2 * mp.mpf("17.418") * mp.exp(a * q_bound * r) * (q_bound * r) ** p,
        2 * mp.mpf("2.077") * mp.exp(-q_bound * r) * (q_bound * r) ** 2,
        mp.mpf("4.3128") * log_x ** mp.mpf("0.6") * mp.exp(-q_bound * r),
    )


def normalized_target(log_x: mp.mpf, decay: mp.mpf) -> mp.mpf:
    return AMPLITUDE * log_x**LOG_POWER * mp.exp(-decay * r_value(log_x))


def tail_check() -> tuple[mp.mpf, mp.mpf, mp.mpf, mp.mpf]:
    q_bound = q_ceiling(C)
    theta_star = (
        (mp.mpf(2) / (3 * C)) ** mp.mpf("0.6")
        * (mp.mpf(5) / 3) ** mp.mpf("0.2")
    )
    assert abs(q_bound - mp.mpf("2.5") * theta_star) < mp.mpf("1e-70")

    theta_exact, exact_coefficient = exact_t_coefficient(L0, C)
    assert theta_exact > theta_star
    assert exact_coefficient > q_bound

    r0 = r_value(L0)
    trial_theta = TRIAL_THETA
    trial_coefficient = (
        trial_theta
        + (mp.log(L0) / mp.log(trial_theta * r0)) ** (mp.mpf(1) / 3)
        / (C * trial_theta ** (mp.mpf(2) / 3))
    )
    assert trial_coefficient < B3

    native_decay = q_bound * (8 * SIGMA - 5) / 3
    assert native_decay > DECAY
    components = tail_components(L0, q_bound, SIGMA, B3)
    ratio = sum(components) / normalized_target(L0, DECAY)
    assert ratio < mp.mpf("0.951")

    previous = ratio
    for index in range(1, 401):
        log_x = L0 * mp.power(10, mp.mpf(index) / 20)
        current = (
            sum(tail_components(log_x, q_bound, SIGMA, B3))
            / normalized_target(log_x, DECAY)
        )
        assert current < previous
        previous = current
    return q_bound, exact_coefficient, native_decay, ratio


def published_repair_check() -> tuple[mp.mpf, mp.mpf, mp.mpf, mp.mpf]:
    old_c = mp.mpf("57.54")
    old_q = q_ceiling(old_c)
    old_sigma = mp.mpf("0.9999932")
    old_l0 = mp.mpf("2.8e10")
    old_b3 = mp.mpf("0.20680")
    printed_native = mp.mpf("0.18525") * (8 * old_sigma - 5) / 3
    defect = mp.mpf("0.1853") - printed_native
    assert defect > mp.mpf("0.000053")

    r0 = r_value(old_l0)
    trial_theta = mp.mpf("0.082")
    trial_coefficient = (
        trial_theta
        + (mp.log(old_l0) / mp.log(trial_theta * r0)) ** (mp.mpf(1) / 3)
        / (old_c * trial_theta ** (mp.mpf(2) / 3))
    )
    assert trial_coefficient < old_b3

    ratio = (
        sum(tail_components(old_l0, old_q, old_sigma, old_b3))
        / normalized_target(old_l0, mp.mpf("0.1853"))
    )
    assert ratio < mp.mpf("1e-2000")

    previous = ratio
    for index in range(1, 401):
        log_x = old_l0 * mp.power(10, mp.mpf(index) / 20)
        current = (
            sum(tail_components(log_x, old_q, old_sigma, old_b3))
            / normalized_target(log_x, mp.mpf("0.1853"))
        )
        assert current < previous
        previous = current
    return defect, ratio, trial_coefficient, previous


def main() -> None:
    pin_sources()
    finite_margin = finite_check()
    buthe_margin, buthe_derivative = buthe_check()
    middle_minimum, middle_turning_point, middle_end = middle_check()
    q_bound, exact_coefficient, native_decay, tail_ratio = tail_check()
    (
        printed_defect,
        old_ratio,
        old_trial_coefficient,
        old_final_ratio,
    ) = published_repair_check()

    print("finite minimum margin:", mp.nstr(finite_margin, 30))
    print("Buthe log margin:", mp.nstr(buthe_margin, 30))
    print("Buthe derivative:", mp.nstr(buthe_derivative, 30))
    print("middle sampled minimum:", mp.nstr(middle_minimum, 30))
    print("middle turning point:", mp.nstr(middle_turning_point, 30))
    print("middle end margin:", mp.nstr(middle_end, 30))
    print("uniform Q(c):", mp.nstr(q_bound, 40))
    print("exact T coefficient at L0:", mp.nstr(exact_coefficient, 40))
    print("native decay:", mp.nstr(native_decay, 40))
    print("tail ratio:", mp.nstr(tail_ratio, 40))
    print("printed B2 defect:", mp.nstr(printed_defect, 30))
    print("published trial T coefficient:", mp.nstr(old_trial_coefficient, 30))
    print("repaired published tail ratio:", mp.nstr(old_ratio, 12))
    print("sampled published final ratio:", mp.nstr(old_final_ratio, 12))
    print("PASS: independent check confirms decay 0.2043325")


if __name__ == "__main__":
    main()
