#!/usr/bin/env python3
"""Independent source-table, SymPy, and mpmath reconstruction."""

from __future__ import annotations

import ast
import re
from fractions import Fraction
from math import factorial
from pathlib import Path

import mpmath as mp
import sympy as sp


ROOT = Path(__file__).resolve().parents[2]
C0050_PATH = ROOT / "canon" / "witnesses" / "C-0050" / "verify.py"
FKS_TEXT_PATH = (
    ROOT
    / "scratch"
    / "adjacent-unconditional--fks-corollary22-pi-li-amplitude-9-221056"
    / "sources"
    / "fiori-kadiri-swidinsky-2206.12557v2.txt"
)

A = Fraction(44_053, 1_000)
C = Fraction(1_673_823_191_040_000, 23)


def decimal_scientific(mantissa: str, exponent: str) -> Fraction:
    power = int(exponent)
    value = Fraction(mantissa)
    return value * (10**power if power >= 0 else Fraction(1, 10 ** (-power)))


def parse_fks_table() -> tuple[tuple[int, Fraction], ...]:
    text = FKS_TEXT_PATH.read_text()
    start = text.index("Table 4. Values for")
    end = text.index("Table 5. Sample of values", start)
    table_text = text[start:end]

    pattern = re.compile(
        r"(?<![\d.])(\d{1,3}(?: \d{3})*)\s+"
        r"(\d+\.\d+)\s+e([−-]\d+)"
    )
    parsed = []
    artifacts = []
    for log_text, mantissa, exponent in pattern.findall(table_text):
        log_x = int(log_text.replace(" ", ""))
        epsilon = decimal_scientific(mantissa, exponent.replace("−", "-"))
        if (log_x, mantissa, exponent.replace("−", "-")) in {
            (106, "3.8635", "-359"),
            (107, "1.0364", "-1153"),
        }:
            artifacts.append((log_x, epsilon))
            continue
        if 44 <= log_x <= 3_000:
            parsed.append((log_x, epsilon))

    rows = tuple(sorted(parsed))
    assert len(artifacts) == 2
    assert len(rows) == 164
    assert rows[0] == (44, Fraction(17_893, 1_000_000_000_000))
    assert rows[-1] == (3_000, Fraction(24_849, 5_000_000_000_000_000_000))
    assert len({log_x for log_x, _ in rows}) == len(rows)
    return rows


def literal_c0050_rows() -> tuple[tuple[int, Fraction], ...]:
    tree = ast.parse(C0050_PATH.read_text())
    assignment = next(
        node
        for node in tree.body
        if isinstance(node, ast.Assign)
        and any(
            isinstance(target, ast.Name) and target.id == "FKS_TABLE_ROWS"
            for target in node.targets
        )
    )

    def convert(node: ast.AST):
        if isinstance(node, ast.Tuple):
            return tuple(convert(element) for element in node.elts)
        if isinstance(node, ast.Constant) and isinstance(node.value, int):
            return node.value
        if (
            isinstance(node, ast.Call)
            and isinstance(node.func, ast.Name)
            and node.func.id == "Fraction"
        ):
            args = [convert(argument) for argument in node.args]
            return Fraction(*args)
        raise TypeError(ast.dump(node))

    return convert(assignment.value)


def normalized_majorant(log_x: int) -> Fraction:
    return sum(
        Fraction(factorial(index), log_x**index) for index in range(16)
    ) + C / log_x**16


def requirement(log_x: int, epsilon: Fraction) -> Fraction:
    reciprocal_log = Fraction(1, log_x)
    return log_x**4 * (
        1
        - reciprocal_log
        - reciprocal_log**2
        - 3 * reciprocal_log**3
        - Fraction(1, normalized_majorant(log_x) + epsilon)
    )


def symbolic_checks() -> tuple[sp.Expr, sp.Expr, sp.Expr]:
    t = sp.symbols("t")
    coefficient = sp.Rational(A.numerator, A.denominator)
    remainder = sp.Rational(C.numerator, C.denominator)
    scale = 1 - t - t**2 - 3 * t**3 - coefficient * t**4

    low_majorant = 1 + t + 2 * t**2 + 6 * t**3 + 55 * t**4
    low_residual = sp.factor(1 - scale * low_majorant)
    assert low_residual == (
        t**4
        * (
            53
            + 111_053 * t
            + 161_106 * t**2
            + 429_318 * t**3
            + 2_422_915 * t**4
        )
        / 1_000
    )

    high_majorant = sum(sp.factorial(index) * t**index for index in range(16))
    high_majorant += remainder * t**16
    monotonicity_numerator = sp.factor(
        sp.diff(1 / scale - high_majorant, t) * scale**2 / t**3
    )
    assert sp.degree(monotonicity_numerator, t) == 20
    assert monotonicity_numerator.subs(t, 0) == sp.Rational(31_053, 250)

    reserve = (coefficient - 13) / 2
    tail_factor = sp.factor(
        (1 - scale * (high_majorant + reserve * t**4)) / t**4
    )
    assert sp.degree(tail_factor, t) == 16
    assert tail_factor.subs(t, 0) == sp.Rational(31_053, 2_000)
    return low_residual, monotonicity_numerator, tail_factor


def approximate_low_requirement() -> tuple[mp.mpf, mp.mpf]:
    mp.mp.dps = 80
    left = mp.log(mp.mpf("29.53"))
    right = mp.mpf(40)
    maximum = mp.ninf
    maximum_log = None
    cells = 400_000
    for index in range(cells + 1):
        log_x = left + (right - left) * index / cells
        normalized_li = log_x * mp.e ** (-log_x) * mp.ei(log_x)
        required = log_x**4 * (
            1
            - 1 / log_x
            - 1 / log_x**2
            - 3 / log_x**3
            - 1 / normalized_li
        )
        if required > maximum:
            maximum = required
            maximum_log = log_x
    assert maximum_log is not None
    assert maximum < mp.mpf("32.51")
    return maximum, maximum_log


def main() -> None:
    parsed_rows = parse_fks_table()
    literal_rows = literal_c0050_rows()
    assert parsed_rows == literal_rows

    requirements = [
        (requirement(right, epsilon), left, right)
        for (left, epsilon), (right, _) in zip(parsed_rows, parsed_rows[1:])
    ]
    active = max(requirements)
    assert active[1:] == (2_100, 2_200)
    assert active[0] == Fraction(
        2_709_976_113_860_085_547_103_395_243_942_841_453_423_728_245_400,
        61_517_068_085_059_134_496_927_148_872_832_330_483_236_342_381,
    )

    low_residual, monotonicity_numerator, tail_factor = symbolic_checks()
    low_maximum, low_maximum_log = approximate_low_requirement()

    print("parsed FKS Table 4 rows:", len(parsed_rows))
    print("source tuple matches C-0050:", parsed_rows == literal_rows)
    print("active interval:", active[1:])
    print("active requirement exact:", active[0])
    print("active requirement decimal:", float(active[0]))
    print("44.053 minus requirement:", float(A - active[0]))
    print("low residual:", low_residual)
    print("monotonicity factor degree:", sp.degree(monotonicity_numerator))
    print("tail factor:", tail_factor)
    print("sampled low-range maximum requirement:", mp.nstr(low_maximum, 30))
    print("sampled maximizing log(x):", mp.nstr(low_maximum_log, 30))
    print("PASS")


if __name__ == "__main__":
    main()
