"""
Cota INCONDICIONAL del margen de Fejer  R_N(t) - C_N(t)  (Parte 20, seccion 8, Proposicion 'fejer').

Identidad (Proposicion 'fejer', ventana de Fejer, L = log N):
    R_N(t) - C_N(t) = sum_gamma G^(gamma),   G^(z) = (F(z - t) + F(z + t))/2,   F(z) = 4 sin^2(L z/2)/z^2.
Ceros en la recta (gamma real): G^(gamma) >= 0.
Ceros fuera de la recta: por Platt-Trudgian (2021) todos los ceros con 0 < Im rho <= H = 3*10^12 estan en la recta, asi
que un cero fuera tiene |gamma_0| > H.  Para |t| <= H/2:
    |F(z)| <= (2 + N^{|d|} + N^{-|d|})/|z|^2 <= (3 + N^{|d|})/|z|^2,   |z| >= |gamma_0| - |t| >= |gamma_0|/2,
    |G^(gamma)| <= 4 (3 + N^{|d|})/gamma_0^2,   d = Re rho - 1/2   (G^ real en los cuartetos +-gamma_0 +- i d).
(zeta no tiene ceros reales en 0 < sigma < 1, asi que todo cero con |Im rho| <= H esta en la recta.)
Sumando sobre los posibles ceros fuera (gamma_0 > H y gamma_0 < -H):
    R_N(t) - C_N(t) >= -E_N,   E_N = 8 sum_{Im rho > H} (3 + N^{d(rho)})/(Im rho)^2.
Cota simple:  d < 1/2,  E_N <= 8 (3 + sqrt N) S_H,  S_H = sum_{Im rho > H} (Im rho)^{-2}.
Cota con la region sin ceros de Mossinghoff-Trudgian-Yang:  d < 1/2 - 1/(5.558691 log Im rho).
Sumas sobre ceros:  para f >= 0 decreciente,  sum_{Im rho > H} f(Im rho) <= int_H^inf (-f'(u)) N^+(u) du,
con N^+(u) = (u/2pi) log(u/(2 pi e)) + 7/8 + 0.112 log u + 0.278 log log u + 2.510 + 1/(5u)   (Trudgian 2014).
Este script calcula S_H y E_N (simple y con la region sin ceros), la longitud N* hasta la que E_N es menor que el margen
medido de zeta, y dibuja E_N frente al margen.
Salidas en resultados_incondicional/: cota_incondicional.log y cota_incondicional.png
"""
import os
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

BASE = os.path.dirname(os.path.abspath(__file__))
CARPETA = os.path.join(BASE, "resultados_incondicional")
H = 3.0e12
C_MTY = 5.558691
MARGEN_MEDIDO = 0.0584      # minimo del margen de Fejer de zeta medido (0 <= t < 120, N <= e^12)
MARGEN_MEDIO = 0.0924       # 4 sum_{gamma > 0} gamma^{-2}: valor medio en t = 0 bajo RH


def N_mas(u):
    """cota superior explicita de N(u) (Trudgian 2014, u >= e)."""
    return u / (2 * np.pi) * np.log(u / (2 * np.pi * np.e)) + 7 / 8 + 0.112 * np.log(u) + 0.278 * np.log(np.log(u)) + 2.510 + 1 / (5 * u)


def suma_sobre_ceros(f, df):
    """cota de sum_{Im rho > H} f(Im rho) para f >= 0 decreciente: int_H^inf (-f'(u)) N^+(u) du (cambio u = H e^s)."""
    def integrando(s):
        u = H * np.exp(s)
        return -df(u) * N_mas(u) * u
    valor, error = quad(integrando, 0.0, 200.0, limit=400, epsabs=0.0, epsrel=1e-12)
    return valor + abs(error)


def S_H():
    return suma_sobre_ceros(lambda u: u ** -2.0, lambda u: -2.0 * u ** -3.0)


def E_simple(N):
    return 8.0 * (3.0 + np.sqrt(N)) * S_H()


