cyclotop.py

"""cyclotop.py -- the 2-power tower's own floor: does K = Q(zeta_(2^k))
lose at its own degree D = 2^(k-1)? And the chain over 2 there.

QUESTION. At K = C = Q(zeta_(2^k)), k >= 2, K_ab = K, so L* = L(2^k)
(tower.py's (1)). (A) To what depth J is the place over 2 held by
L(2^k), i.e. the largest j with lambda(P^j) | 2^k, where e = D, outside
the depth step's e > 2^(k-1)? (B) Is the margin's sign at the F2 point
ever negative, at any k?

THE DERIVATION (by hand, before this code).
  (A) The chain. Z_2[zeta] = Z_2[pi], pi = 1 - zeta, Eisenstein of
  degree e = D: (1 - t)^D + 1, lower coefficients c_0 = 2, c_j = (-1)^j
  C(D, j); v(2) = e. U^1/U^j is generated by 1 + pi^i, i < j, so its
  exponent divides 2^k iff j <= a_i for every i < j, a_i = v((1 +
  pi^i)^(2^k) - 1); as a_i > i, J = min over i >= 1 of a_i. Squaring
  1 + y, v(y) = t: v((1 + y)^2 - 1) = v(2y + y^2) = 2t if t < e, t + e
  if t > e. (i) i > e: a_i = i + k e >= (k + 1) e + 1. (ii) i < e not a
  power of 2, R = floor(log2 i): t doubles exactly while below e, first
  exceeds it at 2^(k-1-R) i in (e, 2e), never meeting e, then gains e
  at each of the R + 1 squarings left: a_i = 2^(k-1-R) i + (R + 1) e >
  (R + 2) e. At i = 3 (R = 1, needing 3 < e, k >= 3) that is 1.5 e + 2 e
  = 3.5 e; every other such i has R >= 2, a_i > 4 e. (iii) i = 2^r, 1 <=
  r <= k - 2: zeta^(2^r) = (1 - pi)^(2^r) = 1 + pi^(2^r) + z, each other
  binomial term of valuation (r - v2(j)) e + j > e, so 1 + pi^(2^r) =
  zeta^(2^r) w, v(w - 1) > e, and a_i = v(w - 1) + k e > (k + 1) e. i =
  1: (1 + pi)/(1 - pi) = 1 + 2 pi/(1 - pi), so 1 + pi = zeta w, v(w - 1)
  = e + 1, a_1 = (k + 1) e + 1. i = e: (1 - pi)^e = zeta^e = -1 gives
  pi^e = -2 - sum_(0<j<e) C(e, j) (-pi)^j, the least term at j = e/2 of
  valuation 1.5 e, so 1 + pi^e = -(1 + y), v(y) >= 1.5 e, a_e >= (k +
  1.5) e. So at k >= 3 the least is a_3 = 3.5 e alone ((k + 1) e >= 4 e),
  and at k = 2 (e = 2) a_1 = a_2 = a_3 = 7 = 3.5 e, a_i = i + 4 past.
  So J = 7 * 2^(k-2) at every k >= 2: Q(i)'s 7 and Q(zeta_8)'s 14
  (tower.py's C2) are its k = 2, 3, and J >= 2^k = 2e, as the general
  step (lambda(Q^j) = 2^ceil(log2 j) to j = 2e) already gives.
  (B) The margin. Every odd norm is 1 mod 2^k, so an odd place counts in
  L(2^k) only at norm n_c = 2^k + 1, to depth 1: L* = P_J S_c, P_J = 1 +
  x + ... + x^J. The F2 column holds the places of Fermat norm (n - 1 a
  power of 2) at depth 1 and no other place (there is one place over 2,
  of degree 1): sigma_F2 = S_c W_k, W_k over the norms past n_c, which
  are Fermat primes F_m with 2^m > k at k >= 3 (9 = n_c at k = 3, and 9
  is no norm at k >= 4, 9 != 1 mod 16), D places each, and at k = 2
  also the one place of norm 9 (3 inert in Q(i)). So log(1 + margin) =
  M = log P_J - w_k, S_c cancelling, at the F2 point x, which is below
  1/2 (the point's equation l + log S_3 + w + log(x + rho') = 0, S_3 >=
  1, w, rho' >= 0, as tower.py's (5)). Bound w_k by counting
  every Fermat number past n_c, prime or not: the least, F_n, has 2^n
  >= k + 1, so log2 F_n > k + 1 and D x^(log2 F_n) < 2^(k-1) x^(k+1),
  and F_(n+j) adds D x^(log2 F_(n+j)) < D x^((j+1)(k+1)); so
  w_k < 2^(k-1) x^(k+1)/(1 - x^(k+1)) = x^2 (2x)^(k-1)/(1 - x^(k+1)).
  With P_J >= P_(2^k) (the theorem below does not lean on (A)), l(x) >=
  x + x^2/2 and l_c = log(1/(1 - x^(2^k+1))) <= 1.002 x^9 at k >= 3, x
  <= 1/2:
    M >= x + x^2/2 - (16/15) x^2 - 1.002 x^9 = x (1 - (17/30) x - 1.002
    x^8) > x (1 - 17/60 - 0.004) > 0 on (0, 1/2],
  (2x)^(k-1) <= 1 and 1/(1 - x^(k+1)) <= 16/15 at k >= 3. At k = 2,
  Q(i), the quadratic case already says no loss. So the F2
  place of Q(zeta_(2^k)) tops all but finitely many columns at every k
  >= 2, with no read of the point: M > 0 wherever the point can lie.
  Where it holds and why, against the shape K != C: there d = D m'
  places of each Fermat norm, so w_k ~ m' x^2 (2x)^(k-1), m' unbounded,
  which is what lets it lose; at K = C, m' = 1.

PREDICTIONS (frozen before the engine).
  P1 the a_i engine (Z/2^P[t]/(E), E = (1 - t)^D + 1, P = k + 6) at k =
     2..7 and every i = 1..2e + 2: J = min a_i = 7 * 2^(k-2); a_i = i +
     k e at every i > e; a_i = 2^(k-1-R) i + (R + 1) e at every i < e
     not a power of 2; a_3 = 3.5 e the unique least at k >= 3; a_1 = (k +
     1) e + 1; a_(2^r) > (k + 1) e for 1 <= r <= k - 2; a_e >= (k + 1.5)
     e.
  P2 f_k(x) = log P_(2^k)(x) - w_k^+(x), w_k^+ over every Fermat number
     F_m > 2^k + 1, D places each (and log(1 + x^(log2 9)) at k = 2), is
     positive at every k = 2..64 and every x of a grid of 4001 points
     log-spaced on [1e-8, 1/2], 1/2 included; and the hand bound x (1 -
     (17/30) x - 1.002 x^8) lies at or below f_k and above 0 there at
     every k = 3..64.
  P3 w_k^+(1/2) < 1/4 at every k = 2..64 (largest at k = 2^n - 1, near
     2^(2^n - 2)/(2^(2^n) + 1)), and f_k(1/2) > 0.4.
  KILL (frozen): some k in 2..64 and some grid x with f_k(x) <= 0, read
  on the arm a confound (the largest D against the least Fermat number
  above n_c) inflates: k = 2^n - 1 (3, 7, 15, 31, 63) at x = 1/2, where
  w_k^+ is largest against log P; or the hand bound above f_k at some
  point (the derivation then not carrying the theorem).

CONTROLS (frozen).
  C1 POSITIVE (the kill can fire): the same f with d = 1072 places of
     each Fermat norm in place of D at k = 4 (the shape K != C at D(4),
     where a field loses) is negative at some grid x.
  C2 the a_i engine's J against tower.chain at Q(i) (chain([2, -2], 8):
     1, 2, 4, 4, 4, 4, 4, 8, so J = 7) and Q(zeta_8) (chain([2, -4, 6,
     -4], 15): 8 through 14, 16 at 15, so J = 14).
  C3 the engine on a known K != C chain: E = t^8 - 2 (e = 8) at 2^k = 8
     gives J = 8 (tower.py's C2: 1, 2, 4, 4, 8, 8, 8, 8, 16).
  C4 the engine's precision: at k = 4 every a_i read at P = k + 6 equals
     its read at P = k + 10.

FINDINGS (entered after the run, from its printed output).
  CONTROLS. C1: with 1072 places at k = 4, f is negative at 61 of the
  4001 grid points, from x = 0.3833, so the kill can fire. C2: the
  engine's J is 7 and 14, tower.chain's 7 and 14 (Q(i) 1, 2, 4, 4, 4,
  4, 4, 8; Q(zeta_8) 8 from j = 5 through 14, 16 at 15). C3: 8. C4: the
  18 values at k = 4 agree at P = 10 and 14.
  P1 HOLDS: J = 7, 14, 28, 56, 112, 224 at k = 2..7, the least at i =
  3 alone from k = 3 (at k = 2 at i = 1, 2, 3); every a_i read (i <=
  2e + 2; those past it, i + k e > 3.5 e from k = 3, by (A)(i) alone) on
  its formula; a_1 = a_2 = (k + 1) e + 1 (17, 41, 97, 225, 513 from k = 3);
  a_e = 7, 18, 44, 104, 240, 544 at k = 2..7, (k + 1.5) e exactly.
  P2 HOLDS: f_k > 0 at all 63 x 4001 points, the hand bound in (0, f_k]
  at every point from k = 3; the least f_k/log P_(2^k) 0.639326, at k
  = 63, x = 1/2 (f 0.443147).
  P3 HOLDS: w_k^+(1/2) = 0.227475 at k = 2, 0.244229 at 3, 0.249521 at
  7, 0.249998 at 15, then 1/4 less 2.91e-11 at 31 and 6.78e-21 at 63
  (2^-67, as the count past F_n gives), the largest at k = 63; f_k(1/2)
  0.433924 at k = 2, 0.446963 at 3, 0.443627 at 7, 0.443149 at 15 and
  0.443147 at 31 and 63.
  KILL MISSED: with the derivation's bound, the F2 place of
  Q(zeta_(2^k)) tops all but finitely many columns at every k >= 2, the
  margin positive wherever the point lies in (0, 1/2]; and the place
  over 2 is held by L(2^k) to depth exactly J = 7 * 2^(k-2), a theorem
  by (A), read k = 2..7.

RUN RECORD. `python cyclotop.py` under a 2048 MB cap: 7 of 7 checks,
1.0 s, peak commit 57.6 MB. The engine was appended to the frozen
docstring in one write and run once; P3's reads at x = 1/2 are in
mpmath (at k = 63, 2^64 + 1 is not a float), the grid in floats.
At review: C2 and C3 re-implement tower.chain's ring arithmetic (its
mul and val), so they read the reduction to J, not the arithmetic; the hand
chain at Q(i) (U^1 = mu_4 x U^3) is the independent read. w_plus stops
at F_12, dropping a mass below 2^63 x^8192 < 10^-2400 at x <= 1/2,
under every margin read. Repaired at review, no print moved: three doc
citations in the derivation restated (a script cites no doc); its step
"each later F_m adding at most x^(k+1) times the previous term", false
at k = 2^n - 1, replaced by the bound the count uses; KILL MISSED now
printed only when the hand bound also holds (a scratch copy with hand()
doubled, P2 failing, printed it before the fix and not after).
"""

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")

