"""Free-chiral squashing response, with explicit error diagnostics.

Python 3.13.7, mpmath 1.3.0. No repository imports or network calls.
u=log(b), r=1/2 fixed, radius-independent real determinant representative.
P v3 (8.11) instead holds alpha=1/2 fixed; only its round Hessian is used.
"""
import argparse
import json
from functools import lru_cache

import mpmath as mp


@lru_cache(maxsize=None)
def sinhc_coefficients(dps):
    # log(sinh(t)/t), convergent for |t|<pi. Evaluation uses |t|<0.01.
    with mp.workdps(dps):
        return tuple(2**(2*n-1)*mp.bernoulli(2*n)
                     / (n*mp.factorial(2*n)) for n in range(1, 25))


def log_sinhc(t):
    if t == 0:
        return mp.mpf(0)
    if abs(t) < mp.mpf('.01'):
        return mp.fsum(c*t**(2*n) for n, c in
                       enumerate(sinhc_coefficients(mp.mp.dps), 1))
    return t + mp.log1p(-mp.exp(-2*t)) - mp.log(2*t)


def response_kernel(x):
    """Stable J=(sinh(x)cosh(x)-x)/sinh(x)^3 for real x>=0."""
    if x == 0:
        return mp.mpf(2)/3
    if x < mp.mpf('.1'):
        return (mp.expm1(log_sinhc(2*x))
                * mp.exp(-3*log_sinhc(x)) / x**2)
    e = mp.exp(-2*x)
    inverse = 2*mp.exp(-x)/(-mp.expm1(-2*x))
    return inverse*(1+e)/(1-e) - x*inverse**3


def response(cutoff=None):
    end = mp.inf if cutoff is None else mp.mpf(cutoff)
    return mp.quad(response_kernel, [0, mp.mpf('.1'), 1, 4, end])


def free_energy(u, physical=True):
    """Evaluate F before differentiating; real u, alpha strictly inside strip."""
    b, h = mp.exp(u), mp.cosh(u)
    alpha = h/2 if physical else mp.mpf('.5')

    def small(x):
        if x == 0:
            return -alpha*((2*alpha)**2-b*b-b**-2)/6
        log_ratio = (log_sinhc(2*alpha*x)-log_sinhc(b*x)
                     - log_sinhc(x/b))
        return -alpha*mp.expm1(log_ratio)/x**2

    def large(x):
        # Positive hyperbolic arguments: no complex logarithm branch.
        log_ratio = (2*alpha*x + mp.log1p(-mp.exp(-4*alpha*x))
                     -(b+b**-1)*x
                     - mp.log1p(-mp.exp(-2*b*x))
                     - mp.log1p(-mp.exp(-2*x/b)) + mp.log(2))
        return -mp.exp(log_ratio)/(2*x)

    # Integrate the subtraction alpha/x^2 exactly from 1 to infinity.
    return (mp.quad(small, [0, mp.mpf('.01'), 1])
            + mp.quad(large, [1, 4, mp.inf]) + alpha)


def literal_shell(k):
    """Differentiate each of the k unpaired mode factors in the shell."""
    a = mp.mpf('.5')
    w = mp.mpf(k)
    f_prime = 2*a/(w*w-a*a)
    f_second = -4*a*w/(w*w-a*a)**2
    return mp.fsum(f_prime*w + f_second*(2*m-k+1)**2
                   for m in range(k))


def shell_interval(K, omit_finite_part=False):
    finite = (mp.mpf(0) if omit_finite_part else -mp.mpf(1)/6)
    finite += mp.fsum(literal_shell(k)-mp.mpf(1)/3
                      for k in range(1, K+1))
    correction = 2/(mp.mpf(3)*(2*K+1))-1/(mp.mpf(2)*(2*K+1)**2)
    lo = finite+correction+1/(mp.mpf(2)*(2*K+1))
    hi = finite+correction+1/(mp.mpf(2)*(2*K-1))
    return finite, lo, hi