def E_region(N):
    """con d < 1/2 - 1/(C log u):  f(u) = (1 + N^{1/2 - 1/(C log u)})/u^2."""
    L = np.log(N)

    def f(u):
        return (3.0 + np.exp(L * (0.5 - 1.0 / (C_MTY * np.log(u))))) / u ** 2

    def df(u):
        e = np.exp(L * (0.5 - 1.0 / (C_MTY * np.log(u))))
        return -2.0 * (3.0 + e) / u ** 3 + e * L / (C_MTY * u * np.log(u) ** 2) / u ** 2
    # f es decreciente para u > H en todo el rango de N usado (se comprueba en una malla)
    malla = H * np.exp(np.linspace(0.0, 50.0, 2000))
    assert np.all(df(malla) < 0), "f no es decreciente, la cota de Stieltjes no se aplica"
    return 8.0 * suma_sobre_ceros(f, df)


def escribir(texto, archivo):
    print(texto, flush=True)
    archivo.write(texto + "\n")
    archivo.flush()


def main():
    os.makedirs(CARPETA, exist_ok=True)
    log = open(os.path.join(CARPETA, "cota_incondicional.log"), "w")
    s = S_H()
    escribir(f"H = {H:.0e} (Platt-Trudgian), S_H = sum_(Im rho > H) (Im rho)^-2 <= {s:.6e}", log)
    escribir(f"   comparacion: log(H/2pi)/(pi H) = {np.log(H / (2 * np.pi)) / (np.pi * H):.6e}", log)
    escribir(f"Teorema (incondicional): para todo N > 1 y |t| <= {H / 2:.1e}:   R_N(t) - C_N(t) >= -8 (3 + sqrt N) S_H = -{8 * s:.4e} (3 + sqrt N)", log)
    escribir("        N        E_N simple     E_N con region sin ceros", log)
    exponentes = np.arange(2, 41, 2)
    simple = []
    region = []
    for k in exponentes:
        N = 10.0 ** k
        simple.append(E_simple(N))
        region.append(E_region(N))
        escribir(f"   10^{k:<3d}   {simple[-1]:.4e}     {region[-1]:.4e}", log)
    # longitud hasta la que la cota es menor que el margen medido y que el medio
    for nombre, margen in (("medido minimo", MARGEN_MEDIDO), ("medio en t = 0", MARGEN_MEDIO)):
        N_estrella = ((margen / (8 * s)) - 3) ** 2
        escribir(f"E_N simple < margen {nombre} ({margen}) mientras N < {N_estrella:.3e}", log)
    logN = np.linspace(np.log(10.0), np.log(1e40), 400)
    Es = np.array([E_simple(np.exp(x)) for x in logN])
    Er = np.array([E_region(np.exp(x)) for x in logN])
    i = int(np.argmax(Er > MARGEN_MEDIDO))
    escribir(f"E_N con region sin ceros < margen medido minimo mientras N < {np.exp(logN[i]):.3e}", log)

    fig, eje = plt.subplots(figsize=(10, 5.5))
    eje.loglog(np.exp(logN), Es, color="tab:red", label=r"$E_N=8(3+\sqrt{N})\,S_H$ (unconditional)")
    eje.loglog(np.exp(logN), Er, "--", color="tab:red", label="with the zero-free region of Mossinghoff-Trudgian-Yang")
    eje.axhline(MARGEN_MEDIDO, color="tab:blue", lw=1.0, label=r"least Fejer margin of $\zeta$ measured (up to $N=e^{12}$, $t<120$)")
    eje.axhline(MARGEN_MEDIO, color="tab:blue", ls=":", lw=1.0, label=r"mean margin at $t=0$ under RH, $4\sum\gamma^{-2}$")
    eje.axvspan(1.0, np.exp(12.0), color="0.9", zorder=0, label=r"computed range $N\leq e^{12}$")
    eje.set_xlabel("$N$")
    eje.set_ylabel("size")
    eje.set_title(r"Unconditional bound for the off-line contribution, against the margin of $\zeta$")
    eje.legend(fontsize=8, loc="upper left")
    fig.tight_layout()
    fig.savefig(os.path.join(CARPETA, "cota_incondicional.png"), dpi=120)
    log.close()
    plt.show()


if __name__ == "__main__":
    main()