import mpmath as mp  # noqa: E402

import tower  # noqa: E402

CHECKS = []
KS_CHAIN = range(2, 8)
KS_BOUND = range(2, 65)
MMAX = 13


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)


# ------------------------------------------------------------ (A) chain

def eis(D):
    """Lower coefficients of (1 - t)^D + 1, monic of degree D."""
    return [(2 if j == 0 else (-1) ** j * math.comb(D, j)) for j in range(D)]


class Ring:
    """Z/2^P [t]/(E), E Eisenstein given by its lower coefficients."""

    def __init__(self, E, P):
        self.E, self.e, self.mod = E, len(E), 2 ** P

    def mul(self, a, b):
        e, mod = self.e, self.mod
        r = [0] * (2 * e)
        for i, ai in enumerate(a):
            if ai:
                for j, bj in enumerate(b):
                    if bj:
                        r[i + j] += ai * bj
        for i in range(2 * e - 1, e - 1, -1):
            c = r[i]
            if c:
                r[i] = 0
                for j, Ej in enumerate(self.E):
                    if Ej:
                        r[i - e + j] -= c * Ej
        return [c % mod for c in r[:e]]

    def tpow(self, i):
        out = [1] + [0] * (self.e - 1)
        t = [0, 1] + [0] * (self.e - 2) if self.e > 1 else [-self.E[0]]
        for _ in range(i):
            out = self.mul(out, t)
        return out

    def val(self, a):
        """v(a) in pi, capped at P e (a = 0 mod 2^P reads the cap)."""
        v = math.inf
        for i, c in enumerate(a):
            c %= self.mod
            if c:
                v = min(v, self.e * ((c & -c).bit_length() - 1) + i)
        return min(v, self.e * (self.mod.bit_length() - 1))


