"""The doubler term of the one-loop anisotropy and its logarithm: Proposition 6.5 and the rate of Theorem 6.8.

Paper: "Linearized Einstein gravity from Lorentzian Regge calculus on the Kuhn lattice, and the divergence of its one-loop
graviton limit", J. Washburn (2026), https://doi.org/10.5281/zenodo.23764491. Section, theorem and proposition numbers below
are the paper's.

What this program does
  * Proposition 6.5 reduces the model double pole at each axis corner to one integral,
        I(eps) = -pi^2 int_0^inf exp(-eps tau) tau (1 + i tau) (1 + tau^2)^(-3/2) d tau.
    The program evaluates its real and imaginary parts in 30-digit quadrature at eps = 10^(-n/4), n = 0, ..., 48, and checks that
    Re I(eps) tends to -pi^2 and that Im I(eps) = -pi^2 J(eps) with J(eps) = log(1/eps) + log 2 - 1 - gamma + o(1), gamma
    Euler's constant (the integrand of J falls like 1/tau until the regulator cuts it off near tau = 1/eps).
  * Section 6.4 writes the doubler term as D(eps) = rho Im I(eps) / (2 |det L| (2 pi)^4) with rho = -3/7688 and
    |det L| = 3 sqrt(1013375)/1240^2. Since 1013375 = 55^2 * 335, the program checks in exact rational arithmetic that this is
    D(eps) = (5/(44 sqrt335 pi^2)) J(eps): the doubler term grows by 5/(44 sqrt335 pi^2), about 6.3e-4, per e-fold of 1/eps,
    the rate of Theorem 6.8.
  * Proposition 6.13 gives each corner line a constant dressing z, replacing (r + i eps)^(-1) by (r + z + i eps)^(-1). Then
        I(eps; z) = -pi^2 int_0^inf exp((i z - eps) tau) tau (1 + i tau) (1 + tau^2)^(-3/2) d tau,
    and the dressed doubler term is D_z(eps) = (5/(44 sqrt335 pi^2)) J_z(eps) with J_z = -Im I(eps; z)/pi^2. The program
    computes J_z by quadrature up to tau = 200 and the exact large-tau series beyond (each term an exponential integral E_n),
    and checks: a width z = i mu is the regulator shifted, J_{i mu}(eps) = J(eps + mu); for a mass (z > 0), a tachyonic shift
    (z < 0) and a width alike, J_z(eps) stops changing once eps is far below |z|, while J(eps) keeps growing; and the value it
    stops at lies within a few |z| log(1/|z|) of J(|z|) (numerical).
  * With --json FILE.json, writes J(eps) and D(eps) at the 49 regulators: the data of the figure "The doubler term D(eps)" on
    the gravity page. With --dressed as well, also writes D_z(eps) for a mass z = 1e-3, 1e-6 and 1e-9: the figure's thin curves.

Requires Python 3.9+ and mpmath. Run: python3 doubler_log.py   [--json log.json [--dressed]]
License: MIT. Copyright (c) 2026 Recognition Physics Institute.
"""
import argparse
import json
import sys
from fractions import Fraction

from mpmath import mp, mpf, mpc, quad, exp, expint, log, sqrt, pi, euler

mp.dps = 30
NSTEPS = 48                      # eps = 10^(-n/4), n = 0..48
DRESSINGS = ('1e-3', '1e-6', '1e-9')   # the masses z of the figure's thin curves
TAIL_T = mpf(200)                # J_z: quadrature on [0, T], exact series beyond
TAIL_TERMS = 16


def J(eps):
    """int_0^inf exp(-eps tau) tau^2 (1 + tau^2)^(-3/2) d tau, in the variable s = log tau (integrand tau^3 (1 + tau^2)^(-3/2) ...)."""
    L = log(1 / eps)
    f = lambda s: exp(-eps * exp(s) + 3 * s) / (1 + exp(2 * s)) ** mpf(1.5)
    return quad(f, [-40, -5, 0, 5, L - 3, L, L + 3, L + 6])


def re_part(eps):
    """-pi^2 int_0^inf exp(-eps tau) tau (1 + tau^2)^(-3/2) d tau, the real part of I(eps)."""
    L = log(1 / eps)
    f = lambda s: exp(-eps * exp(s) + 2 * s) / (1 + exp(2 * s)) ** mpf(1.5)
    return -pi ** 2 * quad(f, [-40, -5, 0, 5, L - 3, L, L + 3, L + 6])


