"""quadtop.py -- the place with residue field F_2 tops the thermal place
spectrum over every quadratic field that has one.

QUESTION. fieldtop.py finds the F_2 place losing the top at three cyclic
octics and one cyclic quartic, and topping at every quadratic field it
reads (least gap 0.030 at D = 264). A quadratic field holds at most two
places of each norm. Does the F_2 place top over every quadratic field,
or is there one where it loses?

THE ARGUMENT. K quadratic with a place of residue field F_2, so 2 splits
or ramifies; x = 2^-b.
  (1) THE COMPARISON (fieldtop.py's (1)). An odd column P sits below
      the F_2 place iff its headroom sum sigma_P exceeds the F_2
      column's sigma_2 at b = beta_F2.
  (2) THE F_2 COLUMN. sigma_2 = T2 prod_R (1 + N(R)^-b), T2 =
      (1 - x)^-(g - 1), g the number of F_2 places (2 at a split 2, 1
      at a ramified one), R over the places of Fermat norm 2^k + 1. A
      norm in K is p or p^2, and p^2 - 1 = 2^k only at p = 3, so these
      norms are 3, 5, 9, 17, 257, 65537 and the Fermat primes past
      65537. Those are at least F_5 = 2^32 + 1, and the sum of F_m^-b
      over m >= 5 is below 2 F_5^-b, so their factor is below
      exp(2 n F_5^-b), n = 2.
  (3) AN ODD COLUMN FROM BELOW. P of norm N, v = v_2(N - 1). Its
      headroom holds (a) each F_2 place to depth v + 2 where 2 is
      unramified (the chain of (Z/2^j)^*), and where 2 ramifies
      (e = 2) to depth 2 at v = 1 and at least 4 at v >= 2
      (fieldtop.py's wild chain (3)); (b) every place R != P of
      Fermat norm 2^k + 1 with k <= v at depth >= 1, since
      lambda(R) = 2^k divides N - 1; (c) where N is a Fermat norm, the
      conjugate of P without bound. Every factor is at least 1, so
      dropping all other places lowers sigma_P, and sigma_P / sigma_2
      >= A X / (T2 W): A = prod over the F_2 places of
      (1 + x + ... + x^depth), X = 1/(1 - y^2) for a conjugate of
      Fermat norm (y = N^-b; sigma_2 holds it at depth 1), W = P's own
      factor 1 + y where N is a Fermat norm, times the places of Fermat
      norm with k > v, times the bound of (2) on the Fermat primes past
      65537.
  (4) THE FLOOR AND THE CEILING. zeta_K >= its local factors at 2, 3,
      5, 17, 257, 65537 times (1 - p^-2b)^-1 at every other p (the
      inert factor, the least of the three), and Psi_2 = zeta_K -
      sigma_2 is at least that less the bound on sigma_2. Psi_2
      decreases, so every b where this floor is >= 1 lies at or below
      beta_F2. And Psi_2 <= zeta_K - 1 <= zeta^2 - 1, so beta_F2 <
      b_top, the least b with zeta(b)^2 < 2.
  (5) THE CLASSES. A X and T2 W read K only through the types of 2
      (split, ramified) and of 3, 5, 17, 257, 65537 (split, ramified,
      inert): 486 configurations. They read the column only through v
      and whether N is a Fermat norm. W moves with v only where v
      crosses k in {1, 2, 3, 4, 8, 16}, and A grows with v, so the
      generic classes v = 1, 2, 3, 4, 8, 16 bound the rest (v = 5..7
      by 4, 9..15 by 8, past 16 by 16). A column of Fermat norm 2^k +
      1 is its own class at v = k; one over a Fermat prime past 65537
      is bounded by v = 16, whose tail bound counts its own place and
      its conjugate.
  (6) THE BRACKET. A X and T2 W are products of terms that decrease in
      b, so A X(b_(i+1)) > T2 W(b_i) at consecutive points of a grid
      covering [b_low, b_top] leaves no b there with A X <= T2 W. Where
      it holds at every class of every configuration, every odd column
      sits below the F_2 place, by (1).

DESIGN. One grid, step 1e-4. b_low: the largest grid point of
[1.0001, 4] where the floor of (4) is >= 1, read by a 0.01 scan refined
in its last bracket; any grid point where the floor holds is a floor of
beta_F2, so the scan's exactness is not needed. b_top: the grid point
past the last one where zeta^2 >= 2 (zeta decreases, so the scan is
exact). R, the bracket of (6) on a grid from b_low past b_top, at every
class of every configuration, printing the count of failing columns
and the least A X(b_(i+1)) / T2 W(b_i) per (2's type, class). At the least
column, b_c: the largest point of the grid on [1.01, 12] where
A X <= T2 W, beside b_low (a print, not the proof). Every comparison is
in double precision, the least ratio far above its rounding.
CONTROLS. (C) the engine at a field where the F_2 place loses, the
cyclic quartic of conductor 272 (fieldtop.py's O): one F_2 place,
e = 4; four places of norm 5, two of norm 9, one of 17, four of norm 47.
The column of norm 47 has v = 1, holds the place over 2 to depth 2 and
its three conjugates; its A X < T2 W at that field's F_2 point
1.54337, equal to the product written out by hand, and the floor sits
below that point. (C2) the floor of each of the twelve quadratic fields
fieldtop.py reads sits at or below its printed F_2 point.

PREDICTIONS, fixed before the first run. (h1) The least column is the
one of norm 3 where 2 and 3 ramify and 5, 17, 257, 65537 split, with
b_low - b_c negative, between -0.03 and 0. (h2) Every configuration
with 2 split passes. (h3) Every column with v >= 2 passes.

FINDINGS. C: A X = 1.472359 against T2 W = 1.490681 at 1.54337, each
equal to the hand product; floor 1.4793. C2: floors 1.4458 to 1.5659,
each below its field's F_2 point (D = -84, at the least configuration,
1.4641 under 1.65813). beta_F2 < 2.3353. R: 0 failing columns of 486
configurations. The least ratio per (2's type, class): the column of
norm 3 at 1.005847 (2 ramified) and 1.062621 (2 split), each where 3
ramifies and 5, 17, 257, 65537 split; generic v = 1 at 1.129157 and
1.177085, where 3 is inert; every other class 1.2075 or more. At the
least column b_low 1.4641, b_c 1.4424, b_low - b_c = +0.0217. (h1) holds
as to the column and misses as to the sign: the margin is positive.
(h2) holds. (h3) holds.
READING. Over every quadratic field with a place of residue field F_2,
that place's point is the largest of the spectrum. The race is closest
at the column over 3 where 2 and 3 ramify: there the bound counts
nothing but the F_2 place to depth 2, against the F_2 column's two
places each of norm 5, 17, 257, 65537 and the one of norm 3. A real
such field adds the places of norm 2 3^i + 1 (7, 19, 163, 487, ...)
over the primes of that form that split or ramify.

RUN RECORD. `python quadtop.py`: 3 of 3 checks, 23.6 s, peak commit
47.6 MB. With the configurations counted: 3 of 3, 24.1 s, 47.2 MB,
0 failing columns of 486.
"""