def a_list(E, k, P, imax):
    """a_i = v((1 + pi^i)^(2^k) - 1), i = 1..imax."""
    R = Ring(E, P)
    out = []
    tp, t1 = R.tpow(0), R.tpow(1)
    for i in range(1, imax + 1):
        tp = R.mul(tp, t1)
        u = list(tp)
        u[0] = (u[0] + 1) % R.mod
        for _ in range(k):
            u = R.mul(u, u)
        u[0] -= 1
        out.append(R.val(u))
    return out


def j_of_chain(ch, m):
    """The largest j with chain[j - 1] | m."""
    j = 0
    while j < len(ch) and m % ch[j] == 0:
        j += 1
    return j


# ----------------------------------------------------------- (B) margin

def log2_fermat(m):
    if m >= 11:
        return 2.0 ** m
    return 2.0 ** m + math.log1p(2.0 ** -(2 ** m)) / math.log(2)


def w_plus(k, x, places=None):
    """Every Fermat number past 2^k + 1, D (or `places`) places each, and
    at k = 2 the one place of norm 9."""
    D = 2 ** (k - 1) if places is None else places
    lx = math.log(x)
    s = 0.0
    for m in range(MMAX):
        if 2 ** m > k:
            s += math.log1p(math.exp(log2_fermat(m) * lx))
    w = D * s
    if k == 2:
        w += math.log1p(math.exp(math.log2(9) * lx))
    return w


