# VERIF_MULT3.py -- ancillary file for the paper "An optimal multiplicity estimate for interpolation determinants in three logarithms"
# David Niedbala Giraudin, ORCID 0009-0009-1526-1178.
# Technical assistance: Claude (Anthropic).
# Plain Python + numpy (no scipy). No argument, no external file. Run in plain Python, not under Sage.
# Size parameter: FULL = False keeps the lightest tests (a few minutes at most on a slow machine).
# Sections:
#   [0] arithmetic self-test
#   [1] Proposition 1: exact order at x = 0 of det(z to the power mu times exp(l x th.z)) computed in F_Q[[x]] (Q a 61-bit prime),
#       generic points and points of the lattice of the kit; it equals Theta_V, beyond Waldschmidt's bound
#   [2] closed forms: Theta_V(N) - sum(k+m) = N(2K+1)(L-1)/6, max D_V = KL-1, Theta_V >= Theta(K-1, .), for 1 <= K, L < 40
#   [3] Theorem 4 (p-adic): exact valuations of integer determinants against the refined bound, p = 2, 3, 5, 7, 11, 13
#   [4] p-adic table: (Z1)-(Z3), refined criterion, cases (2.14) and (2.15), Yu's constant
FULL = False
import math, random, time
import numpy as np
from fractions import Fraction
from math import comb, gcd
def say(*a):
    print(*a, flush=True)

def selftest():
    assert 6 * 7 == 42 and pow(3, 4) == 81
    assert int(np.sum(np.arange(5, dtype=np.int64))) == 10
    assert comb(10, 3) == 120 and gcd(12, 18) == 6
    assert math.isqrt(10 ** 12) == 10 ** 6
    say("[0] arithmetic self-test: OK")
# Exact test of the p-adic analytic bound (no factorial loss):
#   3 logs: v_p(Delta) >= Theta(K-1, N) - sum(k_i + m_i)      [exact relation theta3 = theta1 to the power beta1 times theta2 to the power beta2]
#   2 logs: v_p(Delta) >= N K (L-1) / 2                        [exact relation theta2 = theta1 to the power (b1/b2)]

def theta(K0, I):
    tot = 0; d = 0; left = I
    while left > 0:
        c = min(d, K0) + 1; t = min(c, left); tot += t * d; left -= t; d += 1
    return tot

def vp(n, p):
    if n == 0: return 10 ** 9
    v = 0
    while n % p == 0: n //= p; v += 1
    return v

def det_int(M):
    # Bareiss fraction-free determinant over Z
    M = [row[:] for row in M]; n = len(M); sign = 1; prev = 1
    for k in range(n - 1):
        if M[k][k] == 0:
            sw = next((i for i in range(k + 1, n) if M[i][k] != 0), None)
            if sw is None: return 0
            M[k], M[sw] = M[sw], M[k]; sign = -sign
        for i in range(k + 1, n):
            for j in range(k + 1, n):
                M[i][j] = (M[i][j] * M[k][k] - M[i][k] * M[k][j]) // prev
        prev = M[k][k]
    return sign * M[n - 1][n - 1]

def zp(fr, p, mod):
    # element of Z_(p) as residue mod 'mod'
    return fr.numerator * pow(fr.denominator, -1, mod) % mod

def pw(theta, e, p, M):
    # theta a principal unit modulo p to the power M, e in Z_(p): theta to the power e modulo p to the power M
    mod = p ** M
    return pow(theta, zp(Fraction(e), p, p ** (M - 1)), mod)

