# VERIF_FERME3P.py -- verification de la borne en forme close pour trois logarithmes p-adiques
# d'unites principales rationnelles (note NOTE_FERME3P). Python pur, aucun argument, quelques secondes.
# [1] majoration de log F_K ; [2] constantes de la recette et condition lineaire en B ;
# [3] grille : conditions (Z1)-(Z3) et critere verifies directement, toutes les branches majorees ;
# [4] exemple p = 7, alpha = (2, 3, 5) apres changement de base, comparaison avec Yu.
import math, itertools
from fractions import Fraction

M, KAP, ELL = 520, 0.65, 7.2          # recette : K = floor(M L A), chi = KAP sqrt(L), L = ceil(ELL Bc)
B0 = 10.0                             # Bc = max(log B / log p, B0)
LN2 = math.log(2.0)

def say(*a):
    print(*a, flush=True)

def logF_exact(K):
    s = math.fsum((K - k) * math.lgamma(k + 1) for k in range(1, K))
    return -12.0 * s / (K ** 3 - K)

def logF_bound(K):
    return -2.0 * math.log(K) + 11.0 / 3.0 + 4.0 * math.log(K) / K

def tau(L):
    return KAP ** (2 / 3) * M ** (-1 / 3) + (M * L) ** (-1 / 3) + (3.0 / M) ** (1 / 3)

def recipe(a, Bc):
    A = a[0] * a[1] * a[2]
    L = math.ceil(ELL * Bc)
    K = math.floor(M * L * A)
    chi = KAP * math.sqrt(L)
    c1 = (M * KAP) ** (2 / 3) * L
    c2 = (M * L) ** (2 / 3)
    c3 = (3.0 * M * M) ** (1 / 3) * L
    def bx(c):
        return (math.floor(c * a[1] * a[2]), math.floor(c * a[0] * a[2]), math.floor(c * a[0] * a[1]))
    return A, L, K, chi, bx(c1), bx(c2), bx(c3), (c1, c2, c3)

def bl_cor2(p, H1, H2, logB):
    # Bugeaud-Laurent, Corollaire 2 (g = 1, D = 1), log A_i = H_i >= log p
    P = math.log(p)
    assert H1 >= P - 1e-12 and H2 >= P - 1e-12
    lb = logB + math.log(1.0 / H2 + 1.0 / H1)
    mx = max(lb + math.log(P) + 0.4, 10.0 * P, 10.0)
    return 24.0 * p / ((p - 1) * P ** 4) * mx * mx * H1 * H2

def constants():
    L0 = math.ceil(ELL * B0); K0 = math.floor(M * L0)
    t0 = tau(L0)
    th0 = 2.25 * t0 * (1.0 + 1.0 / K0)
    beta0 = 2 * math.log(t0) - 2 * math.log(1 - 1.0 / (M * L0)) + 11.0 / 3.0 + 4 * math.log(K0) / K0
    extra = (max(beta0, 0.0) + 4.5 * math.log(K0 + 1) / K0) / LN2
    slope = ELL * (1 - th0) - 2.0
    c1min = (M * KAP) ** (2 / 3) * L0
    eta = (1 + 1.0 / c1min) ** 2 / (1 - 1.0 / (KAP ** (4 / 3) * M ** (1 / 3) * L0))
    mu = M ** (1 / 3) * KAP ** (-2 / 3)
    C1 = M * (ELL + 1.0 / B0) ** 2 + 0.1
    C2inf = 96.0 * mu * mu * eta * eta
    assert C1 <= C1_NOTE and C2inf + 0.1 <= C2_NOTE
    return dict(L0=L0, K0=K0, th0=th0, beta0=beta0, extra=extra, slope=slope, eta=eta, mu=mu, C1=C1_NOTE, C2inf=C2inf)

C1_NOTE, C2_NOTE = 27711.0, 11107.0    # constantes annoncees dans la note

