# VERIF_KIT3P.py (version 3) -- verificateur de la note "Un kit pour les formes lineaires en trois logarithmes p-adiques"
# Python pur + numpy (aucun scipy). Aucun argument, aucun fichier externe. A lancer en Python pur, pas dans Sage.
# Parametre de taille : FULL = False garde les tests les plus legers (quelques minutes au plus sur une machine lente).
# Sections :
#   [0] auto-test arithmetique
#   [1] Proposition (multiplicite) : ordre exact en x = 0 calcule dans F_Q[[x]] (Q premier de 61 bits),
#       points generiques et points du reseau du kit ; il vaut Theta_V
#   [2] formes closes : Theta_V(N) - somme(k+m) = N(2K+1)(L-1)/6, max D_V = KL-1, Theta_V >= Theta(K-1, .), 1 <= K, L < 40
#   [3] Theoreme 1 : valuations exactes de determinants entiers, trois logarithmes (p = 2..13) et deux logarithmes (p = 2..11)
#   [5] Proposition (lemme de zeros dans une classe) : rang N des matrices restreintes a chaque classe
#   [6] Lemme 10 de Bugeaud-Laurent sur les translates de tous les sous-reseaux d'indice g (choix adverse des points)
#   [7] Mode reseau : base exacte du reseau des exposants de partie de Teichmuller triviale (gamma_j = 1 mod p),
#       indice, bornes B', et la table complete pour les trois orientations du terme isole
#   [4] Table : (Z1)-(Z3), critere, cas (2.14) ; cas (2.15) par le meilleur de BL Thm 3 et du Corollaire 3
#       (avec le Lemme 10 de BL sur un translate de sous-reseau d'indice g) ; constante de Yu
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] auto-test arithmetique : OK")
# Tests exacts de la minoration analytique p-adique sans perte (fonctions de base ; en section [3], la
# borne a trois logarithmes est appelee avec Theta_V a la place de Theta(K-1, .)) :
#   3 logarithmes : v_p(Delta) >= Theta_V(N) - somme(k_i + m_i)   [relation exacte theta3 = theta1 puissance beta1 fois theta2 puissance beta2]
#   2 logarithmes : v_p(Delta) >= N K (L-1) / 2                    [relation exacte theta2 = theta1 puissance (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 unite principale modulo p puissance M, e dans Z_(p) : theta puissance e modulo p puissance 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 puissance b1 fois a2 puissance b2, moins a3 puissance b3, a_i integers >= 2 prime to p, b_i <= B.  Output: upper bound for v_p(Lambda).

def E(K, L):
    # multiplicite exacte (Proposition 1) : (2/N)(Theta_V(N) - somme(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, 119, 82607, 5.0, [('C', 27, 385765), ('C', 28, 72021)]), (5, (2, 3, 7), 30, 263, 352066, 12.0, [('C', 82, 985383), ('C', 82, 181204)]), (7, (2, 3, 5), 10, 106, 60054, 5.0, [('C', 23, 276599), ('C', 23, 64272)]), (7, (2, 3, 5), 30, 370, 142897, 8.0, [('C', 68, 814225), ('C', 68, 180514)]), (11, (2, 3, 7), 10, 133, 65669, 5.0, [('C', 19, 478879), ('C', 19, 92659)]), (11, (2, 3, 7), 30, 235, 281543, 12.0, [('C', 54, 1245331), ('C', 61, 201993)]), (13, (2, 3, 7), 10, 60, 160400, 8.0, [('C', 18, 513419), ('C', 18, 99294)]), (13, (2, 3, 7), 30, 235, 286069, 12.0, [('C', 54, 1300444), ('C', 54, 238768)]), (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, swap=False):
    # Corollaire 3 aux parametres (L, K) : renvoie (marge du critere, borne KL-2, borne de la branche non injective)
    # swap=True : roles de gamma_1 et gamma_2 echanges (cas ou p divise le second exposant), boites echangees
    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
    if swap:
        R1, S1, R2, S2, h1, h2 = S1, R1, S2, R2, h2, h1
    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
    c1 = min(lam * (R - 1), L * (R + g - 1) / 2.0 - g * N * L / (6.0 * (S + g - 1)))
    c2 = min(lam * (S - 1), L * (S + g - 1) / 2.0 - g * N * L / (6.0 * (R + g - 1)))
    m = K * (L - 1) * lp - (math.log(N) + (K - 1) * logb + c1 * h1 + c2 * 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 puissance o fois une unite ; inverse de l'unite
        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]
            # facteur = e / piv ; e est d'ordre au moins 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) puissance mu1 fois (x z2) puissance mu2 fois exp(l x (th1 z1 + th2 z2)), modulo x puissance 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