def test3(p, K, L, trials, M, rng):
    N = K * (K + 1) * L // 2
    rows = [(k, m, l) for l in range(L) for k in range(K) for m in range(K - k)]
    assert len(rows) == N
    bound = theta(K - 1, N) - sum(k + m for (k, m, l) in rows)
    worst = None
    for _ in range(trials):
        mod = p ** M
        th1 = 1 + p * rng.randrange(1, p ** (M - 1)); th2 = 1 + p * rng.randrange(1, p ** (M - 1))
        b3 = rng.choice([q for q in range(1, 12) if q % p]); b1 = rng.randrange(1, 15); b2 = rng.randrange(1, 15)
        be1, be2 = Fraction(b1, b3), Fraction(b2, b3)
        # N distinct points (r, s, t) in a box
        side = 2
        while side ** 3 < N: side += 1
        pts = rng.sample([(r, s, t) for r in range(side) for s in range(side) for t in range(side)], N)
        Mx = []
        for (k, m, l) in rows:
            row = []
            for (r, s, t) in pts:
                X = r * b3 + t * b1; Y = s * b3 + t * b2           # = b3 z1, b3 z2
                z1 = r + t * be1; z2 = s + t * be2
                val = comb(X, k) * comb(Y, m) * pw(th1, l * z1, p, M) * pw(th2, l * z2, p, M) % mod
                row.append(val)
            Mx.append(row)
        d = det_int(Mx) % mod
        v = vp(d, p) if d else M
        if v < bound:
            return ("FAIL", p, K, L, v, bound)
        worst = v - bound if worst is None else min(worst, v - bound)
    return ("ok", p, K, L, N, bound, worst)

def test2(p, K, L, trials, M, rng):
    N = K * L
    rows = [(k, l) for l in range(L) for k in range(K)]
    bound = N * K * (L - 1) // 2
    worst = None
    for _ in range(trials):
        mod = p ** M
        th1 = 1 + p * rng.randrange(1, p ** (M - 1))
        b2 = rng.choice([q for q in range(1, 12) if q % p]); b1 = rng.randrange(1, 15)
        pts = rng.sample([(r, s) for r in range(N) for s in range(N)], N)
        Mx = []
        for (k, l) in rows:
            row = []
            for (r, s) in pts:
                X = r * b2 + s * b1; z = r + Fraction(s * b1, b2)
                row.append(comb(X, k) * pw(th1, l * z, p, M) % mod)
            Mx.append(row)
        d = det_int(Mx) % mod
        v = vp(d, p) if d else M
        if v < bound:
            return ("FAIL", p, K, L, v, bound)
        worst = v - bound if worst is None else min(worst, v - bound)
    return ("ok", p, K, L, N, bound, worst)


def mixbound3(K, L, e):
    N = K * (K + 1) * L // 2
    return min((N - I) * e + theta(K - 1, I) for I in range(N + 1)) - L * (K - 1) * K * (K + 1) // 3
