"""
classes.py -- the real support read by residue classes mod P = d - 1:
L(d) = W(d) + E(d), W the class sums' norm in closed form, E the mass
against each class's sign.

QUESTION R. With h the bounded inverse of Delta(z) = 1 - z - ... - z^d
(real_support.py), H_r = sum over m = r mod P of h_m and W = sum_r |H_r|,
is E = L - W at most A/d, A small enough that the real support closes
the carry bound past a computed d, and where in m does E sit?

BY HAND.
  (D1) At w^P = 1, w != 1, the w^j, j = 1..P, sum to 0 and w^d = w, so
       Delta(w) = 1 - w; Delta(1) = -P.
  (D2) h is in l1, so H_r = (1/P) sum_k w^(-kr)/Delta(w^k), w =
       e^(2 pi i/P); the k = 0 term is -1/P^2 and the DFT of r -> r gives
       the rest: H_r = ((P - 1)/2 - r)/P - 1/P^2, r = 0..P-1.
  (D3) W = P/4 at P even; W = (P^2 - 1)/(4P) + 1/P^2 at P odd. So W <=
       (d - 1)/4 at every d >= 3 but d = 4 (P odd needs 1/P^2 <= 1/(4P)),
       with equality at odd d; at d = 4, W = 7/9 > 3/4.
  (D4) With s_r the sign of H_r, E = sum_m e_m, e_m = |h_m| - s_r(m) h_m
       = 2 (s_r(m) h_m)^- >= 0: E is a sum over m, so its share by band
       of m is defined.
  (D5) At d = 0 mod 4 (P = 3 mod 4), C_d + 1 - W = 1/4 + 1/(4P) - 1/P^2,
       so past a computed d the bound needs E < 1/4 less the sup over t's
       excess.

PREDICTIONS AND KILL (fixed before the runs).
  R1 CONTROL: the FFT class sums equal D2's sawtooth to 1e-12 at d =
     3..100, and W equals the triangle wave's carry P/4 at d = 9, 13, 21.
  R2 d E(d) lies in [4, 7] at every d = 21..100 and falls at every step
     within each class of d mod 4. KILL of the route: d E rising over
     d = 60..100 in some class mod 4.
  R3 the share of E from |m| < 2P, printed at d = 21, 40, 60, 100; above
     1/2 at d = 100 means the edge carries E.
  R4 (fixed after R1-R3's run, before its print) at d = 40, 60, 100 more
     than 90% of E lies in the classes within 3 of r = 0 mod P or of
     (P - 1)/2, and E's shares over m/P^3 in [0, 1/16), [1/16, 1/8),
     [1/8, 1/4), [1/4, inf) agree between d = 60 and 100 to 0.02.

FINDINGS (copied from the runs' prints).
  - R1 holds: the class sums equal D2's sawtooth to 5.7e-14 and W D3's
    closed form to 1.2e-12 at d = 3..100; W = 2, 3, 5 at d = 9, 13, 21,
    the wave's P/4.
  - R2 holds: d E runs 4.6739 .. 5.6594 over d = 21..100 and falls at
    every step of 4; the kill does not fire. It splits by parity: at d =
    99, 100, E = 4.721e-02, 4.914e-02, d E = 4.6739, 4.9140 (L =
    24.547210794, 24.796717085, W = 24.5, 24.747576778).
  - R3: the edge carries at most 0.036 of E at d = 21 and none to the
    third decimal at d = 40, 60, 100. The share of E at 0 <= m < 2P is 0.032
    at d = 21 and 0.000 at d = 40, 60, 100; 2P <= m < P^3/8 holds 0.783,
    0.826, 0.828, 0.829 and m >= P^3/8 the rest, m < 0 at most 0.004.
  - R4 holds: within 3 classes of r = 0 or (P - 1)/2 lie 0.9634, 0.9661,
    0.9678 of E at d = 40, 60, 100; E's shares by m/P^3 over the four
    bands are 0.3994, 0.4297, 0.1664, 0.0046 at d = 100, within 0.0024
    of d = 60's. All three d are 0 mod 4 (P odd, 14 classes in the
    near set, its first half r = 0..3 and P - 3..P - 1 about the jump);
    d = 2 mod 4 and odd d are not printed. So at d = 0 mod 4 E reads as
    a scaling profile in m/P^3 at the sawtooth's two zeros (observation,
    three d).

HEAT FLOW (BY HAND after R1-R3, before R4; a reading, not a proof, its
root expansion a TRANSPLANT from k much below P). For m >= 0, h_m = sum_k
c_k x_k^m over the d - 1 roots inside the disc, x^d (x - 2) = -1, x_k
near e^(i theta_k), theta_k = 2 pi k/P: (d - 1) phi = 2 pi k - phi^3 +
..., |x|^d = (5 - 4 cos phi)^(-1/2), so log x_k ~ -theta_k^2/P +
i(theta_k - theta_k^3/P). At m = qP + r, x_k^m ~ e^(i theta_k r)
e^(-4 pi^2 k^2 m/P^3): in each class the leading term is the sawtooth H_r
under the heat flow for time m/P^3. The sawtooth is odd about r = -1/2
(its jump, between classes P - 1 and 0) and about r = (P - 1)/2 (its
zero), so the flow pins both zeros and adds none, and the class's sign
holds; E comes from the O(1/P) corrections (the drift theta^3/P,
c_k's corrections, the -1/P^2) moving a zero by O(1) classes: E ~ c/P,
split by parity (a class sits on the zero (P - 1)/2 only when P is odd,
d even, its H_r = -1/P^2). R4 reads it at P odd only.

RUN RECORD: the bare run 7.5 s, peak commit 198.4 MB.
"""

import os
import sys

os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import numpy as np  # noqa: E402
from real_support import inverse_fft  # noqa: E402

