towerpast.py

"""towerpast.py -- the 2-power tower past its reads: a lemma that every d
= 0 mod 4 from some d0 <= 1000 loses at k = 3, and the turn of the bound
at k = 8, where D(k) as the first degree past d* is untested.

QUESTION. (A) With Q(zeta_8) in K, K_ab = Q(zeta_8), type (1, 0), does
the law fixed below Y = 10^5 (tower.py's arm Z) lose at every d = 0 mod
4 from a d0 <= 1000, by a lemma as sillclass's at k = 2? (B) At k = 8
(Q(zeta_256), d = 0 mod 128), does the bound at count 0 turn at the
first such d past d* = l(x0)/u_8(x0), and does a certificate there read
the loss?

THE DERIVATION (by hand, before this code).
  (A) At k = 3 the places of norm 17 number d, weight 1 (at k = 2 those
  of norm 9 numbered d/2): x1 solving d 17^-b = lam = log(1/x) would
  give c_z >= 1/(1 - y) > 1, the slope 1 - c_z - c_s negative, so the
  lemma would say nothing (a TRANSPLANT from k = 2 that fails here). So
  x1 solves d x^a17 = h lam, h in (0, 1) a free parameter, a17 = log2
  17, y = 17^-b = x^a17. With d = h lam 17^b at x1:
    lz <= h lam [1/(1 - y) + sum_{n = 257, 65537} (17/n)^b/(1 - n^-b)
          + (17/Y)^b (1 + Y/(b - 1))/(1 - Y^-b)] + E,
  the last the tail (d places per prime past Y, each of norm >= p), E
  over the places of norm (p^f)^m' below Y off W_3: at most 4 per p,
  their norms distinct prime powers >= 9^m' (3 has order 2 mod 8, its
  two primes of Z[zeta_8] of norm 9, n_c's), so E <= 4 9^-s (1 + 9/(s -
  1))/(1 - 9^-s), s = m' b, m' = d/4 (integral test).
    s_/x <= h lam [x^(a17 - 1)/(1 - y^2) + sum_n (34/n^2)^b/(1 - n^-2b)
          + (34/Y)^b (1 + Y/(b - 1))/(1 - Y^-b)] + E/x,
  s_ = lz - sw + lt (a W place adds -log(1 - y^2) <= y^2/(1 - y^2), the
  others whole; 1/x = 2^b folded in). g_up = l + lz + lt + log(1 - (1 -
  x)e^-s_) <= l + lz + lt + log x + s_/x. Every coefficient rises with x
  (17/n < 1, 34/n^2 < 1, the tail's factors), so on x1 <= X they are
  bounded at X, l(x1) <= l(X), lam >= lam0 = log(1/X): g_up(x1) <= l(X)
  - (1 - c_z - c_s) lam0 + E' (c_z, c_s at X, h included), needing 1 -
  c_z - c_s > 0. w >= d log(1 + y) >= h lam (1 - y/2): M(x1) <= l(X) - h
  lam0 (1 - y0/2), y0 = X^a17; M'(x1) <= 1/(1 - X) - a17 h lam0/(X(1 +
  y0)), the l_c and other W terms of M' negative. x1 <= X iff d >= d0 =
  h lam0 X^-a17. Where g_up(x1) < 0 the point is above x1, at most x_c;
  M rises then falls once on [0, x_c], so M'(x1) < 0 puts x1 past the
  peak and the margin is at most M(x1) < 0. Hand sizes at X = 1/4, h =
  1/3: d0 = 134, c_z ~ 0.336, c_s ~ 0.009, g_up <= -0.66, M <= -0.17,
  M' <= -6.
  (B) At k = 8, n_c = 257 cancels, W_8 = {65537} and tower.wk's composite
  F_5..F_9; x0^16 ~ 2.0e-7, d* ~ 2.35e6. At count 0, M(x_c) = h_c(x_c)
  turns where x_c crosses h_c's root, x0 shifted by l_c = -log(1 -
  x^257) (below 1e-100): the turn is d* up to that shift, so the kill
  can fire only through an engine fault or a rounding at the step. The
  certificate at m' ~ 18,000: every place off W has norm at least
  9^m', so the law is the bound's but for the tail past Y; at Y = 10^30
  that tail is about d Y^(1 - b)/(b - 1), b = log2(1/x0) ~ 1.389, about
  2e-5, against a bound at the first d past d* of order (d - d*)
  times M's slope in d, at most 128 u_8 ~ 3e-5: the certificate may
  fail at D(8) for the tail alone, read at Y = 10^60 beside it.

PREDICTIONS (frozen before the engine).
  P1 (A) at X = 1/4, h = 1/3: d0 = ceil(h lam0 X^-a17) in [120, 150],
     slope 1 - c_z - c_s >= 0.6, g_up <= -0.5, M <= -0.1, M' <= -4.
  P2 (A) at every d = 0 mod 4 from d0 to 1000, at the arm-Z law: x1 <=
     X, x1 <= x_low, C <= M(x1) <= the M bound, g_up(x1) and M'(x1)
     below theirs.
  P3 (B) bound_at(8, d, 0) over d = 0 mod 128, 40 steps either side of
     d*, is positive below the first d past d* and at most 0 from it.
  P4 (B) the ARM Y certificate (tower.BigLaw at Y = 10^30) at k = 8 is
     negative with M' < 0 and x_c < 1/2 at the first d past d* + 128 and
     at every printed d past it (a guess, 50%, at the first d itself).
  KILL of (A): g_up >= 0, M >= 0, M' >= 0 or slope <= 0 printed at X, or
  d0 > 1000, or P2 failing. KILL of the first-degree pattern at k = 8: a
  turn of P3's count-0 bound more than one step from the first d past
  d*.

CONTROLS (frozen).
  C1 POSITIVE: the same lemma code at k = 2's constants (a9, h = 1/2,
     X = 1/7, the Fermat norms 17, 257, 65537, base 5, factor 2)
     reproduces sillclass.lemma()'s c_z, c_s, g_up, M and M' to 1e-12.
  C2 PLANT: at h = 1 (the transplant) the slope at X = 1/4 is <= 0, so
     the lemma can say no.
  C3 the k = 3 coefficients against a brute sum over the arm-Z law's
     places at d = 1000, x = x1: lz + lt <= c_z lam + E and (lz - sw +
     lt)/x <= c_s lam + E/x, coefficients at X and at x1.
  C4 at k = 8, bound_at(8, d, 0) equals TowerLaw(8, law of exactly
     those counts, no tail).m(x_c) at d = d*'s step, to 1e-12.

FINDINGS (entered after the run, from its printed output).
  CONTROLS. C1: 0.0, d0 930. C2: slope -0.036139 (c_z 1.010738), so the
  lemma can say no. C3: at d = 1000, lz + lt 0.605076 against c_z lam
  0.610882 (at X) and 0.605259 (at x1); (lz - sw + lt)/x 2.2708e-03
  against 1.5352e-02 and 2.2708e-03 (the norm-17 term all but exact at
  x1). C4: 4.4e-16 at d = 2344064. C1 passes sillclass's own E (base 5,
  factor 4: the frozen "factor 2" misnamed it), so e3 is read by no
  control; e3(d0, X) bounds E(d, x1)(1 + 1/x1) since (2/9^m')^b falls in
  b and in m'; it prints 2.644e-63.
  (A) P1 HOLDS: d0 = 134, c_z 0.336913, c_s 0.008467, slope 0.654620,
  E' 2.644e-63; g_up <= -0.619814, M <= -0.173617, M' <= -6.1958. P2
  HOLDS at the 217 d = 0 mod 4 from 136 to 1000, the largest M(x1)
  -0.178151 at d = 136 (x1 0.249055, x_low 0.287237, C -0.495215).
  KILL (A) MISSED: the arm-Z law loses at every d = 0 mod 4 >= 136, so
  with tower.py's certificates from 28, a field containing Q(zeta_8),
  K_ab = Q(zeta_8), type (1, 0), loses at every d = 0 mod 4 >= 28.
  (B) d* = 2343958.59; the first d = 0 mod 128 past it is 2344064 (0.8235
  of a step past d*). P3 HOLDS: the count-0 bound is positive at every
  read d from 2338944 to 2343936 (+1.113e-06 there) and at most 0 from
  2344064 (-5.193e-06) through the read's last d, 40 steps past. KILL
  (B) MISSED. P4 FAILS as frozen:
  at Y = 10^30 the certificate is positive at each printed d from
  2344064 (+2.161e-05) to 2344320 (+8.986e-06), the tail past Y as the
  derivation priced it, and negative at 2345344 (-4.150e-05), the d
  between unread (the P4 print's bracket names the first d, not where P4
  fails, which is 2344192); at Y = 10^60 it is negative with M' <
  0 and x_c < 1/2 at every printed d from 2344064 (-4.973e-06) and
  positive at 2343936 (+1.333e-06). So at k = 8 the loss is read at the
  printed degrees 2344064, +128, +256 and +1280 at Y = 10^60, and the
  no-loss below the turn at count 0 alone (the full scan over the count
  of 257 unread).

RUN RECORD. `python towerpast.py` under a 2048 MB cap: 4 of 4 checks,
1.1 s, peak commit 58.9 MB. The engine was appended to the frozen
docstring in one write and run once. Moved in from the lab unchanged;
every line printed as before.
"""