def G_BL(N, L, R, S, g):
    return N * L * (R + g - 1) / 8.0 - g * N * N * L / (24.0 * (S + g - 1))
def maxdev(rs, K, L):
    # max over subsets of size N=KL and permutations of |sum l_i r_sigma(i) - (L-1)/2 sum r| ; exact for chosen multiset by sorting;
    N = K * L
    rs = sorted(rs)
    best = 0.0
    # adversarial subsets: take a smallest, N-a largest, a = 0..N
    for a in range(0, N + 1):
        sub = sorted(rs[:a] + rs[len(rs) - (N - a):]) if N - a > 0 else sorted(rs[:a])
        if len(sub) != N: continue
        l = [i // K for i in range(N)]
        M = (L - 1) / 2.0 * sum(sub)
        hi = sum(li * ri for li, ri in zip(l, sub)); lo = sum(li * ri for li, ri in zip(l, sub[::-1]))
        best = max(best, hi - M, M - lo)
    return best


# ---------------------------------------------------------------- mode reseau (changement de base, g = 1)
LATTABLE = [(5, (2, 3, 7), 10, [(-1, -1, 0), (-1, 0, 1), (-3, 1, 0)], [((1, 2, 0), 85, 53195, 2.0), ((0, 2, 1), 85, 53137, 2.0), ((0, 1, 2), 95, 48605, 2.0)]), (5, (2, 3, 7), 30, [(-1, -1, 0), (-1, 0, 1), (-3, 1, 0)], [((1, 2, 0), 331, 111216, 3.0), ((0, 2, 1), 235, 169302, 5.0), ((0, 1, 2), 235, 169257, 5.0)]), (7, (2, 3, 5), 10, [(-3, 0, 0), (-1, 2, 0), (-1, 1, -1)], [((1, 2, 0), 85, 30295, 1.5), ((0, 2, 1), 67, 40054, 2.0), ((0, 1, 2), 67, 40034, 2.0)]), (7, (2, 3, 5), 30, [(-3, 0, 0), (-1, 2, 0), (-1, 1, -1)], [((1, 2, 0), 235, 84680, 3.0), ((0, 2, 1), 235, 84659, 3.0), ((0, 1, 2), 263, 80150, 3.0)]), (11, (2, 3, 7), 10, [(-2, -1, 0), (-1, 1, -1), (-4, 3, 0)], [((1, 2, 0), 60, 31789, 1.5), ((0, 2, 1), 67, 29356, 1.5), ((0, 1, 2), 54, 38868, 2.0)]), (11, (2, 3, 7), 30, [(-2, -1, 0), (-1, 1, -1), (-4, 3, 0)], [((1, 2, 0), 167, 88458, 3.0), ((0, 2, 1), 167, 88434, 3.0), ((0, 1, 2), 210, 76116, 3.0)]), (13, (2, 3, 7), 10, [(-1, 0, -1), (-4, 1, 0), (-1, 3, -1)], [((1, 2, 0), 54, 29331, 1.5), ((0, 2, 1), 60, 26833, 1.5), ((0, 1, 2), 54, 32465, 2.0)]), (13, (2, 3, 7), 30, [(-1, 0, -1), (-4, 1, 0), (-1, 3, -1)], [((1, 2, 0), 149, 81221, 3.0), ((0, 2, 1), 149, 81195, 3.0), ((0, 1, 2), 167, 73192, 3.0)]), (101, (2, 3, 7), 10, [(-2, -3, 1), (-7, 3, 0), (-5, 1, 4)], [((1, 2, 0), 43, 40705, 5.0), ((0, 2, 1), 43, 40613, 5.0), ((0, 1, 2), 38, 61185, 8.0)]), (101, (2, 3, 7), 30, [(-2, -3, 1), (-7, 3, 0), (-5, 1, 4)], [((1, 2, 0), 85, 97521, 5.0), ((0, 2, 1), 85, 97497, 5.0), ((0, 1, 2), 133, 74136, 5.0)])]

def lat_is_prime(n):
    return n > 1 and all(n % q for q in range(2, int(n ** 0.5) + 1))

def lat_prim_root(p):
    fs = [q for q in range(2, p) if (p - 1) % q == 0 and lat_is_prime(q)]
    for a in range(2, p):
        if all(pow(a, (p - 1) // q, p) != 1 for q in fs):
            return a

def lat_dlog(a, p, gen):
    x = 1
    for e in range(p - 1):
        if x == a % p:
            return e
        x = x * gen % p

def lat_det3(M):
    return (M[0][0] * (M[1][1] * M[2][2] - M[1][2] * M[2][1]) - M[0][1] * (M[1][0] * M[2][2] - M[1][2] * M[2][0])
            + M[0][2] * (M[1][0] * M[2][1] - M[1][1] * M[2][0]))

def lat_numden(al, e):
    num, den = 1, 1
    for a, k in zip(al, e):
        if k > 0: num *= a ** k
        elif k < 0: den *= a ** (-k)
    return num, den

def lat_inv_l1(E):
    D = lat_det3(E)
    adj = [[0] * 3 for _ in range(3)]
    for r in range(3):
        for c in range(3):
            mi = [[E[i][j] for j in range(3) if j != c] for i in range(3) if i != r]
            adj[c][r] = (-1) ** (r + c) * (mi[0][0] * mi[1][1] - mi[0][1] * mi[1][0])
    inv = [[Fraction(adj[i][j], D) for j in range(3)] for i in range(3)]
    return [sum(abs(inv[k][j]) for k in range(3)) for j in range(3)]

def lat_two_log(p, g, h1, h2, B):
    lp = math.log(p); lB = math.log(B); best = None
    for L in list(range(2, 40)) + list(range(40, 400, 7)):
        def ok(K):
            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)
            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
            lamL = 2.0 * ((L * L) // 4) / L
            c1 = min(lamL * (R - 1), L * (R + g - 1) / 2.0 - g * N * L / (6.0 * (S + g - 1)))
            c2 = min(lamL * (S - 1), L * (S + g - 1) / 2.0 - g * N * L / (6.0 * (R + g - 1)))
            return K * (L - 1) * lp - (math.log(N) + (K - 1) * logb + c1 * h1 + c2 * h2), R2, S2
        lo, hi = 3, 3
        while ok(hi)[0] <= 0:
            hi *= 2
            if hi > 10 ** 12: break
        if hi > 10 ** 12: continue
        while hi - lo > 1:
            mid = (lo + hi) // 2
            if ok(mid)[0] > 0: hi = mid
            else: lo = mid
        m, R2, S2 = ok(hi)
        if m <= 0: continue
        tot = max(hi * L - 2, ((p - 1) * ((R2 - 1) * h1 + (S2 - 1) * h2) + math.log(2) + lB) / lp)
        if best is None or tot < best: best = tot
    return best

def lat_bl3(p, g, h1, h2, B):
    lp = math.log(p); A1, A2 = max(h1, lp), max(h2, lp)
    m = max(math.log(B / A2 + B / A1) + math.log(lp) + 0.4, 10 * lp, 10.0)
    return 24.0 * p * g / (p - 1) / lp ** 4 * m * m * A1 * A2

def lat_check():
    out = []
    for (p, al, eB, Est, per) in LATTABLE:
        g = p - 1; gen = lat_prim_root(p)
        n = [lat_dlog(a, p, gen) for a in al]
        for a, k in zip(al, n):
            assert pow(gen, k, p) == a % p
        d = math.gcd(math.gcd(math.gcd(n[0], n[1]), n[2]), g); index = g // d
        E = [list(e) for e in Est]
        for e in E:
            assert (n[0] * e[0] + n[1] * e[1] + n[2] * e[2]) % g == 0
            num, den = lat_numden(al, e)
            assert (num - den) % p == 0                       # gamma_j = num/den = 1 mod p : unite principale
        assert abs(lat_det3(E)) == index                      # base du reseau Lambda_c (indice exact)
        hs = []
        for e in E:
            num, den = lat_numden(al, e); hs.append(math.log(max(num, den)))
        cb = lat_inv_l1(E); B = 10.0 ** eB; Bp = float(max(cb)) * B
        lp = math.log(p); worst = 0.0
        for (order, L, K, chi) in per:
            la = [hs[j] for j in order]
            m, U, triv = parts(p, 1, la, math.log(Bp), K, L, chi)
            assert m > 0, (p, eB, order, m)
            U1, U2, U3 = U; dg = 0.0
            for (h1, h2, extra) in ((U3 * la[0] + U1 * la[2], U3 * la[1] + U2 * la[2], math.log(max(U3, 1)) / lp),
                                    (U2 * la[0] + U1 * la[1], U1 * la[2], math.log(max(U1, 1)) / lp)):
                t = lat_two_log(p, 1, h1, h2, Bp)
                v = lat_bl3(p, 1, h1, h2, Bp) if t is None else min(lat_bl3(p, 1, h1, h2, Bp), t)
                dg = max(dg, v + extra + math.log(Bp) / lp)
            tot = max(K * L - 2, dg, triv + math.log(2) / lp) + math.log(Bp) / lp
            worst = max(worst, tot)
        # un b'_j nul : forme a deux logarithmes (hauteurs h_i, h_j) ; deux nuls : un seul logarithme
        for (i, j) in ((0, 1), (0, 2), (1, 2)):
            t = lat_two_log(p, 1, hs[i], hs[j], Bp)
            v2 = lat_bl3(p, 1, hs[i], hs[j], Bp) if t is None else min(lat_bl3(p, 1, hs[i], hs[j], Bp), t)
            assert v2 + 2 * math.log(Bp) / lp <= worst, (p, eB, i, j, v2, worst)
        for h in hs:
            assert (h + math.log(2) + math.log(Bp)) / lp <= worst
        y = yu(p, al, B)
        out.append((p, eB, [round(h, 3) for h in hs], Bp, worst, y))
        say("[7] mode reseau p=%3d B=1e%d : base %s (indice %d, gamma_j = 1 mod p), hauteurs %s, B' = %.3g ; total %.4g ; Yu %.4g ; rapport %.1f"
            % (p, eB, E, index, [round(h, 3) for h in hs], Bp, worst, y, y / worst))
    return out


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, ordre exact (points generiques) %d, (points du reseau) %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] formes closes et comparaison avec Waldschmidt verifiees pour 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] Theoreme 1, trois logarithmes : %d series de determinants exacts, aucune valuation sous la borne ; egalite dans %d series  [%.0f s]" % (n, eq, time.time() - t0))
    # [3 bis] deux logarithmes
    rng2 = random.Random(7); n2 = 0
    for p in (2, 3, 5, 7, 11):
        for (K, L) in ((2, 2), (2, 3), (3, 2), (3, 3), (2, 5)):
            r = test2(p, K, L, 3, K * L * K * (L - 1) // 2 + 8, rng2)
            assert r[0] == "ok", r
            n2 += 1
    rng2 = random.Random(11)
    for p in (2, 3, 5, 7):
        for (K, L) in ((2, 3), (3, 3)):
            for e in (1, 2, 3, 5):
                r = test2e(p, K, L, e, 3, max(mixbound2(K, L, e), 2) + 6, rng2)
                assert r[0] == "ok", r
                n2 += 1
    say("[3] deux logarithmes : %d series, aucune valuation sous la borne  [%.0f s]" % (n2, time.time() - t0))
    # [5]
    for (p, K, L, al, bs) in ((5, 3, 5, (2, 3, 7), (97, 131, 211)), (3, 4, 5, (2, 5, 7), (97, 131, 211))):
        N, boxes, R, ncols, ranks = test(p, K, L, al, bs)
        assert all(rk == N for rk in ranks), (p, K, L, ranks)
        say("[5] p=%d K=%d L=%d : N=%d, boites %s, grande boite %d, %d colonnes par classe ; rang %s  [%.0f s]" % (p, K, L, N, boxes, R, ncols, ranks, 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)
                m2s, main2s, br2s = cor3_at(p, h1, h2, B, L2, K2, swap=True)
                assert m2 > 0 and m2s > 0 and abs(m2 - m2s) < 1e-6 * max(1.0, abs(m2)) + 1e-3, (p, eB, L2, K2, m2, m2s)
                vals.append(max(main2, br2, main2s, br2s))
            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) marge %.3g | non deg. %.3g | (2.15) %.3g | (2.14) %.3g | total %.3g | Yu %.3g | rapport %.2f"
            % (p, eB, L, K, chi, m, nd, deg, triv, tot, y, y / tot))
    # [6] Lemme 10 de BL sur les translates de tous les sous-reseaux d'indice g (forme normale de Hermite)
    rng6 = random.Random(1); n6 = 0; w6 = 0.0
    for g in (2, 3, 4, 6, 8, 12):
        for gp in [d for d in range(1, g + 1) if g % d == 0]:
            gpp = g // gp
            for x in range(gpp):
                for (K, L) in ((2, 2), (3, 2), (2, 3), (3, 3)):
                    N = K * L
                    for (R, S) in ((5, 7), (8, 4), (10, 10), (6, 15)):
                        cr, cs = rng6.randrange(gp), rng6.randrange(gpp)
                        pts = [(r, s_) for r in range(R) for s_ in range(S)
                               if (r - cr) % gp == 0 and (s_ - cs - x * ((r - cr) // gp)) % gpp == 0]
                        if len(pts) < N: continue
                        d = maxdev([r for r, s_ in pts], K, L); Gb = G_BL(N, L, R, S, g)
                        assert d <= Gb + 1e-9, (g, gp, x, K, L, R, S, d, Gb)
                        n6 += 1
                        if Gb > 0: w6 = max(w6, d / Gb)
    say("[6] Lemme 10 de BL, translates de sous-reseaux d'indice g : %d cas, rapport ecart/borne maximal %.4f  [%.0f s]" % (n6, w6, time.time() - t0))
    lat_check()
    say("TOUT EST VERIFIE.  [%.0f s]" % (time.time() - t0))