def Jz(eps, z):
    """J_z(eps) = Im int_0^inf tau (1 + i tau) exp((i z - eps) tau) (1 + tau^2)^(-3/2) d tau, Proposition 6.13's dressed integral.

    On [0, T] by quadrature (subdivided by the oscillation when z has a real part); beyond T the integrand is
    exp(w tau) g(tau) with w = i z - eps and g(tau) = sum_n c_n tau^(-n), the series of (u^2 + i u)(1 + u^2)^(-3/2) in u = 1/tau,
    and int_T^inf exp(w tau) tau^(-n) d tau = T^(1-n) E_n(-w T) exactly."""
    z = mpc(z)
    w = 1j * z - eps
    g = lambda t: t * (1 + 1j * t) / (1 + t * t) ** mpf(1.5)
    pts = [mpf(0), mpf(1), mpf(10), mpf(50), TAIL_T]
    if abs(z.real) * TAIL_T > 2:
        n = int(abs(z.real) * TAIL_T / pi) + 1
        pts = sorted(set(pts + [TAIL_T * k / n for k in range(1, n)]))
    head = quad(lambda t: exp(w * t) * g(t), pts)
    tail = sum(c * TAIL_T ** (1 - n) * expint(n, -w * TAIL_T) for n, c in enumerate(_SERIES) if c != 0)
    return (head + tail).imag


def tail_series(nterms):
    """c_n with g(tau) = sum_n c_n tau^(-n) for tau > 1: g = (u^2 + i u)(1 + u^2)^(-3/2) with u = 1/tau, and
    (1 + u^2)^(-3/2) = sum_k binom(-3/2, k) u^(2k), so c_(2k+1) = i binom(-3/2, k) and c_(2k+2) = binom(-3/2, k)."""
    c = [mpc(0)] * (nterms + 1)
    b = mpf(1)
    for k in range(nterms):
        if 2 * k + 1 <= nterms:
            c[2 * k + 1] += 1j * b
        if 2 * k + 2 <= nterms:
            c[2 * k + 2] += b
        b *= (mpf(-3) / 2 - k) / (k + 1)
    return c


_SERIES = tail_series(TAIL_TERMS)


def dressed_curves(ns, rate):
    """D_z(eps) = rate J_z(eps) at eps = 10^(-n/4) for each mass z in DRESSINGS: the figure's thin curves."""
    return {zs: [rate * Jz(mpf(10) ** (mpf(-n) / 4), mpf(zs)) for n in ns] for zs in DRESSINGS}


def report(name, ok, detail):
    print(f"  [{'PASS' if ok else 'FAIL'}] {name}: {detail}")
    return ok


