"""tailsplit.py -- which loss makes the dual tail's start J_t outgrow
r^(3/2)?

QUESTION. dualtail.py's tail passes at order s when m_s > I_r sum_m
|b_m| Lh_(s+m), Lh_n = 2n I_n (n >= 1), 2^J (n = 0); the exact
certificate passes when m_s > H_ch ||Q_s c*||_1. So the tail's test is
T = A/(l1 l2 l3): A = DUAL_s/H_ch, l1 = I_r/H_ch, l2 = sum |b_m|
L_(s+m)/||Q_s c*||_1 (the triangle inequality), l3 = sum |b_m|
Lh_(s+m)/sum |b_m| L_(s+m) (the variation lemma), L_n = ||(1 + x)^J
(1 - x)^n||_1. Read at the binding order s = r - 3.

BY HAND (the Gaussian limit at fixed r, a TRANSPLANT; closed forms by
quadrature, limits(), f = He_s phi + b2 He_(s+2) phi, b1 = 0 by
parity): l1 ->
sqrt(2 pi) 2^(r/2-1) Gamma(r/2)/(pi max|He_(r-1) phi|), 1 at odd r; l2
-> sum |b_m| int|He_(s+m) phi|/int|f|; l3 -> sum |b_m| 2n 2^(n/2-1)
Gamma(n/2)/pi over sum |b_m| int|He_n phi|/sqrt(2 pi). Their product
Lam is 4.39 .. 13.76 at ranks 5..10 against J_t/J_c = 5.11 .. 14.95; at
r = 20, 30, 40 l2 is 7.58, 9.43, 10.96 and l3 3.78, 4.71, 5.48.

PREDICTIONS (frozen before the engine). Q1 the exact l1, l2, l3 at
J = J_t(r) within 15% of their limits at ranks 5..10. Q2 l1 in
[0.95, 1.2] there. Q3 Lam(J_t)/(J_t/J_c) in [0.75, 1.25].
KILL: l1 carries the growth: l1(J_t(10))/l1(J_t(5)) > 1.5, or an l1 at
J_t outside [0.8, 1.5].
CONTROL, read first: the exact T >= 1 at J_t(r) and < 1 at 0.9 J_t(r),
ranks 5..10. THE PLANT: Lh_n = n I_n must lift T(0.9 J_t) above 1 at
some rank.

FINDINGS (copied from the printed output).
  (1) Control: T(J_t) = 1.0030, 1.0021, 1.0018, 1.0009, 1.0006, 1.0004
      and T(0.9 J_t) = 0.8975 .. 0.8983 at ranks 5..10: the exact split
      is the tail's own test, its rational bounds losing under 0.3%. The
      plant reads 1.7950 .. 1.7966, caught at every rank.
  (2) The kill misses; Q1, Q2, Q3 hold. At J_t, ranks 5..10: l1 =
      1.0000, 1.1054, 1.0043, 1.0795, 1.0024, 1.0625; l2 = 3.1758,
      4.0202, 4.0715, 4.7702, 4.8021, 5.4077; l3 = 1.5488, 1.7910,
      2.0063, 2.1986, 2.3728, 2.5336; Lam = 4.919, 7.959, 8.204,
      11.321, 11.422, 14.557 against J_t/J_c = 5.106 .. 14.953 (Lam over
      it 0.949 .. 0.975). Exact over limit: l1 and l3 within 0.4%, l2
      1.116 falling to 1.058. The plant reads nothing: n I_n halves l3
      and doubles T by arithmetic (2 x 0.8983 = 1.7966); the control's
      content is T at J_t and 0.9 J_t.
  (3) So J_t is about Lam J_c: the tail's start is the certificate's
      times the triangle inequality and the variation lemma; the
      champion's mean over its height, l1, stays within 11% of 1. In the
      limit (limits(), r = 5..10, 20, 30, 40) l2 = 2.8448 .. 10.9643 and
      l3 = 1.5431 .. 5.4799, their local slopes in log r falling from
      0.846, 0.715 (r = 5..10) to 0.532, 0.535 (r = 20..40): about
      sqrt(r) each, Lam = 61.011 at r = 40. A tail from O(J_c) needs
      ||Q_s c*||_1 bounded within a constant, not term by term.

RUN RECORD. 2.4 s, every rank exact in rationals but the Beta values
(30 digits). The limits first ran as a one-off quadrature outside the
script; limits() now prints them, the same to the printed digits; 3.2 s
in all.
"""
import math
import os
import sys
import time
from fractions import Fraction
from math import comb

os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import mpmath as mp  # noqa: E402
from bandend import Q, pmul, champ_height, gram, solve  # noqa: E402

mp.mp.dps = 30
JT = {5: 240, 6: 478, 7: 817, 8: 1267, 9: 1839, 10: 2542}  # dualtail.py
JC = {5: 47, 6: 57, 7: 97, 8: 108, 9: 157, 10: 170}        # bandend.py


