#!/usr/bin/env python3
"""Recompute the unique minimum of the exact third-interval formula.

Only Python's standard library is used.  Digamma and log-gamma values
are computed after recurrence to arguments >= 100, using 20 terms of
the Bernoulli expansions.  The absolute remainder is bounded by the
first omitted term for positive real arguments.  Decimal arithmetic
uses 90 significant digits; it is not directed-rounding interval
arithmetic.  The printed decimals are numerical approximations.

Run from any working directory:
    python supplement/critical_point.py
"""

from __future__ import annotations

import argparse
import json
from decimal import Decimal, localcontext
from fractions import Fraction
from pathlib import Path


D = Decimal
TERMS = 20
SHIFT_TARGET = D(100)
PRECISION = 90


def bernoulli_numbers(limit: int) -> list[Fraction]:
    """Return exact Bernoulli numbers; the sign of B_1 is immaterial here."""
    work: list[Fraction] = []
    result: list[Fraction] = []
    for m in range(limit + 1):
        work.append(Fraction(1, m + 1))
        for j in range(m, 0, -1):
            work[j - 1] = j * (work[j - 1] - work[j])
        result.append(work[0])
    return result


def as_decimal(value: Fraction) -> Decimal:
    return D(value.numerator) / D(value.denominator)


def digamma(x: Decimal, bernoulli: list[Fraction]) -> tuple[Decimal, Decimal]:
    y, correction = x, D(0)
    while y < SHIFT_TARGET:
        correction -= 1 / y
        y += 1
    value = correction + y.ln() - 1 / (2 * y)
    for n in range(1, TERMS + 1):
        value -= as_decimal(bernoulli[2 * n]) / (D(2 * n) * y ** (2 * n))
    order = 2 * TERMS + 2
    remainder_bound = abs(as_decimal(bernoulli[order])) / (D(order) * y**order)
    return value, remainder_bound


def critical_equation(
    alpha: Decimal, bernoulli: list[Fraction]
) -> tuple[Decimal, Decimal]:
    p, error_p = digamma(2 * alpha + 1, bernoulli)
    q, error_q = digamma(alpha + 1, bernoulli)
    value = 2 * p - 2 * q - 1 / ((alpha + 1) * (alpha + 2))
    return value, 2 * error_p + 2 * error_q


def arctangent_small(z: Decimal) -> Decimal:
    square = z * z
    term, total = z, z
    j = 1
    while abs(term) > D("1e-88"):
        term *= -square
        total += term / D(2 * j + 1)
        j += 1
    return total


def logarithmic_gamma(
    x: Decimal, pi: Decimal, bernoulli: list[Fraction]
) -> tuple[Decimal, Decimal]:
    y, correction = x, D(0)
    while y < SHIFT_TARGET:
        correction -= y.ln()
        y += 1
    value = correction + (y - D("0.5")) * y.ln() - y + (2 * pi).ln() / 2
    for n in range(1, TERMS + 1):
        value += as_decimal(bernoulli[2 * n]) / (
            D(2 * n) * D(2 * n - 1) * y ** (2 * n - 1)
        )
    order = 2 * TERMS + 2
    remainder_bound = abs(as_decimal(bernoulli[order])) / (
        D(order) * D(order - 1) * y ** (order - 1)
    )
    return value, remainder_bound


def compute() -> dict[str, object]:
    with localcontext() as context:
        context.prec = PRECISION
        bernoulli = bernoulli_numbers(2 * TERMS + 2)
        lo, hi = D(0), D(1)
        # This computes an approximation; the separate fixed bracket below
        # has a much larger sign margin than the analytic remainder bounds.
        for _ in range(200):
            mid = (lo + hi) / 2
            value, _ = critical_equation(mid, bernoulli)
            if value > 0:
                hi = mid
            else:
                lo = mid
        alpha = (lo + hi) / 2
        bracket = (
            D("0.1646512224777764061796936438236800488543"),
            D("0.1646512224777764061796936438236800488544"),
        )
        endpoint_data = []
        for endpoint in bracket:
            value, truncation_error = critical_equation(endpoint, bernoulli)
            endpoint_data.append(
                {
                    "alpha": str(endpoint),
                    "h_value": str(value),
                    "analytic_truncation_error_bound": str(truncation_error),
                }
            )
        left_value, left_error = critical_equation(bracket[0], bernoulli)
        right_value, right_error = critical_equation(bracket[1], bernoulli)
        assert left_value + left_error < 0 < right_value - right_error
        assert bracket[0] < alpha < bracket[1]

        pi = 16 * arctangent_small(D(1) / 5) - 4 * arctangent_small(D(1) / 239)
        numerator, numerator_error = logarithmic_gamma(2 * alpha + 1, pi, bernoulli)
        denominator, denominator_error = logarithmic_gamma(alpha + 1, pi, bernoulli)
        c0 = (numerator - 2 * denominator).exp()
        minimum = 4 * pi * (alpha + 2) / (alpha + 1) * c0
        return {
            "parameter": "lambda = 2 + alpha",
            "alpha_star": format(alpha, ".55f"),
            "lambda_star": format(2 + alpha, ".55f"),
            "F_at_minimum": format(minimum, ".55f"),
            "rational_bracket": [str(x) for x in bracket],
            "endpoint_evaluations": endpoint_data,
            "precision_digits": PRECISION,
            "bernoulli_terms": TERMS,
            "shift_target": str(SHIFT_TARGET),
            "bisection_iterations": 200,
            "bisection_width": str(hi - lo),
            "log_C_analytic_truncation_error_bound": str(
                numerator_error + 2 * denominator_error
            ),
            "arithmetic_status": (
                "High-precision Decimal arithmetic with analytic truncation bounds; "
                "not directed-rounding interval arithmetic."
            ),
            "mathematical_basis": (
                "The critical equation is strictly increasing on (0,1), "
                "with endpoint values -1/2 and 5/6, so its zero is unique."
            ),
        }


def main() -> None:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument(
        "--output",
        type=Path,
        default=Path(__file__).resolve().parent / "results" / "critical_point.json",
    )
    args = parser.parse_args()
    result = compute()
    args.output.parent.mkdir(parents=True, exist_ok=True)
    args.output.write_text(json.dumps(result, indent=2) + "\n", encoding="utf-8")
    print(f"alpha_star = {result['alpha_star']}")
    print(f"F_at_minimum = {result['F_at_minimum']}")
    print(f"Wrote {args.output}")


if __name__ == "__main__":
    main()