import os
import sys
from itertools import product

os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))

import numpy as np  # noqa: E402

from fieldtop import kronecker  # noqa: E402
from spectrum import zeta  # noqa: E402

REL = (3, 5, 17, 257, 65537)
FERMAT_NORMS = {3: 1, 5: 2, 9: 3, 17: 4, 257: 8, 65537: 16}   # N: k
F5 = 2.0 ** 32 + 1             # the least Fermat number past 65537
STEP = 1e-4
CHECKS = []


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


def grid(lo, hi, step=STEP):
    n = int(round((hi - lo) / step))
    return [lo + i * step for i in range(n + 1)]


def last_true(pred, lo, hi):
    """The largest b on the 1e-4 grid of [lo, hi] with pred(b), read by
    a 0.01 scan and a 1e-4 refinement of the last bracket that holds a
    True (exact unless pred flips twice inside one 0.01 cell); None if
    no point is found. Whatever it returns satisfies pred."""
    hits = [b for b in grid(lo, hi, 0.01) if pred(b)]
    b0 = hits[-1] if hits else lo
    best = None
    for b in grid(b0, min(hi, b0 + 0.01)):
        if pred(b):
            best = b
    return best


class Profile:
    """A field seen at its relevant primes: degree n, the places over 2
    as (e, f, count), and per relevant odd p a list of (norm, count)."""

    def __init__(self, n, two, odd):
        self.n, self.two, self.odd = n, two, odd

    def g2(self):
        return sum(c for e, f, c in self.two if f == 1)

    def zeta_floor(self, b):
        """zeta_K >= this: the profile's local factors at 2 and REL, the
        least factor (1 - p^-nb)^-1 at every other prime."""
        z = zeta(self.n * b) * (1 - 2 ** (-self.n * b))
        for e, f, c in self.two:
            z *= (1 - 2 ** (-f * b)) ** -c
        for p, pl in self.odd.items():
            z *= 1 - p ** (-self.n * b)
            for N, c in pl:
                z *= (1 - N ** -b) ** -c
        return z

    def fermat_places(self):
        """[(N, k, count)] of the odd places of Fermat norm."""
        out = []
        for p, pl in self.odd.items():
            for N, c in pl:
                if N in FERMAT_NORMS:
                    out.append((N, FERMAT_NORMS[N], c))
        return out

    def tail(self, b):
        """The Fermat primes past 65537, at most n places each: the sum
        of F_m^-b over m >= 5 is below 2 F_5^-b."""
        return np.exp(self.n * 2 * F5 ** -b)

    def sigma2(self, b):
        x = 2 ** -b
        s = (1 - x) ** -(self.g2() - 1)
        for N, k, c in self.fermat_places():
            s *= (1 + N ** -b) ** c
        return s * self.tail(b)

    def b_low(self):
        return last_true(lambda b: self.zeta_floor(b) - self.sigma2(b)
                         >= 1, 1.0001, 4.0)