def limits(r):
    """(l1, l2, l3) in the Gaussian limit at s = r - 3, by quadrature on
    a grid of step 2e-4 over |y| <= 40: f = He_s phi + b1 He_(s+1) phi +
    b2 He_(s+2) phi least in L2 (b1 = 0 by parity)."""
    import numpy as np
    from numpy.polynomial import hermite_e as He
    y = np.linspace(-40, 40, 400001)
    dy = y[1] - y[0]
    phi = np.exp(-y ** 2 / 2)

    def h(n):
        c = np.zeros(n + 1)
        c[n] = 1
        return He.hermeval(y, c) * phi
    s = r - 3
    fs = [h(s + m) for m in range(3)]
    G = np.array([[np.sum(fs[i] * fs[j]) * dy for j in range(3)]
                  for i in range(3)])
    beta = np.r_[1, np.linalg.solve(G[1:, 1:], -G[1:, 0])]
    f = sum(beta[m] * fs[m] for m in range(3))
    l1s = [np.sum(np.abs(g)) * dy for g in fs]
    l2 = (sum(abs(beta[m]) * l1s[m] for m in range(3))
          / (np.sum(np.abs(f)) * dy))
    l1 = (math.sqrt(2 * math.pi) * 2 ** (r / 2 - 1) * math.gamma(r / 2)
          / (math.pi * np.max(np.abs(h(r - 1)))))

    def ih(n):
        return 2 ** (n / 2 - 1) * math.gamma(n / 2) / math.pi
    l3 = (sum(abs(beta[m]) * 2 * (s + m) * ih(s + m) for m in range(3))
          / sum(abs(beta[m]) * l1s[m] / math.sqrt(2 * math.pi)
                for m in range(3)))
    return float(l1), float(l2), float(l3)


def I(n, J):
    """2^(J+n-1) B((J + 1)/2, n/2)/pi, the mean modulus on the circle."""
    return (mp.mpf(2) ** (J + n - 1) * mp.beta(mp.mpf(J + 1) / 2,
                                               mp.mpf(n) / 2) / mp.pi)


def split(r, J, plant=False):
    s = r - 3
    d = r - 1 - s
    y = solve(gram(J, s, d), [1] * (d + 1))
    tot = sum(y)
    c = [v / tot for v in y]
    m = 1 / tot
    b = [(-1) ** k * sum(comb(i, k) * c[i] for i in range(k, d + 1))
         for k in range(d + 1)]
    l1q = sum(abs(v) for v in pmul(Q(J, s), c))
    H = champ_height(r, J)
    L = [sum(abs(v) for v in Q(J, s + k)) for k in range(d + 1)]
    k2 = 1 if plant else 2
    Lh = [k2 * (s + k) * I(s + k, J) if s + k else mp.mpf(2) ** J
          for k in range(d + 1)]
    fr = lambda q: mp.mpf(q.numerator) / q.denominator  # noqa: E731
    A = fr(m / (l1q * H))
    l1 = I(r, J) / H
    sb = sum(abs(fr(b[k])) * L[k] for k in range(d + 1))
    l2 = sb / fr(l1q)
    l3 = sum(abs(fr(b[k])) * Lh[k] for k in range(d + 1)) / sb
    return float(A), float(l1), float(l2), float(l3), float(A / (l1 * l2 * l3))


def main():
    t0 = time.time()
    ctl, plant = True, False
    rows = {}
    for r in range(5, 11):
        at = split(r, JT[r])
        lo = split(r, int(0.9 * JT[r]))
        pl = split(r, int(0.9 * JT[r]), plant=True)
        ctl &= at[4] >= 1 and lo[4] < 1
        plant |= pl[4] >= 1
        rows[r] = at
        print("r = %2d  T(J_t) %.4f  T(0.9 J_t) %.4f  plant %.4f  [%.1f s]"
              % (r, at[4], lo[4], pl[4], time.time() - t0), flush=True)
    print("CONTROL T >= 1 at J_t and < 1 at 0.9 J_t, ranks 5..10: %s;"
          " the plant lifts T(0.9 J_t) above 1 at some rank: %s"
          % (ctl, plant), flush=True)
    LIM = {}
    for r in (5, 6, 7, 8, 9, 10, 20, 30, 40):
        LIM[r] = limits(r)
        a, b, c = LIM[r]
        print("limit r = %2d  l1 %.4f  l2 %.4f  l3 %.4f  Lam %.3f"
              % (r, a, b, c, a * b * c), flush=True)
    for r1, r2 in ((5, 10), (10, 20), (20, 40)):
        print("local slope over r = %d..%d: l2 %.3f, l3 %.3f" % (
            r1, r2, *(math.log(LIM[r2][i] / LIM[r1][i]) / math.log(r2 / r1)
                      for i in (1, 2))))
    q1 = q2 = q3 = True
    for r, (A, l1, l2, l3, T) in rows.items():
        lam = l1 * l2 * l3
        rel = [x / y for x, y in zip((l1, l2, l3), LIM[r])]
        q1 &= all(0.85 <= v <= 1.15 for v in rel)
        q2 &= 0.95 <= l1 <= 1.2
        q3 &= 0.75 <= lam / (JT[r] / JC[r]) <= 1.25
        print("r = %2d  J_t %4d  A %.3f  l1 %.4f  l2 %.4f  l3 %.4f  Lam %.3f"
              "  J_t/J_c %.3f (Lam over it %.3f)  exact/limit %.3f %.3f %.3f"
              % (r, JT[r], A, l1, l2, l3, lam, JT[r] / JC[r],
                 lam / (JT[r] / JC[r]), *rel))
    kill = rows[10][1] / rows[5][1] > 1.5 or any(
        not 0.8 <= v[1] <= 1.5 for v in rows.values())
    print("Q1 within 15%% of the limits: %s; Q2 l1 in [0.95, 1.2]: %s; Q3"
          " Lam/(J_t/J_c) in [0.75, 1.25]: %s; KILL (l1 carries the"
          " growth) %s  [%.1f s]" % (q1, q2, q3, "fires" if kill else
                                     "misses", time.time() - t0))
    return not (ctl and plant) or kill


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