def log_p(n, x):
    """log(1 + x + ... + x^n)."""
    return math.log1p(-math.exp((n + 1) * math.log(x))) - math.log1p(-x)


def f_k(k, x, places=None):
    return log_p(2 ** k, x) - w_plus(k, x, places)


def hand(x):
    return x * (1 - 17 * x / 30 - 1.002 * x ** 8)


# 4001 points log-spaced from 1e-8 to 1/2, 1/2 itself the last
GRID = [math.exp(math.log(1e-8) + (math.log(0.5) - math.log(1e-8)) * i / 4000)
        for i in range(4001)]
GRID[-1] = 0.5


def w_plus_mp(k):
    """w_k^+ and f_k at x = 1/2 in mpmath: x^log2(F) = 1/F there."""
    mp.mp.dps = 60
    D = mp.mpf(2) ** (k - 1)
    s = mp.mpf(0)
    for m in range(MMAX):
        if 2 ** m > k:
            s += mp.log1p(1 / (mp.mpf(2) ** (2 ** m) + 1))
    w = D * s
    if k == 2:
        w += mp.log1p(mp.mpf(1) / 9)
    fval = mp.log(2 * (1 - mp.mpf(2) ** -(2 ** k + 1))) - w
    return w, fval


# ------------------------------------------------------------- the runs

def controls():
    print("=== CONTROLS")
    neg = [x for x in GRID if f_k(4, x, places=1072) <= 0]
    check("C1 POSITIVE: f with d = 1072 places at k = 4 goes negative",
          bool(neg), (f"{len(neg)} of {len(GRID)} grid points, least x "
                      f"{min(neg):.4f}") if neg else "none")
    ch2 = tower.chain([2, -2], 8)
    ch3 = tower.chain([2, -4, 6, -4], 15)
    J2 = min(a_list(eis(2), 2, 8, 6))
    J3 = min(a_list(eis(4), 3, 9, 10))
    print(f"  tower.chain Q(i): {ch2}; Q(zeta_8): {ch3}")
    check("C2 the engine's J against tower.chain at Q(i), Q(zeta_8)",
          J2 == j_of_chain(ch2, 4) == 7 and J3 == j_of_chain(ch3, 8) == 14,
          f"engine {J2}, {J3}; chain {j_of_chain(ch2, 4)}, "
          f"{j_of_chain(ch3, 8)}")
    J8 = min(a_list([-2] + [0] * 7, 3, 9, 18))
    check("C3 t^8 - 2 at 2^k = 8: J = 8", J8 == 8, f"{J8}")
    lo = a_list(eis(8), 4, 4 + 6, 18)
    hi = a_list(eis(8), 4, 4 + 10, 18)
    check("C4 precision: k = 4 at P = k + 6 against k + 10", lo == hi,
          f"{lo}")
    guard()