def depth2(e, v):
    """The depth to which an odd column with v = v_2(N - 1) holds an
    F_2 place of ramification e (f = 1); a lower bound past 2e."""
    if e == 1:
        return v + 2
    return min(2 ** v, 2 * e)


def column_ratio(P, b, v, own=None, conj=0):
    """(A X, T2 W) at b for an odd column with v; own = its norm when
    it is a Fermat norm (then 2^v + 1), conj its conjugates of that
    norm, held without bound."""
    x = 2 ** -b
    A = 1.0
    for e, f, c in P.two:
        if f == 1:
            A *= sum(x ** j for j in range(depth2(e, v) + 1)) ** c
    X = 1.0
    if conj:
        y = own ** -b
        X = (1 / (1 - y * y)) ** conj
    W = (1 + own ** -b) if own else 1.0
    for N, k, c in P.fermat_places():
        if k > v:
            W *= (1 + N ** -b) ** c
    W *= P.tail(b)
    T2 = (1 - x) ** -(P.g2() - 1)
    return A * X, T2 * W


def b_cross(P, v, own=None, conj=0):
    def lose(b):
        ax, tw = column_ratio(P, b, v, own, conj)
        return ax <= tw
    return last_true(lose, 1.01, 12.0)


def bracket(P, v, own, conj, lo, hi):
    """(6) on [lo, hi]: A X and T2 W both decrease in b, so A X(b_(i+1))
    > T2 W(b_i) on consecutive grid points leaves no b in [lo, hi] with
    A X <= T2 W. The grid runs past hi. Returns (ok, least ratio)."""
    b = lo + STEP * np.arange(int(np.ceil((hi - lo) / STEP)) + 2)
    assert b[-1] > hi
    ax, tw = column_ratio(P, b, v, own, conj)
    r = ax[1:] / tw[:-1]
    return bool((r > 1).all()), float(r.min())


def b_top(n=2):
    """beta_F2 < b_top: Psi_2 <= zeta_K - 1 <= zeta^n - 1 < 1 past it."""
    return last_true(lambda b: zeta(b) ** n >= 2, 1.0001, 4.0) + STEP


def quad_profile(t2, types):
    two = [(1, 1, 2)] if t2 == "s" else [(2, 1, 1)]
    odd = {}
    for p, t in zip(REL, types):
        odd[p] = {"s": [(p, 2)], "r": [(p, 1)], "i": [(p * p, 1)]}[t]
    return Profile(2, two, odd)


