#!/usr/bin/env python3
"""
PROVENANCE: NUMERICAL

Computation 153 -- The L_comm rate is exactly Theta(1/sqrt D), and why it looked better
=======================================================================================
Closes a question comp 150 deliberately left open.

THE QUESTION.  Along the matched-scaling sequence, comp 150's seventeen points looked
like they were trending down: 3.67, 2.26, 2.46, 2.86, 0.38, 1.46, ... 0.59.  If
G = |gap| * N really tends to zero, then gamma_D = o(1/sqrt D), better than the rate
the book claims.  Comp 150 refused to claim it -- seventeen noisy points where m(k)
moves in integer steps are not a proof.  The refusal was right, and this computation
says why.

FINDING 1 -- THE RATE IS EXACTLY Theta(1/sqrt D).  At fixed k, G converges to a
NONZERO constant as N grows, even though m(k)/N -> 0 over the same range:

    k = 1 :  G = 1.0161, 1.0128, 1.0106, 1.0091, 1.0079, 1.0070, 1.0063, 1.0057
             over N = 16..44, with m/N falling 0.062 -> 0.023

Monotone, and plainly converging to about 1.00 rather than to 0.  Same at k = 2 and
k = 4, around 0.87 and 0.95, with no downward drift.  So gap ~ C/N with C = O(1), and

    gamma_D  ~  C / (2 sqrt D)  =  Theta(1 / sqrt D),

not o(1/sqrt D).  The book's rate is correct AND tight; the apparent improvement was
not there.  Note this does not conflict with comp 150: the Englis-Zhu relative deficit
IS O(m/N) and does vanish at rate D^(-1/4).  A vanishing relative deficit on the norms
is consistent with an O(1/N) gap in the ratio, and that is what happens.

FINDING 2 -- THE TREND WAS SAWTOOTH PHASE, NOT DECAY.  G is smooth in k within a fixed
m(k) stratum and jumps at stratum boundaries, so a window that straddles several strata
samples different phases and can look like a trend in either direction.  Comp 150's
seventeen points span m = 4..8.

FINDING 3 -- THE CONTROLLING VARIABLE IS THE TRUNCATION HEADROOM, and G IS UNBOUNDED.
Comp 150 identified m(k)/N.  The sharper variable is

    s  =  N - 2(m(k) + 1),

the room the top bridge monomial w^(m+1) has inside the degree-N truncation.  m/N is a
proxy that works because m/N -> 1/2 is the same event as s -> 0.  As s closes and goes
negative, G does not approach a maximum -- it grows without bound:

    N = 40, m = 20, s = -2 :  G = 3.9, 7.8, 11.6, 15.5, 19.3, 23.1, 26.9, 30.7, 34.4, 38.1
                              for k = 106..115

COROLLARY, and it corrects how comp 150 phrased things.  Comp 150 reported "G = 22.0"
as the worst case outside the operative regime.  That number is a property of how far
that sweep happened to sample, not of the object: this computation reaches 38.1 simply
by looking further.  There is no maximum to quote.  The honest statement is that G is
bounded on the operative regime and unbounded off it, and the boundary is s, not a
number.

WHAT THIS DOES NOT DO.  It does not supply the analytical constant, which stays open.
It removes a false hope and names the right variable.
"""
import math
import sys

import numpy as np
import numpy.linalg as la

sys.path.insert(0, "computations-se")
from computation_59_se import (  # noqa: E402
    bergman_basis, J_matrices, D_alpha, T_bridge, lcomm_ratio, m_of_k,
)

I2 = np.eye(2, dtype=complex)


def G_of(k, N, cache):
    if N not in cache:
        b = bergman_basis(N)
        Jp, Jm, Jz = J_matrices(b, 0.0)
        cache[N] = (b, D_alpha(Jp, Jm, Jz))
    b, D = cache[N]
    T = T_bridge(b, 0.0, k)
    if la.norm(T) < 1e-12:
        return None
    return abs(lcomm_ratio(T, D, len(b), I2) - 2.0 * math.sqrt(k)) * N


def main():
    cache = {}
    print("=" * 94)
    print("  Computation 153 -- the L_comm rate is exactly Theta(1/sqrt D)")
    print("=" * 94)

    print()
    print("  (1) At FIXED k, does G tend to zero or to a nonzero constant as N grows?")
    print("      %4s %6s %8s %10s" % ("k", "m(k)", "m/N", "G"))
    for k in (1, 2, 4):
        vals = []
        for N in (16, 20, 24, 28, 32, 36, 40, 44):
            G = G_of(k, N, cache)
            if G is None:
                continue
            vals.append(G)
            print("      %4d %6d %8.3f %10.5f" % (k, m_of_k(k), m_of_k(k) / N, G))
        print("      -> k=%d settles near %.3f, NOT 0\n" % (k, sum(vals[-3:]) / 3))
    print("      gap ~ C/N with C = O(1)  =>  gamma_D = Theta(1/sqrt D), not o(1/sqrt D).")
    print("      The book's rate is correct and tight; the apparent improvement was phase.")

    print()
    print("  (2) The controlling variable is the truncation headroom s = N - 2(m+1),")
    print("      and G is UNBOUNDED as it closes -- there is no maximum to quote.")
    for N in (30, 40):
        print("      N = %d:" % N)
        print("        %5s %4s %5s %10s" % ("k", "m", "s", "G"))
        shown = 0
        for k in range(1, 400):
            m = m_of_k(k)
            s = N - 2 * (m + 1)
            if s != -2:
                continue
            G = G_of(k, N, cache)
            if G is None:
                continue
            print("        %5d %4d %5d %10.4f" % (k, m, s, G))
            shown += 1
            if shown >= 10:
                break
    print()
    print("      comp 150 reported 22.0 as the worst case off the operative regime.")
    print("      That was a property of its sampling window: looking further reaches 38.")
    print("      The boundary is s, not a number.")

    print()
    print("  CONCLUSION")
    print("    gamma_D = Theta(1/sqrt D) exactly.  o(1/sqrt D) is refuted, not merely")
    print("    unproven.  The analytical constant remains open.")


if __name__ == "__main__":
    main()