import math
import os
import sys
import time

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

from scipy.optimize import brentq  # noqa: E402

import sillclass as sc  # noqa: E402
import splitsharp as ss  # noqa: E402
import tower as tw  # noqa: E402

CHECKS = []
A17 = math.log2(17)
X3, H3 = 0.25, 1 / 3


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


def guard():
    if not all(CHECKS):
        print("a control failed: no verdict")
        sys.exit(1)


# ------------------------------------------------------------ the lemma

def coeffs(x, c, wB, a, B, norms, Y):
    """(c_z, c_s) at x: d = c lam B^b, wB d places of norm B, d of each
    norm in norms, d per prime past Y."""
    b = -math.log2(x)
    y = B ** -b
    fac = (1 + Y / (b - 1)) / (1 - Y ** -b)
    cz = c * (wB / (1 - y) + sum((B / n) ** b / (1 - n ** -b) for n in norms)
              + (B / Y) ** b * fac)
    cs = c * (wB * x ** (a - 1) / (1 - y * y)
              + sum((2 * B / n ** 2) ** b / (1 - n ** (-2 * b))
                    for n in norms)
              + (2 * B / Y) ** b * fac)
    return cz, cs


def lemma(X, c, wB, a, B, norms, Y, E):
    lam0 = -math.log(X)
    y0 = X ** a
    cz, cs = coeffs(X, c, wB, a, B, norms, Y)
    lx = -math.log1p(-X)
    return dict(cz=cz, cs=cs, slope=1 - cz - cs, E=E,
                d0=math.ceil(c * lam0 * X ** -a),
                G=lx - (1 - cz - cs) * lam0 + E,
                M=lx - wB * c * lam0 * (1 - y0 / 2),
                Mp=1 / (1 - X) - a * wB * c * lam0 / (X * (1 + y0)))