def mixbound2(K, L, e):
    N = K * L
    return min((N - I) * e + I * (I - 1) // 2 for I in range(N + 1)) - L * K * (K - 1) // 2
def test3e(p, K, L, e, trials, M, rng):
    N = K * (K + 1) * L // 2
    rows = [(k, m, l) for l in range(L) for k in range(K) for m in range(K - k)]
    bound = mixbound3(K, L, e); worst = None
    for _ in range(trials):
        mod = p ** M
        th1 = 1 + p * rng.randrange(1, p ** (M - 1)); th2 = 1 + p * rng.randrange(1, p ** (M - 1))
        eps = p ** e * rng.choice([u for u in range(1, 50) if u % p])     # v(eps) = e exactly
        b3 = rng.choice([q for q in range(1, 12) if q % p]); b1 = rng.randrange(1, 15); b2 = rng.randrange(1, 15)
        be1, be2 = Fraction(b1, b3), Fraction(b2, b3)
        side = 2
        while side ** 3 < N: side += 1
        pts = rng.sample([(r, s, t) for r in range(side) for s in range(side) for t in range(side)], N)
        Mx = []
        for (k, m, l) in rows:
            row = []
            for (r, s, t) in pts:
                X = r * b3 + t * b1; Y = s * b3 + t * b2
                val = comb(X, k) * comb(Y, m) * pw(th1, l * (r + t * be1), p, M) * pw(th2, l * (s + t * be2), p, M)
                val = val * pw(1 + eps, Fraction(l * t, b3), p, M) % mod
                row.append(val)
            Mx.append(row)
        d = det_int(Mx) % mod; v = vp(d, p) if d else M
        if v < bound: return ("FAIL", p, K, L, e, v, bound)
        worst = v - bound if worst is None else min(worst, v - bound)
    return ("ok", p, K, L, "e=%d" % e, bound, worst)
def test2e(p, K, L, e, trials, M, rng):
    N = K * L; rows = [(k, l) for l in range(L) for k in range(K)]
    bound = mixbound2(K, L, e); worst = None
    for _ in range(trials):
        mod = p ** M
        th1 = 1 + p * rng.randrange(1, p ** (M - 1))
        eps = p ** e * rng.choice([u for u in range(1, 50) if u % p])
        b2 = rng.choice([q for q in range(1, 12) if q % p]); b1 = rng.randrange(1, 15)
        pts = rng.sample([(r, s) for r in range(N) for s in range(N)], N)
        Mx = []
        for (k, l) in rows:
            row = []
            for (r, s) in pts:
                X = r * b2 + s * b1
                val = comb(X, k) * pw(th1, l * (r + Fraction(s * b1, b2)), p, M) * pw(1 + eps, Fraction(l * s, b2), p, M) % mod
                row.append(val)
            Mx.append(row)
        d = det_int(Mx) % mod; v = vp(d, p) if d else M
        if v < bound: return ("FAIL", p, K, L, e, v, bound)
        worst = v - bound if worst is None else min(worst, v - bound)
    return ("ok", p, K, L, "e=%d" % e, bound, worst)
# Rank test of the class-restricted interpolation matrix of the p-adic 3-log kit.

Q = (1 << 61) - 1   # prime modulus for rank computation (rank mod Q = N implies rank over Q = N)

def rank_mod(rows, q):
    A = [r[:] for r in rows]; n = len(A); m = len(A[0]); rk = 0; col = 0
    for col in range(m):
        piv = next((i for i in range(rk, n) if A[i][col] % q), None)
        if piv is None: continue
        A[rk], A[piv] = A[piv], A[rk]
        inv = pow(A[rk][col], q - 2, q)
        for i in range(n):
            if i != rk and A[i][col] % q:
                f = A[i][col] * inv % q
                A[i] = [(x - f * y) % q for x, y in zip(A[i], A[rk])]
        rk += 1
        if rk == n: break
    return rk

def dlog(a, p):
    # index of a w.r.t. a primitive root mod p
    for gen in range(2, p):
        if all(pow(gen, (p - 1) // q, p) != 1 for q in range(2, p) if (p - 1) % q == 0 and all(q % d for d in range(2, q))):
            break
    x = 1
    for e in range(p - 1):
        if x == a % p: return e
        x = x * gen % p

def minbox(V):
    a = 0
    while (a + 1) ** 3 <= V: a += 1
    return a

def test(p, K, L, alphas, bs, chi=0.1):
    g = p - 1
    N = K * (K + 1) * L // 2
    # box 1
    a1 = 0
    while not ((a1 + 1) ** 3 > g * K * max(2 * a1 + 1, chi * (a1 + 1) ** 1.5) and (a1 + 1) ** 3 > g * L): a1 += 1
    a2 = 0
    while not ((a2 + 1) ** 3 > 2 * g * K * L and (a2 + 1) ** 3 > g * K * K): a2 += 1
    a3 = 0
    while not ((a3 + 1) ** 3 > 3 * g * K * K * L): a3 += 1
    R = a1 + a2 + a3 + 1
    b1, b2, b3 = bs
    d1, d2 = gcd(b1, b3), gcd(b2, b3)
    m = [dlog(a, p) for a in alphas]
    rows = [(k, mm, l) for l in range(L) for k in range(K) for mm in range(K - k)]
    out = []
    for c in range(g):
        cols = [(r, s, t) for r in range(R) for s in range(R) for t in range(R) if (m[0] * r + m[1] * s + m[2] * t - c) % g == 0]
        M = []
        for (k, mm, l) in rows:
            M.append([comb(r * (b3 // d1) + t * (b1 // d1), k) * comb(s * (b3 // d2) + t * (b2 // d2), mm)
                      * pow(alphas[0], l * r, Q) * pow(alphas[1], l * s, Q) * pow(alphas[2], l * t, Q) % Q for (r, s, t) in cols])
        out.append(rank_mod(M, Q))
    return N, (a1, a2, a3), R, len(cols), out

# p-adic three-log kit (non-degenerate case) versus Yu's explicit bound (as in Alvarado et al. 2019, Thm 3.1).
# Lambda = a1 to the power b1 times a2 to the power b2, minus a3 to the power b3, a_i integers >= 2 prime to p, b_i <= B.  Output: upper bound for v_p(Lambda).

def E(K, L):
    # exact multiplicity (Proposition 1): (2/N)(Theta_V(N) - sum(k+m))
    return (2 * K + 1) * (L - 1) / 3.0

def box(V, la):
    # integer box (R,S,T) with (R+1)(S+1)(T+1) > V and R*la1 ~ S*la2 ~ T*la3 ; returns R,S,T
    G = (la[0] * la[1] * la[2]) ** (1.0 / 3.0)
    s = V ** (1.0 / 3.0)
    while True:
        R = max(0, math.ceil(s * G / la[0]) - 1); S = max(0, math.ceil(s * G / la[1]) - 1); T = max(0, math.ceil(s * G / la[2]) - 1)
        if (R + 1) * (S + 1) * (T + 1) > V:
            return R, S, T
        s *= 1.001

def box1(g, K, L, la):
    s = 2.0
    while True:
        R, S, T = box(s, la)
        V = (R + 1) * (S + 1) * (T + 1)
        M = max(R + S + 1, S + T + 1, R + T + 1)
        if V > g * K * M and V > g * L:
            return R, S, T
        s *= 1.01

def crit(p, g, la, logB, K, L):
    N = K * (K + 1) * L // 2
    R1, S1, T1 = box1(g, K, L, la)
    R2, S2, T2 = box(max(2.0 * g * K * L, g * K * K), la)
    R3, S3, T3 = box(3.0 * g * K * K * L, la)
    R, S, T = R1 + R2 + R3 + 1, S1 + S2 + S3 + 1, T1 + T2 + T3 + 1
    logb = math.log((R + T - 2) / 2.0) + math.log((S + T - 2) / 2.0) + 2 * logB - 2 * math.log(K) + 11.0 / 3.0
    lam = 2.0 * ((L * L) // 4) / L
    lhs = E(K, L) * math.log(p)
    rhs = math.log(N) + 2.0 * (K - 1) / 3.0 * logb + lam * ((R - 1) * la[0] + (S - 1) * la[1] + (T - 1) * la[2])
    return lhs - rhs

def yu(p, alphas, B):
    # Yu (Alvarado et al. Thm 3.1 / Cor 3.3), n = 3, base field Q ; p = 1 mod 4 directly, p = 3 mod 4 via Q(i) (p inert, f = 2)
    n = 3
    if p % 4 == 1:
        d, f, w = 1, 1, 2
        c1 = 1473.0
    else:
        d, f, w = 2, 2, 4
        c1 = 1282.0
    u = 1 if w == 2 else 2
    a1 = 8.0 * (p - 1) / (p - 2)
    k2 = c1 * a1 * n ** n * (n + 1) ** (n + 1) / math.factorial(n)
    k3 = p ** f / 2.0 ** u * (d / (f * math.log(p))) ** (n + 2) * math.log(max(d, math.e))
    k4 = max(math.log(math.e ** 4 * (n + 1) * d), 1.0, f * math.log(p))
    C = (n + 1) * k2 * k3 * k4
    hp = max(math.log(alphas[0]), f / (10.0 * (n + 4) * d))
    Om = hp * math.log(alphas[1]) * math.log(alphas[2])
    return C * Om * math.log(max(B, 3))


def check34(K, L):
    N = K * (K + 1) * L // 2
    best = K * L * N; th = 0; I = 0; d = 0; argI = 0
    while I < N:
        c = min(d, K - 1) + 1
        for _ in range(c):
            if I == N: break
            I += 1; th += d
            v = K * L * (N - I) + th
            if v <= best: best, argI = v, I
        d += 1
    bound = N * N / (2.0 * K) * (1 + 2.0 / L - 6.0 / (K * L) - 1.0 / (3.0 * L * L))
    return best, bound, argI, N

TABLE = [(5, (2, 3, 7), 10, 133, 114445, 8.0, [('C', 27, 560264), ('C', 27, 107905)]), (5, (2, 3, 7), 30, 521, 259645, 12.0, [('C', 82, 1788651), ('C', 82, 327479)]), (7, (2, 3, 5), 10, 106, 91045, 8.0, [('C', 23, 424378), ('C', 23, 97837)]), (7, (2, 3, 5), 30, 415, 196589, 12.0, [('C', 68, 1301202), ('C', 68, 287300)]), (11, (2, 3, 7), 10, 133, 102041, 8.0, [('C', 19, 745874), ('C', 19, 143039)]), (11, (2, 3, 7), 30, 210, 553429, 20.0, [('C', 54, 2198788), ('C', 61, 355363)]), (13, (2, 3, 7), 10, 133, 106032, 8.0, [('C', 17, 896020), ('C', 19, 154270)]), (13, (2, 3, 7), 30, 235, 522896, 20.0, [('C', 54, 2188619), ('C', 54, 399914)]), (101, (2, 3, 7), 10, 149, 3384594, 30.0, [('B', 0, 0), ('B', 0, 0)]), (101, (2, 3, 7), 30, 295, 4144076, 30.0, [('B', 0, 0), ('B', 0, 0)])]


def box1chi(g, K, L, la, chi):
    s = 2.0
    while True:
        R, S, T = box(s, la)
        V = (R + 1) * (S + 1) * (T + 1)
        M = max(R + S + 1, S + T + 1, R + T + 1, chi * math.sqrt(V))
        if V > g * K * M and V > g * L:
            return R, S, T, M
        s *= 1.01

def parts(p, g, la, logB, K, L, chi):
    N = K * (K + 1) * L // 2
    R1, S1, T1, M = box1chi(g, K, L, la, chi)
    V1 = (R1 + 1) * (S1 + 1) * (T1 + 1)
    assert V1 > g * K * M and V1 > g * L
    R2, S2, T2 = box(max(2.0 * g * K * L, g * K * K), la)
    assert (R2 + 1) * (S2 + 1) * (T2 + 1) > g * max(2 * K * L, K * K)
    R3, S3, T3 = box(3.0 * g * K * K * L, la)
    assert (R3 + 1) * (S3 + 1) * (T3 + 1) > 3 * g * K * K * L
    R, S, T = R1 + R2 + R3 + 1, S1 + S2 + S3 + 1, T1 + T2 + T3 + 1
    logb = math.log((R + T - 2) / 2.0) + math.log((S + T - 2) / 2.0) + 2 * logB - 2 * math.log(K) + 11.0 / 3.0
    lam = 2.0 * ((L * L) // 4) / L
    margin = E(K, L) * math.log(p) - (math.log(N) + 2.0 * (K - 1) / 3.0 * logb + lam * ((R - 1) * la[0] + (S - 1) * la[1] + (T - 1) * la[2]))
    U1 = (S1 + 1) * (T1 + 1) / (M - max(S1, T1))
    U2 = (R1 + 1) * (T1 + 1) / (M - max(R1, T1))
    U3 = (R1 + 1) * (S1 + 1) / (M - max(R1, S1))
    triv = (max(R1, R2) * la[0] + max(S1, S2) * la[1] + max(T1, T2) * la[2]) / math.log(p)
    return margin, (U1, U2, U3), triv

def bl3(p, h1, h2, B):
    lp = math.log(p)
    A1, A2 = max(h1, lp), max(h2, lp)
    bprime = B / A2 + B / A1
    m = max(math.log(bprime) + math.log(lp) + 0.4, 10 * lp, 10.0)
    return 24.0 * p / lp ** 4 * m * m * A1 * A2

def cor3_at(p, h1, h2, B, L, K):
    # two-logarithm criterion at parameters (L, K): returns (margin of the criterion, bound KL-2, bound of the non-injective branch)
    g = p - 1 if p > 2 else 1
    lp = math.log(p)
    def bx(V):
        s = math.sqrt(V * h2 / h1); R = max(1, int(s)); S = int(V / R) + 1
        while R * S <= V: S += 1
        return R, S
    R1, S1 = bx(g * L); R2, S2 = bx(g * (K - 1) * L)
    assert R1 * S1 > g * L and R2 * S2 > g * (K - 1) * L
    R, S = R1 + R2 - 1, S1 + S2 - 1; N = K * L
    logb = math.log(max(R + S - 2, 1) * B / 2.0) - math.log(K - 1) + 1.5 + math.log(K - 1) / K + 1e-6
    lam = 2.0 * ((L * L) // 4) / L
    m = K * (L - 1) * lp - (math.log(N) + (K - 1) * logb + lam * ((R - 1) * h1 + (S - 1) * h2))
    branch = ((p - 1) * ((R2 - 1) * h1 + (S2 - 1) * h2) + math.log(2) + math.log(B)) / lp
    return m, K * L - 2, branch

QM = (1 << 61) - 1


def DV(K, L):
    return sorted(a + c for a in range(K) for c in range((K - a) * L))

def thetaW(K0, I):
    tot = 0; d = 0; left = I
    while left > 0:
        c = min(d, K0) + 1; t = min(c, left); tot += t * d; left -= t; d += 1
    return tot

def inv(a): return pow(a, QM - 2, QM)

def ord_det(rows, D):
    # rows: list of power series (list of length D over F_QM) matrix entries M[i][j]; compute ord of det via valuation pivoting
    M = [[e[:] for e in r] for r in rows]
    n = len(M); total = 0
    for k in range(n):
        # find pivot with minimal order in submatrix
        best = None
        for i in range(k, n):
            for j in range(k, n):
                e = M[i][j]
                o = next((t for t in range(D) if e[t]), D)
                if best is None or o < best[0]:
                    best = (o, i, j)
        o, pi, pj = best
        if o >= D: return None
        M[k], M[pi] = M[pi], M[k]
        for r in M: r[k], r[pj] = r[pj], r[k]
        total += o
        piv = M[k][k]
        # piv = x to the power o times a unit; inverse of the unit
        unit = piv[o:] + [0] * o
        uinv = [0] * D; uinv[0] = inv(unit[0])
        for t in range(1, D):
            s = 0
            for q in range(1, t + 1): s += unit[q] * uinv[t - q]
            uinv[t] = (-s * uinv[0]) % QM
        for i in range(k + 1, n):
            e = M[i][k]
            # factor = e / piv; e has order at least o
            eo = e[o:] + [0] * o
            f = [0] * D
            for t in range(D):
                if eo[t]:
                    for q in range(D - t): f[t + q] = (f[t + q] + eo[t] * uinv[q]) % QM
            for j in range(k, n):
                a = M[k][j]; b = M[i][j]
                for t in range(D):
                    if f[t]:
                        for q in range(D - t): b[t + q] = (b[t + q] - f[t] * a[q]) % QM
    return total

def series_entry(mu1, mu2, l, z1, z2, th, D):
    # (x z1) to the power mu1 times (x z2) to the power mu2 times exp(l x (th1 z1 + th2 z2)), modulo x to the power D
    w = (th[0] * z1 + th[1] * z2) % QM
    e = [0] * D; fac = 1; pw = 1
    for n in range(D):
        if n > 0: fac = fac * n % QM; pw = pw * (l * w) % QM
        e[n] = pw * inv(fac) % QM
    c = pow(z1, mu1, QM) * pow(z2, mu2, QM) % QM
    s = mu1 + mu2
    out = [0] * D
    for n in range(D - s): out[n + s] = e[n] * c % QM
    return out


def thetaV(K, L, I):
    return sum(DV(K, L)[:I])

def mixV(K, L, e):
    N = K * (K + 1) * L // 2
    return min((N - I) * e + thetaV(K, L, I) for I in range(N + 1)) - L * (K - 1) * K * (K + 1) // 3

if __name__ == "__main__":
    t0 = time.time()
    selftest()
    # [1]
    rng = random.Random(3)
    for (K, L) in ((2, 2), (2, 3), (3, 2), (2, 4), (3, 3), (4, 2)):
        N = K * (K + 1) * L // 2
        rows_idx = [(k, m, l) for l in range(L) for k in range(K) for m in range(K - k)]
        tv = thetaV(K, L, N); tw = thetaW(K - 1, N)
        D = tv + 4
        th = (rng.randrange(1, QM), rng.randrange(1, QM))
        pts = [(rng.randrange(1, QM), rng.randrange(1, QM)) for _ in range(N)]
        og = ord_det([[series_entry(k, m, l, z1, z2, th, D) for (z1, z2) in pts] for (k, m, l) in rows_idx], D)
        i3 = inv(3); side = 2
        while side ** 3 < N: side += 1
        lat = rng.sample([(r, s_, t) for r in range(side) for s_ in range(side) for t in range(side)], N)
        pl = [((r + t * 5 * i3) % QM, (s_ + t * 7 * i3) % QM) for (r, s_, t) in lat]
        ol = ord_det([[series_entry(k, m, l, z1, z2, th, D) for (z1, z2) in pl] for (k, m, l) in rows_idx], D)
        assert og == tv and ol == tv, (K, L, og, ol, tv)
        say("[1] K=%d L=%d N=%d : Waldschmidt %d, Theta_V %d, exact order (generic points) %d, (lattice points) %d  [%.0f s]" % (K, L, N, tw, tv, og, ol, time.time() - t0))
    # [2]
    for K in range(1, 40):
        for L in range(1, 40):
            N = K * (K + 1) * L // 2
            Dv = sorted(DV(K, L)); skm = L * (K - 1) * K * (K + 1) // 3
            assert len(Dv) == N and 6 * (sum(Dv) - skm) == N * (2 * K + 1) * (L - 1) and max(Dv) == K * L - 1
            pref = 0
            for I in range(1, N + 1):
                pref += Dv[I - 1]
                if I % max(1, N // 50) == 0 or I == N:
                    assert pref >= thetaW(K - 1, I)
    say("[2] closed forms and comparison with Waldschmidt checked for 1 <= K, L < 40  [%.0f s]" % (time.time() - t0))
    # [3]
    rng = random.Random(5); n = 0; eq = 0
    th_old, mb_old = theta, mixbound3
    for p in (2, 3, 5, 7, 11, 13):
        for (K, L) in ((2, 2), (2, 3), (3, 2), (2, 4), (3, 3)):
            N = K * (K + 1) * L // 2
            theta = (lambda K_, L_: (lambda K0, I: thetaV(K_, L_, I)))(K, L)
            b = thetaV(K, L, N) - L * (K - 1) * K * (K + 1) // 3
            r = test3(p, K, L, 3, b + 6, rng); n += 1; eq += (r[-1] == 0)
            assert r[0] == 'ok', r
            mixbound3 = lambda K_, L_, e_: mixV(K_, L_, e_)
            for e in (1, 2, 4, 7):
                r = test3e(p, K, L, e, 2, max(mixV(K, L, e), 2) + 6, rng); n += 1; eq += (r[-1] == 0)
                assert r[0] == 'ok', r
            theta, mixbound3 = th_old, mb_old
    say("[3] Theorem 4: %d series of exact determinants, no valuation below the refined bound; equality in %d series  [%.0f s]" % (n, eq, time.time() - t0))
    # [4]
    for (p, al, eB, L, K, chi, meth) in TABLE:
        g = p - 1 if p > 2 else 1
        la = [math.log(a) for a in al]; B = 10.0 ** eB; lp = math.log(p)
        m, U, triv = parts(p, g, la, eB * math.log(10.0), K, L, chi)
        assert m > 0, (p, eB, m)
        U1, U2, U3 = U
        hh = ((U3 * la[0] + U1 * la[2], U3 * la[1] + U2 * la[2]), (U2 * la[0] + U1 * la[1], U1 * la[2]))
        degs = []
        for (h1, h2), (kind, L2, K2) in zip(hh, meth):
            vals = [bl3(p, h1, h2, B)]
            if kind == 'C':
                m2, main2, br2 = cor3_at(p, h1, h2, B, L2, K2)
                assert m2 > 0, (p, eB, L2, K2, m2)
                vals.append(max(main2, br2))
            degs.append(min(vals) + math.log(B) / lp)
        deg = max(degs); nd = K * L - 2
        tot = max(nd, deg, triv); y = yu(p, al, B)
        say("[4] p=%3d B=1e%d : (L,K,chi)=(%d,%d,%g) margin %.3g | non-deg. %.3g | (2.15) %.3g | (2.14) %.3g | total %.3g | Yu %.3g | ratio %.2f"
            % (p, eB, L, K, chi, m, nd, deg, triv, tot, y, y / tot))
    say("ALL CHECKS PASSED.  [%.0f s]" % (time.time() - t0))