CHECKS = []
READS = {}


def check(name, ok, detail=""):
    CHECKS.append(bool(ok))
    print(f"  [{'PASS' if ok else 'FAIL'}] {name}" +
          (f" -- {detail}" if detail else ""), flush=True)


def closed_H(P):
    r = np.arange(P)
    return ((P - 1) / 2 - r) / P - 1 / P ** 2


def closed_W(P):
    return P / 4 if P % 2 == 0 else (P * P - 1) / (4 * P) + 1 / P ** 2


def read(d):
    """L, W (numeric), the class sums' error against D2, E and its
    pointwise split e_m by band, kept in READS."""
    if d in READS:
        return READS[d]
    P = d - 1
    hf, hb, alias = inverse_fft(d)
    m = np.concatenate([-np.arange(len(hb) - 1, 0, -1), np.arange(len(hf))])
    h = np.concatenate([hb[:0:-1], hf])
    del hf, hb
    r = m % P
    H = np.bincount(r, weights=h, minlength=P)
    err = np.abs(H - closed_H(P)).max()
    s = np.sign(H)
    e = np.abs(h) - s[r] * h
    L = np.abs(h).sum()
    W = np.abs(H).sum()
    E = e.sum()
    c = np.minimum(r, P - r)                    # distance to 0 mod P
    c = np.minimum(c, np.abs(r - (P - 1) / 2))   # or to (P - 1)/2
    near = e[c <= 3].sum()
    tau = m / P ** 3
    prof = [e[(tau >= a) & (tau < b)].sum()
            for a, b in ((0, 1 / 16), (1 / 16, 1 / 8), (1 / 8, 1 / 4),
                         (1 / 4, np.inf))]
    bands = {
        "m < 0": e[m < 0].sum(),
        "0 <= m < 2P": e[(m >= 0) & (m < 2 * P)].sum(),
        "2P <= m < P^3/8": e[(m >= 2 * P) & (m < P ** 3 / 8)].sum(),
        "m >= P^3/8": e[m >= P ** 3 / 8].sum(),
    }
    READS[d] = (L, W, err, E, bands, alias, near, prof)
    return READS[d]


def section_control():
    """R1: the class sums against D2's sawtooth and W against D3, d =
    3..100; W against the triangle wave's carry at d = 9, 13, 21."""
    print("R1/R2 THE CLASSES MOD d - 1, d = 3..100")
    errs, Wgap = [], []
    for d in range(3, 101):
        L, W, err, E, bands, alias, near, prof = read(d)
        P = d - 1
        errs.append(err)
        Wgap.append(abs(W - closed_W(P)))
        print(f"  d = {d:3d}: L {L:.9f}  W {W:.9f}  E {E:.3e}  d E "
              f"{d * E:.4f}  class err {err:.1e}  alias {alias:.1e}",
              flush=True)
    check("R1 class sums = ((P - 1)/2 - r)/P - 1/P^2 at d = 3..100",
          max(errs) < 1e-12, f"max error {max(errs):.1e}")
    check("R1 W = D3's closed form at d = 3..100", max(Wgap) < 1e-11,
          f"max gap {max(Wgap):.1e}")
    check("R1 control: W = the triangle wave's P/4 at d = 9, 13, 21",
          all(abs(read(d)[1] - (d - 1) / 4) < 1e-11 for d in (9, 13, 21)),
          " ".join(f"{read(d)[1]:.12f}" for d in (9, 13, 21)))


def section_e():
    """R2-R4: d E over d = 21..100, where in m E sits, and the classes
    and scale that carry it at d = 40, 60, 100."""
    dE = {d: d * read(d)[3] for d in range(21, 101)}
    band = [d for d in range(21, 101) if not 4 <= dE[d] <= 7]
    check("R2 d E in [4, 7] at d = 21..100", not band,
          f"outside at {band}; range {min(dE.values()):.4f}"
          f" .. {max(dE.values()):.4f}")
    rises = [d for d in range(25, 101) if dE[d] >= dE[d - 4]]
    check("R2 d E falls at every step of 4, d = 21..100", not rises,
          f"rises at {rises}")
    late = {c: [d for d in range(64, 101) if d % 4 == c and dE[d] >= dE[d - 4]]
            for c in range(4)}
    check("R2 KILL not fired: no class mod 4 rises over d = 60..100",
          not any(late.values()), f"{late}")
    print("R3 WHERE E SITS")
    for d in (21, 40, 60, 100):
        L, W, err, E, bands, alias, near, prof = read(d)
        print(f"  d = {d}: E {E:.6f}; " + "; ".join(
            f"{k} {v / E:.3f}" for k, v in bands.items()))
    print("R4 THE CLASSES AND THE SCALE THAT CARRY E")
    for d in (40, 60, 100):
        L, W, err, E, bands, alias, near, prof = read(d)
        print(f"  d = {d}: share within 3 of 0 or (P - 1)/2 {near / E:.4f};"
              f" by m/P^3 " + " ".join(f"{v / E:.4f}" for v in prof))
    check("R4 more than 90% of E within 3 classes of r = 0 or (P - 1)/2, "
          "d = 40, 60, 100", all(read(d)[6] > 0.9 * read(d)[3]
                             for d in (40, 60, 100)))
    gap = max(abs(a / read(60)[3] - b / read(100)[3])
              for a, b in zip(read(60)[7], read(100)[7]))
    check("R4 E's profile in m/P^3 agrees at d = 60 and 100 to 0.02",
          gap < 0.02, f"largest gap {gap:.4f}")


def main():
    section_control()
    section_e()
    print(f"{sum(CHECKS)}/{len(CHECKS)} checks passed")
    sys.exit(0 if all(CHECKS) else 1)


if __name__ == "__main__":
    main()