def C2(p, cst):
    return C2_NOTE * p / (p - 1.0) * (1 + 0.31 / math.log(p)) ** 2

def check_point(p, a, logB, cst):
    # verifie directement les hypotheses du critere et majore chaque branche ; renvoie (v_max, borne annoncee)
    P = math.log(p)
    Bt = logB / P; Bc = max(Bt, B0)
    A, L, K, chi, (R1, S1, T1), (R2, S2, T2), (R3, S3, T3), cs = recipe(a, Bc)
    h = [x * P for x in a]
    V1 = (R1 + 1) * (S1 + 1) * (T1 + 1)
    Mm = max(R1 + S1 + 1, S1 + T1 + 1, R1 + T1 + 1, chi * math.sqrt(V1))
    assert min(R1, S1, T1, R2, S2, T2, R3, S3, T3) >= 1 and K >= 2 and L >= 2
    assert V1 > K * Mm and V1 > L, "Z1"
    assert (R2 + 1) * (S2 + 1) * (T2 + 1) > max(2 * K * L, K * K), "Z2"
    assert (R3 + 1) * (S3 + 1) * (T3 + 1) > 3 * K * K * L, "Z3"
    R, S, T = R1 + R2 + R3 + 1, S1 + S2 + S3 + 1, T1 + T2 + T3 + 1
    N = K * (K + 1) * L // 2
    logb = math.log((R + T - 2) / 2.0) + math.log((S + T - 2) / 2.0) + 2 * logB + logF_bound(K)
    lam = 2.0 * ((L * L) // 4) / L
    E = (2 * K + 1) * (L - 1) / 3.0
    rhs = math.log(N) + 2.0 * (K - 1) / 3.0 * max(logb, 0.0) + lam * ((R - 1) * h[0] + (S - 1) * h[1] + (T - 1) * h[2])
    assert E * P > rhs, "critere"
    assert Mm > max(R1, S1, T1), "M"
    nondeg = K * L - 2
    U1 = (S1 + 1) * (T1 + 1) / (Mm - max(S1, T1))
    U2 = (R1 + 1) * (T1 + 1) / (Mm - max(R1, T1))
    U3 = (R1 + 1) * (S1 + 1) / (Mm - max(R1, S1))
    for j, U in enumerate((U1, U2, U3)):
        assert U <= cst["mu"] * cst["eta"] * a[j] * (1 + 1e-12), "U"
    deg = 0.0
    if U3 >= 1:
        deg = max(deg, bl_cor2(p, U3 * h[0] + U1 * h[2], U3 * h[1] + U2 * h[2], logB))
    if U1 >= 1 and U2 >= 1:
        deg = max(deg, bl_cor2(p, U2 * h[0] + U1 * h[1], U1 * h[2], logB))
    triv = (math.log(2) + max(R1, R2) * h[0] + max(S1, S2) * h[1] + max(T1, T2) * h[2]) / P
    two = max(bl_cor2(p, h[i], h[j], logB) for (i, j) in ((0, 1), (0, 2), (1, 2)))
    one = max((h[j] + math.log(2) + logB) / P for j in range(3))
    v = max(nondeg, deg, triv, two, one) + Bt
    claim = A * Bc * Bc * max(cst["C1"], C2(p, cst) * max(a))
    assert v <= claim, ("borne", p, a, logB, v, claim)
    return v, claim, nondeg, deg

def yu(p, hs, logB):
    n = 3
    if p % 4 == 1:
        d, f, w, c1 = 1, 1, 2, 1473.0
    else:
        d, f, w, c1 = 2, 2, 4, 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(hs[0], f / (10.0 * (n + 4) * d))
    return C * hp * hs[1] * hs[2] * max(logB, math.log(3.0))

EXEMPLES = [   # p, alpha, base de Lambda_0, B = 10 puissance eB, valeurs du tableau : forme close, Yu, rapport
    (5, (2, 3, 7), [(-1, -1, 0), (-1, 0, 1), (-3, 1, 0)], [(10, 1.01e7, 3.08e9, 306), (30, 8.94e7, 9.24e9, 103)]),
    (7, (2, 3, 5), [(-3, 0, 0), (-1, 2, 0), (-1, 1, -1)], [(10, 5.54e6, 4.27e9, 771), (30, 4.99e7, 1.28e10, 257)]),
    (11, (2, 3, 7), [(-2, -1, 0), (-1, 1, -1), (-4, 3, 0)], [(10, 4.35e6, 4.16e9, 956), (30, 3.61e7, 1.25e10, 345)]),
    (13, (2, 3, 7), [(-1, 0, -1), (-4, 1, 0), (-1, 3, -1)], [(10, 3.97e6, 6.37e8, 160), (30, 2.90e7, 1.91e9, 66)]),
    (101, (2, 3, 7), [(-2, -3, 1), (-7, 3, 0), (-5, 1, 4)], [(10, 5.69e6, 2.43e8, 42), (30, 1.25e7, 7.29e8, 58)]),
]

def det3(E):
    return (E[0][0] * (E[1][1] * E[2][2] - E[1][2] * E[2][1]) - E[0][1] * (E[1][0] * E[2][2] - E[1][2] * E[2][0])
            + E[0][2] * (E[1][0] * E[2][1] - E[1][1] * E[2][0]))

def teich_index(p):
    fs = [q for q in range(2, p) if (p - 1) % q == 0 and all(q % r for r in range(2, int(q ** 0.5) + 1))]
    gen = next(x for x in range(2, p) if all(pow(x, (p - 1) // q, p) != 1 for q in fs))
    tab = {}; y = 1
    for e in range(p - 1):
        tab[y] = e; y = y * gen % p
    return lambda r: tab[r % p]

def closed(p, a, logBp, cst):
    Bc = max(logBp / math.log(p), B0)
    return a[0] * a[1] * a[2] * Bc * Bc * max(cst["C1"], C2(p, cst) * max(a))

if __name__ == "__main__":
    say("[1] log F_K <= -2 log K + 11/3 + 4 log K / K")
    worst = max(logF_exact(K) - logF_bound(K) for K in range(2, 3001))
    for K in (37440, 200000):
        worst = max(worst, logF_exact(K) - logF_bound(K))
    say("    max(log F_K - borne) sur 2 <= K <= 3000 et K = 37440, 200000 : %.3e (doit etre < 0)" % worst)
    assert worst < 0

    say("[2] constantes de la recette (M, kappa, ell) = (%d, %.2f, %.1f)" % (M, KAP, ELL))
    cst = constants()
    say("    L0 = %d, K0 = %d, tau0 = %.6f, theta0 = %.6f, beta0 = %.6f" % (cst["L0"], cst["K0"], tau(cst["L0"]), cst["th0"], cst["beta0"]))
    assert cst["L0"] == 72 and cst["K0"] == 37440
    assert tau(72) <= 0.30256 and cst["th0"] <= 0.68078 and cst["beta0"] <= 1.27688 and cst["beta0"] > 0
    pente = ELL * (1 - 0.68078) - 2.0
    rhs10 = 1.0 + cst["extra"]
    say("    avec theta0 <= 0,68078 : ell(1-theta0)-2 >= %.5f ; condition en Bc = 10 : %.4f > %.4f" % (pente, B0 * pente, rhs10))
    assert pente >= 0.2983 and rhs10 <= 2.8440 and 10 * 0.2983 > 2.8440
    assert KAP ** (4 / 3) * M ** (1 / 3) * 72 >= 326 and (M * KAP) ** (2 / 3) <= 48.53
    say("    mu = %.6f, eta = %.6f, 96 mu mu eta eta = %.3f, C1 = %.1f" % (cst["mu"], cst["eta"], cst["C2inf"], cst["C1"]))
    assert cst["mu"] <= 10.7167 and cst["eta"] <= 1.00366 and cst["C2inf"] <= 11106 and M * (ELL + 0.1) ** 2 + 0.1 <= 27711
    for p, val in ((2, 46528), (3, 27390), (5, 19748), (7, 17416), (101, 12776)):
        say("    C2(%d) = %.3f <= %d" % (p, C2(p, cst), val))
        assert C2(p, cst) <= val

    say("[3] grille de controle direct (toutes les permutations des a_j)")
    n = 0; rmax = 0.0
    avals = (1.0, 1.0001, 1.37, 2.0, 3.5, 10.0, 77.0)
    for p in (2, 3, 5, 7, 13, 101, 10007, 1000000007):
        for a in itertools.product(avals, repeat=3):
            for Bt in (0.3, 5.0, 9.99, 10.0, 10.01, 13.7, 31.0, 100.0, 1000.0, 100000.0):
                v, claim, nd, dg = check_point(p, a, Bt * math.log(p), cst)
                n += 1; rmax = max(rmax, v / claim)
    say("    %d points verifies ; max (majoration effective / borne annoncee) = %.4f" % (n, rmax))

    say("[4] exemples du tableau (changement de base, Corollaire 2)")
    for p, alpha, Em, rows in EXEMPLES:
        P = math.log(p); ti = teich_index(p)
        nn = [ti(x) for x in alpha]
        g = p - 1
        index = g // math.gcd(math.gcd(math.gcd(nn[0], nn[1]), nn[2]), g)
        D = det3(Em)
        assert abs(D) == index
        gam = []
        for e in Em:
            assert (nn[0] * e[0] + nn[1] * e[1] + nn[2] * e[2]) % g == 0
            x = Fraction(1)
            for b, k in zip(alpha, e): x *= Fraction(b) ** k
            assert (x.numerator - x.denominator) % p == 0
            gam.append(x)
        hs = [math.log(max(abs(x.numerator), abs(x.denominator))) for x in gam]
        adj = [[0] * 3 for _ in range(3)]
        for r in range(3):
            for c in range(3):
                mi = [[Em[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])
        cs = [sum(abs(Fraction(adj[k][j], D)) for k in range(3)) for j in range(3)]
        a = [max(x, P) / P for x in hs]
        say("  p = %d, alpha = %s : indice %d, gamma = %s, hauteurs %s, c = %s" % (p, alpha, index, [str(x) for x in gam], [round(x, 3) for x in hs], [str(c) for c in cs]))
        for eB, tab_cf, tab_yu, tab_r in rows:
            logBp = eB * math.log(10) + math.log(float(max(cs)))
            bd = closed(p, a, logBp, cst)
            y = yu(p, [math.log(x) for x in alpha], eB * math.log(10))
            say("    B = 1e%d : Bc' = %.2f, forme close %.4e (tableau %.3g) ; Yu %.4e (tableau %.3g) ; rapport %.2f (tableau, partie entiere : %d)" % (eB, max(logBp / P, B0), bd, tab_cf, y, tab_yu, y / bd, tab_r))
            assert bd <= tab_cf <= 1.01 * bd and abs(y / tab_yu - 1) < 0.006 and tab_r <= y / bd < tab_r + 1
    p, al7 = 7, (2, 3, 5)
    a7 = [math.log(8) / math.log(7), math.log(9) / math.log(7), math.log(10) / math.log(7)]
    lo, hi = 30.0, 1.0e6                       # log B ou les deux bornes se croisent (B' = B pour cette base)
    for _ in range(200):
        mid = 0.5 * (lo + hi)
        if closed(7, a7, mid, cst) < yu(7, [math.log(x) for x in al7], mid): lo = mid
        else: hi = mid
    say("  p = 7 : la forme close reste meilleure que Yu tant que log B < %.4g (B < 10 puissance %.4g)" % (lo, lo / math.log(10)))
    assert lo > 1.77e4
    say("FIN : toutes les verifications sont passees.")