"""
El sandwich por los ceros en 10^8 <= N <= 10^15 (Parte 20, Proposicion 'uncond' y Corolario 'sandwich').

El argumento.  Para todo N > 1, sin suponer RH (Proposicion 'fejer' con t = 0):  M(N) = sum_rho F_L(gamma_rho), L = log N,
F_L(z) = 4 sin^2(L z/2)/z^2.  Por Platt y Trudgian, todos los ceros con 0 < Im rho <= H = 3*10^12 estan en la recta, y
cada uno aporta F_L(gamma) >= 0 (con su conjugado, 4 (1 - cos(gamma L))/gamma^2).  Los ceros por encima de H aportan en
valor absoluto como mucho 8 (3 + sqrt N) S_H, S_H <= 2.854*10^-12 (Proposicion 'uncond').  Por tanto, para cualquier
lista de ceros en la recta por debajo de H (no hace falta que este completa: los ceros que falten solo suman terminos >= 0),
    M(N) >= W(L) - 8 (3 + sqrt N) S_H,      W(L) = 4 sum_{gamma en la lista} (1 - cos(gamma L))/gamma^2.
W es lipschitziana con constante Lip = 4 sum 1/gamma, asi que en cada celda [L_j, L_j + h] de una malla
    W >= min(W(L_j), W(L_j + h)) - Lip h/2.
Entrada: los 4520 ceros por debajo de 5000, calculados como cambios de signo de la funcion Z de Hardy y refinados
(paso2_bloques.ceros_por_trozos), comprobados contra los valores publicados de los cinco primeros.  Este es el unico
dato que no esta certificado con intervalos: el resultado depende de la exactitud de la lista (un error delta en cada
cero cambia W en menos de 4 sum L/gamma^2 delta, unos 3.2 delta en log N = 35).
Salida en resultados_certificado/: certificado_ceros.log
"""
import os
import sys
import time
import numpy as np

BASE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, BASE)
from paso2_bloques import ceros_por_trozos

CARPETA = os.path.join(BASE, "resultados_certificado")
N_INI = 1e8
N_FIN = 1e15
T0 = 5000.0
PASO = 2e-5
S_H = 2.854e-12
PUBLICADOS = np.array([14.134725141734693790, 21.022039638771554993, 25.010857580145688763, 30.424876125859513210, 32.935061587739189691])


def escribir(texto, archivo):
    """imprime y guarda en el log una linea."""
    print(texto, flush=True)
    archivo.write(texto + "\n")
    archivo.flush()


def main():
    os.makedirs(CARPETA, exist_ok=True)
    log = open(os.path.join(CARPETA, "certificado_ceros.log"), "w")
    t0 = time.time()
    gam = ceros_por_trozos(T0)
    error_publicados = float(np.max(np.abs(gam[:5] - PUBLICADOS)))
    escribir(f"{gam.size} zeros of zeta below {T0:.0f} (sign changes of Z),  largest difference with the published values of the first five: {error_publicados:.1e}", log)
    lip = 4 * float(np.sum(1 / gam))
    a = 4 / gam ** 2
    escribir(f"Lipschitz constant of W: 4 sum 1/gamma = {lip:.4f},  grid step {PASO} in log N,  loss per cell Lip h/2 = {lip * PASO / 2:.2e}", log)
    # la malla empieza un paso antes de log 10^8 (el flotante log(1e8) supera al valor exacto en ~1e-15) y acaba despues de log 10^15
    L = np.arange(np.log(N_INI) - PASO, np.log(N_FIN) + PASO, PASO)
    W = np.empty(L.size)
    media = float(np.sum(a))
    for i0 in range(0, L.size, 4000):
        bloque = L[i0:i0 + 4000]
        W[i0:i0 + 4000] = media - np.cos(np.outer(bloque, gam)) @ a
    celda = np.minimum(W[:-1], W[1:]) - lip * PASO / 2
    error_cola = 8 * (3 + np.sqrt(np.exp(L[1:]))) * S_H
    cota = celda - error_cola
    j = int(np.argmin(cota))
    escribir(f"log N in [{L[0]:.4f}, {L[-1]:.4f}], {L.size} grid points,  {time.time() - t0:.0f} s", log)
    escribir(f"least value of W on the grid: {W.min():.6f} at N = {np.exp(L[np.argmin(W)]):.4e}", log)
    escribir(f"off-line bound 8 (3 + sqrt N) S_H at N = 10^15: {8 * (3 + np.sqrt(N_FIN)) * S_H:.2e}", log)
    escribir(f"lower bound for M(N), 10^8 <= N <= 10^15: {cota.min():.6f} (cell at N = {np.exp(L[j]):.4e})", log)
    for x in (1e9, 1e10, 1e11, 1e12, 1e13, 1e14, 1e15):
        sel = L[1:] <= np.log(x)
        escribir(f"   least lower bound up to N = {x:.0e}: {cota[sel].min():.6f}", log)
    if cota.min() > 0:
        escribir(f"VERIFIED (given the zero list): M(N) > 0 for every real N with 10^8 <= N <= 10^15", log)
    log.close()


if __name__ == "__main__":
    main()