def main():
    parser = argparse.ArgumentParser(description='The doubler term D(eps) of Section 6.4 and its logarithm.')
    parser.add_argument('--json', metavar='FILE.json', help='write J(eps) and D(eps) at eps = 10^(-n/4), n = 0..48')
    parser.add_argument('--dressed', action='store_true', help='with --json, also write D_z(eps) for a mass z = 1e-3, 1e-6, 1e-9')
    args = parser.parse_args()
    ok = True

    print('Section 6.4: the prefactor of the doubler term')
    # |rho| pi^2 / (2 |det L| (2 pi)^4) = (3/7688) / (2 * 3 * 55 sqrt335 / 1240^2 * 16) / pi^2 = q / (sqrt335 pi^2)
    q = Fraction(3, 7688) / (2 * Fraction(3 * 55, 1240 ** 2) * 16)
    ok &= report('rho Im I / (2 |det L| (2 pi)^4) = (5/(44 sqrt335 pi^2)) J, exactly', q == Fraction(5, 44) and 55 ** 2 * 335 == 1013375,
                 f'(3/7688) * 1240^2 / (2 * 3 * 55 * 16) = {q}, and 1013375 = 55^2 * 335')
    rate = mpf(5) / (44 * sqrt(335) * pi ** 2)
    print(f'  [INFO] rate 5/(44 sqrt335 pi^2) = {mp.nstr(rate, 12)} per e-fold, {mp.nstr(rate * log(10), 12)} per factor of ten')

    print('Proposition 6.5: the integral I(eps)')
    ns = list(range(NSTEPS + 1))
    Js = [J(mpf(10) ** (mpf(-n) / 4)) for n in ns]
    const = log(2) - 1 - euler
    tail = [Js[n] - (n / mpf(4)) * log(10) - const for n in (36, 40, 44, 48)]
    ok &= report('J(eps) - log(1/eps) -> log 2 - 1 - gamma', abs(tail[-1]) < mpf('1e-10') and all(abs(a) >= abs(b) for a, b in zip(tail, tail[1:])),
                 f'J - log(1/eps) - (log 2 - 1 - gamma) = {", ".join(mp.nstr(t, 3) for t in tail)} at eps = 1e-9, 1e-10, 1e-11, 1e-12')
    per_decade = [Js[n + 4] - Js[n] for n in (36, 40, 44)]
    ok &= report('each factor of ten in eps adds log 10 to J', all(abs(d - log(10)) < mpf('1e-8') for d in per_decade),
                 f'increments {", ".join(mp.nstr(d, 12) for d in per_decade)}; log 10 = {mp.nstr(log(10), 12)}')
    reals = [re_part(mpf(10) ** (mpf(-n) / 4)) for n in (24, 36, 48)]
    ok &= report('Re I(eps) -> -pi^2', abs(reals[-1] + pi ** 2) < mpf('1e-9'),
                 f'Re I = {", ".join(mp.nstr(r, 12) for r in reals)} at eps = 1e-6, 1e-9, 1e-12; -pi^2 = {mp.nstr(-pi ** 2, 12)}')
    Ds = [rate * j for j in Js]
    print(f'  [INFO] D(eps): {mp.nstr(Ds[24], 6)} at eps = 1e-6, {mp.nstr(Ds[48], 6)} at eps = 1e-12')

    print('Proposition 6.13: a constant dressing z of the corner lines')
    same = [abs(Jz(mpf(10) ** (mpf(-n) / 4), 0) - Js[n]) for n in (0, 24, 48)]
    ok &= report('the dressed integral at z = 0 is J(eps)', max(same) < mpf('1e-20'),
                 f'largest difference {mp.nstr(max(same), 3)} at eps = 1, 1e-6, 1e-12')
    width = [abs(Jz(mpf(e), mpc(0, mu)) - J(mpf(e) + mpf(mu))) for e, mu in (('1e-3', '1e-4'), ('1e-8', '1e-6'))]
    ok &= report('a width z = i mu shifts the regulator: J_{i mu}(eps) = J(eps + mu)', max(width) < mpf('1e-20'),
                 f'differences {", ".join(mp.nstr(d, 3) for d in width)}')
    for label, z, zt in (('a mass', mpf('1e-3'), '1e-3'), ('a tachyonic shift', -mpf('1e-3'), '-1e-3'), ('a width', mpc(0, '1e-3'), '1e-3 i')):
        vals = [Jz(mpf(10) ** -k, z) for k in (9, 12)]
        ok &= report(f'{label}, z = {zt}: J_z has stopped growing below eps = 1e-9', abs(vals[1] - vals[0]) < mpf('1e-5'),
                     f'J_z = {mp.nstr(vals[0], 12)} at 1e-9, {mp.nstr(vals[1], 12)} at 1e-12 (J grows by 3 log 10 = {mp.nstr(3 * log(10), 6)})')
    near = []
    for zs in ('1e-3', '1e-6'):
        z = mpf(zs)
        lim = [Jz(mpf('1e-12'), s) for s in (z, -z, mpc(0, z))]
        near.append(max(abs(v - J(z)) for v in lim) / (z * log(1 / z)))
        print(f'  [INFO] |z| = {zs}: J_z(1e-12) = {", ".join(mp.nstr(v, 10) for v in lim)} for a mass, a tachyonic shift and a width;'
              f' J(|z|) = {mp.nstr(J(z), 10)}')
    ok &= report('each stops within a few |z| log(1/|z|) of J(|z|)', max(near) < 3,
                 f'largest |J_z - J(|z|)| / (|z| log(1/|z|)) = {mp.nstr(max(near), 3)} (numerical)')
    dressed = dressed_curves(ns, rate) if args.dressed else {}
    if dressed:
        print(f'  [INFO] D_z(1e-12) for z = {", ".join(DRESSINGS)}: {", ".join(mp.nstr(dressed[zs][-1], 6) for zs in DRESSINGS)}')

    if args.json:
        with open(args.json, 'w') as fh:
            out = {'eps_exponent_quarter_steps': ns, 'J': [f'{float(j):.9f}' for j in Js], 'D': [f'{float(d):.6e}' for d in Ds],
                   'rate': float(rate), 'const': float(const)}
            if dressed:
                out['D_dressed'] = {zs: [f'{float(d):.6e}' for d in v] for zs, v in dressed.items()}
            json.dump(out, fh)
        print(f'  wrote {args.json}')
    print('ALL CHECKS PASS' if ok else 'SOME CHECKS FAILED')
    return ok


if __name__ == '__main__':
    sys.exit(0 if main() else 1)