def run(dps=50, omit_finite_part=False):
    if dps < 40 or dps > 80:
        raise ValueError('Use 40 through 80 decimal digits for this tested series.')
    mp.mp.dps = dps
    fmt = lambda x: mp.nstr(x, 22)
    comparisons = []
    def close(actual, expected):
        error = abs(actual-expected)/max(1, abs(expected))
        assert error < mp.mpf('1e-35'), (actual, expected, error)
        comparisons.append(error)

    target = mp.pi**2/8
    integral = response()
    close(integral, target)
    F0 = free_energy(mp.mpf(0))
    close(F0, mp.log(2)/2)
    close(free_energy(mp.mpf('.1')), free_energy(mp.mpf('-.1')))

    cutoffs = []
    for L in (12, 24, 48):
        truncated = response(L)
        # Analytic tail: 0 < integral_L^infty J < integral cosh/sinh^2.
        bound = 1/mp.sinh(L)
        assert 0 < target-truncated < bound
        cutoffs.append({'L': L, 'integral_to_L': fmt(truncated),
                        'analytic_tail_upper': fmt(bound)})

    shells = []
    for K in (8, 32, 128):
        finite, lo, hi = shell_interval(K, omit_finite_part)
        selected_limit = target + (mp.mpf(1)/6 if omit_finite_part else 0)
        assert lo < selected_limit < hi
        if omit_finite_part:
            assert lo > target  # The wrong normalization is decisively detected.
        shells.append({'K': K, 'finite_sum': fmt(finite),
                       'tail_interval_lower': fmt(lo), 'tail_interval_upper': fmt(hi)})

    values = []
    for delta in map(mp.mpf, ['.1', '.05', '.025', '.0125']):
        D = 2*(free_energy(delta)-F0)/delta**2
        values.append((delta, D))
    differences = []
    for (delta, D), (_, fine) in zip(values, values[1:]):
        R = (4*fine-D)/3
        differences.append({'delta': fmt(delta), 'central': fmt(D),
                            'richardson': fmt(R),
                            'central_error': fmt(abs(D-target)),
                            'richardson_error': fmt(abs(R-target))})
    for j in range(2):
        ratio = mp.mpf(differences[j]['richardson_error'])/mp.mpf(
            differences[j+1]['richardson_error'])
        assert 15 < ratio < 17

    source_path = []
    for u in map(mp.mpf, ['.2', '.1']):
        difference = free_energy(u)-free_energy(u, physical=False)
        assert abs(difference) > mp.mpf('1e-8')
        source_path.append({'log_b': fmt(u), 'fixed_R_minus_fixed_shift': fmt(difference)})

    # Independent arithmetic/quadrature stability test. Finite steps stay fixed.
    stable_values = [integral, free_energy(mp.mpf('.025'))]
    with mp.workdps(dps+20):
        higher_values = [response(), free_energy(mp.mpf('.025'))]
    for low, high in zip(stable_values, higher_values):
        close(low, high)

    return {
        'mode': 'omitted finite part diagnostic' if omit_finite_part else 'correct subtraction',
        'working_decimal_digits': dps,
        'stability_decimal_digits': dps+20,
        'r': '1/2 fixed', 'parameter': 'u=log(b), real b>0',
        'background': 'positive-H homogeneous S3 family; zero flavor mass',
        'integral_Hessian': fmt(integral),
        'tau_RR': fmt(2*integral/mp.pi**2),
        'high_precision_comparisons': len(comparisons),
        'comparison_tolerance': '1e-35',
        'all_high_precision_comparisons_pass': True,
        'finite_integral_cutoffs': cutoffs,
        'literal_mode_shells': shells,
        'physical_finite_differences': differences,
        'finite_source_path_comparison': source_path,
        'intentional_error_detected': bool(omit_finite_part),
        'precision_note': 'Tail bounds are analytic; quadrature is precision-tested, not interval-certified. Finite differences retain truncation errors.'
    }


if __name__ == '__main__':
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--dps', type=int, default=50)
    parser.add_argument('--omit-finite-part', action='store_true')
    options = parser.parse_args()
    print(json.dumps(run(options.dps, options.omit_finite_part), indent=2))