def part_a():
    print("\n=== (A) the chain at K = Q(zeta_(2^k)), k = 2..7")
    ok = True
    for k in KS_CHAIN:
        e = 2 ** (k - 1)
        a = a_list(eis(e), k, k + 6, 2 * e + 2)
        J = min(a)
        bad = []
        for i, ai in enumerate(a, 1):
            R = int(math.floor(math.log2(i)))
            pow2 = i & (i - 1) == 0
            if i > e and ai != i + k * e:
                bad.append((i, ai, "i > e"))
            if i < e and not pow2 and ai != 2 ** (k - 1 - R) * i + (R + 1) * e:
                bad.append((i, ai, "formula"))
            if i == 1 and ai != (k + 1) * e + 1:
                bad.append((i, ai, "i = 1"))
            if pow2 and 2 <= i <= 2 ** (k - 2) and not ai > (k + 1) * e:
                bad.append((i, ai, "2^r"))
            if i == e and not ai >= (k + 1.5) * e:
                bad.append((i, ai, "i = e"))
        argmin = [i for i, ai in enumerate(a, 1) if ai == J]
        uniq = k == 2 or argmin == [3]
        okk = J == 7 * 2 ** (k - 2) and not bad and uniq
        ok &= okk
        print(f"  k = {k}: e {e}, J = {J} (7 * 2^(k-2) = {7 * 2 ** (k - 2)}),"
              f" argmin {argmin}, a_1..a_8 {a[:8]}, a_e {a[e - 1]}; "
              f"off formula {bad[:3]}")
    check("P1 J = 7 * 2^(k-2) and every a_i as derived, k = 2..7", ok)


def part_b():
    print("\n=== (B) the margin's floor, k = 2..64, 4001 x on [1e-8, 1/2]")
    bad, hbad, worst = [], [], None
    for k in KS_BOUND:
        for x in GRID:
            f = f_k(k, x)
            if f <= 0:
                bad.append((k, x))
            if k >= 3:
                h = hand(x)
                if not (0 < h <= f):
                    hbad.append((k, x))
            r = f / log_p(2 ** k, x)
            if worst is None or r < worst[0]:
                worst = (r, k, x, f)
    print(f"  f_k <= 0 at {bad[:4]}; hand bound off at {hbad[:4]}; least "
          f"f/log P {worst[0]:.6f} at k = {worst[1]}, x = {worst[2]:.4g} "
          f"(f {worst[3]:.6f})")
    check("P2 f_k > 0 at every point; the hand bound in (0, f_k] from k = 3",
          not bad and not hbad)
    rows, okw, okf, top = [], True, True, None
    for k in KS_BOUND:
        w, fv = w_plus_mp(k)
        okw &= w < mp.mpf(1) / 4
        okf &= fv > mp.mpf("0.4")
        if top is None or w > top[0]:
            top = (w, k)
        if k in (2, 3, 7, 15, 31, 63):
            rows.append(f"{k}: w {mp.nstr(w, 12)}, 1/4 - w "
                        f"{mp.nstr(mp.mpf(1) / 4 - w, 4)}, f {mp.nstr(fv, 8)}")
    print("  at x = 1/2: " + "; ".join(rows))
    print(f"  largest w_k^+(1/2) {mp.nstr(top[0], 15)} at k = {top[1]}")
    check("P3 w_k^+(1/2) < 1/4 and f_k(1/2) > 0.4 at every k = 2..64",
          okw and okf)
    if not bad and not hbad:
        print("  KILL MISSED: f_k > 0 at every k = 2..64 on the grid, the arm "
              "k = 2^n - 1 at x = 1/2 included")


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


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