def e3(d0, X):
    """E + E/x at k = 3: at most 4 places per p off W_3, their norms
    distinct prime powers >= 9^m', m' >= d0/4."""
    s = (d0 / 4) * -math.log2(X)
    return 4 * 9.0 ** -s * (1 + 9 / (s - 1)) / (1 - 9.0 ** -s) * (1 + 1 / X)


def lemma3(h=H3, X=X3):
    L = lemma(X, h, 1.0, A17, 17, (257, 65537), ss.Y, 0.0)
    return lemma(X, h, 1.0, A17, 17, (257, 65537), ss.Y, e3(L["d0"], X))


def x1(d, c, a):
    return brentq(lambda x: d * x ** a + c * math.log(x), 1e-12, 0.5,
                  xtol=1e-18, rtol=1e-15, maxiter=500)


# ----------------------------------------------------------------- main

def controls():
    print("=== CONTROLS")
    want = sc.lemma()
    E = sc.lemma_e(sc.lemma_d0() // 2, sc.LX0)
    got = lemma(sc.LX0, 1.0, 0.5, sc.A9, 9, sc.FERMAT, sc.LY, E)
    err = max(abs(got[k] - want[k]) for k in ("cz", "cs", "G", "M", "Mp"))
    check("C1 POSITIVE: k = 2's constants reproduce sillclass.lemma()",
          err < 1e-12 and got["d0"] == sc.lemma_d0(),
          f"{err:.1e}, d0 {got['d0']}")
    P = lemma3(h=1.0)
    check("C2 PLANT: h = 1 at X = 1/4 gives slope <= 0", P["slope"] <= 0,
          f"slope {P['slope']:+.6f}, c_z {P['cz']:.6f}")
    d = 1000
    q = tw.TowerLaw(3, tw.law(3, d))
    x = x1(d, H3, A17)
    lam = -math.log(x)
    lz, sw, w, lt = q.base.parts(x)
    E = e3(d, x)
    ok = True
    for tag, xx in (("X", X3), ("x1", x)):
        cz, cs = coeffs(xx, H3, 1.0, A17, 17, (257, 65537), ss.Y)
        good = (lz + lt <= cz * lam + E
                and (lz - sw + lt) / x <= cs * lam + E)
        ok = ok and good
        print(f"  d = 1000, coefficients at {tag}: lz + lt {lz + lt:.6f} "
              f"<= {cz * lam:.6f}; (lz - sw + lt)/x "
              f"{(lz - sw + lt) / x:.4e} <= {cs * lam:.4e}")
    check("C3 the k = 3 coefficients bound the arm-Z law's sums at d = "
          "1000, x1", ok)
    ds = tw.dstar(8)
    d = 128 * math.ceil(ds / 128)
    v, _, xc = tw.bound_at(8, d, 0)
    w, _ = tw.caps(8, d)
    W = [(math.log2(m), c, 1) for m, c in w] + [(math.log2(257), 0, 1)]
    e4 = abs(v - tw.TowerLaw(8, ss.Law(d, W, tail=False)).m(xc))
    check("C4 k = 8: bound_at(8, d, 0) = M(x_c) on its law", e4 < 1e-12,
          f"{e4:.1e} at d = {d}")
    guard()


def part_a():
    print("\n=== (A) k = 3: the lemma at X = 1/4, h = 1/3")
    L = lemma3()
    print(f"  d0 = ceil(h lam0 X^-a17) = {L['d0']}; c_z {L['cz']:.6f}, "
          f"c_s {L['cs']:.6f}, slope {L['slope']:.6f}, E' {L['E']:.3e}")
    print(f"  g_up <= {L['G']:+.6f}; M <= {L['M']:+.6f}; M' <= "
          f"{L['Mp']:+.4f}")
    p1 = (120 <= L["d0"] <= 150 and L["slope"] >= 0.6 and L["G"] <= -0.5
          and L["M"] <= -0.1 and L["Mp"] <= -4)
    d0 = L["d0"] + (-L["d0"]) % 4
    bad, worst, n = [], None, 0
    for d in range(d0, 1001, 4):
        q = tw.TowerLaw(3, tw.law(3, d))
        C, dm, xl, xc = q.cert()
        x = x1(d, H3, A17)
        ok = (x <= X3 and x <= xl and C <= q.m(x) <= L["M"]
              and q.g_up(x) <= L["G"] and q.dm(x) <= L["Mp"])
        n += 1
        if not ok:
            bad.append(d)
        row = (q.m(x), q.g_up(x), q.dm(x), d, x, xl, C)
        worst = row if worst is None or row[0] > worst[0] else worst
    m, g, dmx, d, x, xl, C = worst
    print(f"  at the {n} d = 0 mod 4 from {d0} to 1000: failing at "
          f"{bad[:6]}; largest M(x1) {m:+.6f} at d = {d} (x1 {x:.6f}, "
          f"x_low {xl:.6f}, C {C:+.6f}, g_up(x1) {g:+.6f}, M'(x1) "
          f"{dmx:+.4f})")
    p2 = not bad
    kill = (L["G"] >= 0 or L["M"] >= 0 or L["Mp"] >= 0 or L["slope"] <= 0
            or L["d0"] > 1000 or not p2)
    print(f"  P1 {'holds' if p1 else 'FAILS'}; P2 "
          f"{'holds' if p2 else 'FAILS'}; KILL (A) "
          f"{'fires' if kill else 'missed'}")


def part_b():
    print("\n=== (B) k = 8: the turn at count 0 and the certificate")
    ds = tw.dstar(8)
    first = 128 * math.ceil(ds / 128)
    print(f"  d* = {ds:.2f}; the first d = 0 mod 128 past it {first} "
          f"({(first - ds) / 128:.4f} steps past d*)")
    vals = {d: tw.bound_at(8, d, 0)[0] for d in
            range(first - 40 * 128, first + 40 * 128 + 1, 128)}
    pos = [d for d in vals if vals[d] > 0]
    neg = [d for d in vals if vals[d] <= 0]
    turn = min(neg) if neg else None
    for d in (first - 256, first - 128, first, first + 128):
        print(f"  d = {d}: bound at count 0 {vals[d]:+.3e}")
    p3 = bool(pos) and max(pos) < first and turn == first
    print(f"  positive at {min(pos) if pos else None}.."
          f"{max(pos) if pos else None}, at most 0 from {turn}; P3 "
          f"{'holds' if p3 else 'FAILS'}")
    kill = turn is None or abs(turn - first) > 128
    print(f"  KILL (B) {'fires' if kill else 'missed'}")
    certs = {}
    for Yc, tag in ((10 ** 30, "10^30"), (10 ** 60, "10^60")):
        for d in (first - 128, first, first + 128, first + 256,
                  first + 1280):
            C, dm, xl, xc = tw.TowerLaw(8, tw.BigLaw(tw.law(8, d), Yc,
                                                     128)).cert()
            certs[(tag, d)] = (C, dm, xl, xc)
            print(f"  Y = {tag}, d = {d}: certificate {C:+.3e} (M' "
                  f"{dm:+.3f}, x_low {xl:.7f}, x_c {xc:.7f})")
    p4 = all(certs[("10^30", d)][0] < 0 and certs[("10^30", d)][1] < 0
             and certs[("10^30", d)][3] < 0.5
             for d in (first + 128, first + 256, first + 1280))
    c0 = certs[("10^30", first)][0]
    print(f"  P4 {'holds' if p4 else 'FAILS'} (at the first d itself: "
          f"{c0:+.3e})")


def main():
    t0 = time.time()
    controls()
    part_a()
    part_b()
    print(f"\n{sum(CHECKS)} of {len(CHECKS)} checks, "
          f"{time.time() - t0:.1f} s")
    return 0 if all(CHECKS) else 1


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