def columns(P):
    """(label, v, own, conj) per column class the profile has."""
    out = [(f"v={v} generic", v, None, 0) for v in (1, 2, 3, 4, 8, 16)]
    for N, k, c in P.fermat_places():
        out.append((f"N={N} own", k, N, c - 1))
    return out


def control():
    print("CONTROL C: the cyclic quartic of conductor 272")
    P = Profile(4, [(4, 1, 1)], {3: [(9, 2)], 5: [(5, 4)], 17: [(17, 1)],
                                 47: [(47, 4)]})
    b2 = 1.54337
    x = 2 ** -b2
    A = 1 + x + x * x
    X = (1 - 47 ** -b2) ** -3
    W = (1 + 5 ** -b2) ** 4 * (1 + 9 ** -b2) ** 2 * (1 + 17 ** -b2)
    ax, tw = column_ratio(P, b2, 1)
    ax *= X
    print(f"  at b = {b2}: A X = {ax:.6f} (by hand {A * X:.6f}), "
          f"T2 W = {tw:.6f} (by hand {W:.6f})")
    bl = P.b_low()
    print(f"  floor b_low = {bl:.4f}")
    check("C the norm-47 column beats the F2 column at its F2 point, and "
          "the floor sits below that point",
          ax < tw and abs(ax - A * X) < 1e-9 and abs(tw / W - 1) < 1e-8
          and bl <= b2)


def control_quadratic():
    print("\nCONTROL C2: the floor under the twelve quadratic fields read")
    pts = {-7: 1.66168, 17: 1.60213, -8: 1.63177, 161: 1.69138,
           24: 1.56837, 156: 1.66461, 204: 1.66621, 264: 1.61826,
           -24: 1.67138, -84: 1.65813, -264: 1.69169, -276: 1.68330}
    ok = True
    for D, b2 in pts.items():
        ty = "".join({1: "s", 0: "r", -1: "i"}[kronecker(D, p)]
                     for p in (2,) + REL)
        bl = quad_profile(ty[0], ty[1:]).b_low()
        ok &= bl is not None and bl <= b2
        print(f"  D={D}: config {ty}, floor {bl:.4f} <= F2 point {b2}")
    check("C2 every floor at or below its field's F2 point", ok)


def main():
    control()
    control_quadratic()
    if not all(CHECKS):
        print("a control failed: no verdict")
        sys.exit(1)
    top = b_top()
    print(f"\nbeta_F2 < {top:.4f} at every quadratic field (zeta^2 < 2)")
    bad, least, worst, nconf = 0, None, {}, 0
    for t2 in ("s", "r"):
        for types in product("sri", repeat=len(REL)):
            nconf += 1
            P = quad_profile(t2, types)
            bl = P.b_low()
            assert bl is not None
            for label, v, own, conj in columns(P):
                ok, r = bracket(P, v, own, conj, bl, top)
                bad += not ok
                key = (t2, label)
                if key not in worst or r < worst[key][0]:
                    worst[key] = (r, "".join(types))
                if least is None or r < least[0]:
                    least = (r, t2, types, label, v, own, conj, bl)
    print(f"R the bracket on [b_low, {top:.4f}] by {STEP}: {bad} failing "
          f"columns of {nconf} configurations")
    print("  least A X(b+) / T2 W(b) per (2's type, class); config = the "
          "types of 3, 5, 17, 257, 65537")
    for key in sorted(worst, key=lambda k: worst[k][0]):
        r, ty = worst[key]
        print(f"  2 {key[0]} {key[1]:13s} {r:.6f}  config {ty}")
    r, t2, types, label, v, own, conj, bl = least
    bc = b_cross(quad_profile(t2, types), v, own, conj)
    print(f"  the least column: 2 {t2}, config {''.join(types)}, {label}: "
          f"ratio {r:.6f}, b_low {bl:.4f}, b_c {bc:.4f}, "
          f"b_low - b_c {bl - bc:+.4f}")
    check("R every column of every configuration beats the F2 column on "
          "[b_low, b_top], bracketed", bad == 0)
    print(f"{sum(CHECKS)} of {len(CHECKS)} checks pass")
    sys.exit(0 if all(CHECKS) else 1)


if __name__ == "__main__":
    main()