""" Certificado: 0 < M(N) < 0.14 para todo N real con 1 < N <= X0 = 10^8, sin suponer RH (Parte 20, Corolario 'sandwich'). M(N) = R(N) - C(N), C(N) = 2 sum_{n=1} N^{-a_m}/a_m^2, a_m = 2m + 1/2, c_1 = log pi + gamma_E + pi/2 + 3 log 2, c_0 = -16 + pi^2/2 + 4G. Como 0.14 < 4 (log xi)''(1/2) = 0.18484, el certificado da el sandwich completo en 1 < N <= 10^8. Aritmetica. Todo se hace con el modulo decimal de Python a 34 cifras: suma, producto, division, sqrt, ln y exp estan correctamente redondeados (especificacion General Decimal Arithmetic, libmpdec), con error relativo <= 5e-34 por operacion. Las sumas acumuladas A y B tienen como mucho 6e6 terminos positivos, asi que su error relativo es < 1e-26, y el de cada valor de M (magnitudes < 1e6) es < 1e-19. Las constantes pi, gamma_E y G se escriben con 49 cifras exactas. Se resta (o suma) a cada cota un margen EPS = 1e-15, mucho mayor que todos esos errores. Tres tramos (L = log N, m(L) = M(e^L), r(L) = R(e^L)): I. 0 < L <= 1/15: lema analitico (demostrado en el texto), con las constantes evaluadas aqui: r(L) >= L (3.5 - c_1 + log 6.6) > 0 (usa 1 - e^{-x} >= x - x^2/2 para a_m L <= 1, y K >= 7 terminos), r(L) <= L (4 e^{1/60} - c_1 + 1.8 + log(0.4/L) + 2L) <= 0.102 para L <= 1/30. Como no hay potencias de primo en (1, 2), C = 0 y M = R en todo 1 < N < 2. II. 1/30 <= L <= log 2: subintervalos de anchura 2e-4. En [L1, L2], por monotonia de cada pieza, r >= 8 e^{L1/2} - c_1 L2 + c_0 - 2 T_sup(L1), r <= 8 e^{L2/2} - c_1 L1 + c_0 - 2 T_inf(L2), con T_inf la suma parcial y T_sup = T_inf + cota de la cola N^{-a_{K+1}}/(a_{K+1}^2 (1 - N^{-2})). III. 2 <= N <= X0: entre potencias de primo consecutivas q < q', m(L) = r(L) - 2 A L + 2 B con A, B constantes, y m'' = r'' = 2 sqrt N (1 - 1/(N (N^2 - 1))) > 0: m es convexa. Por eso, en [log q, log q'], min m >= max( m(q) + min(m'(q+), 0) D, m(q') - max(m'(q'-), 0) D ), D = log(q'/q), max m = max( m(q), m(q') ), con m'(L) = r'(L) - 2A, r'(L) = 4 e^{L/2} - c_1 + 2 S(L), S(L) = sum_m N^{-a_m}/a_m. Si la cota inferior no basta se subdivide el tramo (solo hace falta para N pequeno). Salida en resultados_certificado/: certificado_sandwich.log """ import os import sys import time import numpy as np from decimal import Decimal, getcontext BASE = os.path.dirname(os.path.abspath(__file__)) sys.path.insert(0, BASE) from paso2_bloques import potencias_de_primos CARPETA = os.path.join(BASE, "resultados_certificado") getcontext().prec = 34 # X0 se puede cambiar desde la linea de ordenes (por ejemplo para una prueba rapida con 10^6) X0 = int(float(sys.argv[1])) if len(sys.argv) > 1 else 10 ** 8 TECHO = Decimal("0.14") EPS = Decimal("1e-15") PI = Decimal("3.141592653589793238462643383279502884197169399375") EULER = Decimal("0.577215664901532860606512090082402431042159335939") CATALAN = Decimal("0.915965594177219015054603514932384110774149374281") LN2 = Decimal(2).ln() C1 = PI.ln() + EULER + PI / 2 + 3 * LN2 C0 = Decimal(-16) + PI * PI / 2 + 4 * CATALAN UNO = Decimal(1) DOS = Decimal(2) def escribir(texto, archivo): """imprime y guarda en el log una linea.""" print(texto, flush=True) archivo.write(texto + "\n") archivo.flush() def series(L, K): """(T_inf, T_sup, S_inf, S_sup) en N = e^L con K terminos y la cota de la cola: sum_{m>K} e^{-a_m L}/a_m^j <= e^{-a_{K+1} L}/(a_{K+1}^j (1 - e^{-2L})).""" T = Decimal(0) S = Decimal(0) for m in range(1, K + 1): a = Decimal(2 * m) + Decimal("0.5") e = (-a * L).exp() T = T + e / (a * a) S = S + e / a a = Decimal(2 * (K + 1)) + Decimal("0.5") cola = (-a * L).exp() / (UNO - (-2 * L).exp()) return T, T + cola / (a * a), S, S + cola / a def series_grandes(raiz_q): """para N >= 1000: T en [0, N^{-5/2}/(6.25 (1 - N^{-2}))], S en [0, N^{-5/2}/(2.5 (1 - N^{-2}))].""" N = raiz_q * raiz_q base = UNO / (N * N * raiz_q * (UNO - UNO / (N * N))) return Decimal(0), base / Decimal("6.25"), Decimal(0), base / Decimal("2.5") def tramo_I(log): """las dos desigualdades del lema para 0 < L <= 1/15 (cota inferior) y 0 < L <= 1/30 (cota superior).""" abajo = Decimal("3.5") - C1 + Decimal("6.6").ln() L0 = UNO / 30 arriba = L0 * (4 * (L0 / 2).exp() - C1 + Decimal("1.8") + (Decimal("0.4") / L0).ln() + 2 * L0) escribir(f"I. 0 < log N <= 1/15: R(N) >= {abajo:.6f} log N > 0, and for 0 < log N <= 1/30: R(N) <= {arriba:.6f} < {TECHO}", log) return abajo > 0 and arriba < TECHO def tramo_II(log): """1/30 <= L <= log 2, subintervalos de anchura 2e-4: cotas por monotonia de cada pieza de r.""" paso = Decimal("0.0002") L = UNO / 30 puntos = [] while L < LN2: puntos.append(L) L = L + paso puntos.append(LN2) valores = [] for L in puntos: # K con a_{K+1} L >= 60: la cola es menor que e^{-60} K = int(60 / (2 * float(L))) + 2 T_inf, T_sup, _, _ = series(L, K) valores.append((L, (L / 2).exp(), T_inf, T_sup)) minimo = None maximo = None for (L1, e1, Ti1, Ts1), (L2, e2, Ti2, Ts2) in zip(valores[:-1], valores[1:]): abajo = 8 * e1 - C1 * L2 + C0 - 2 * Ts1 - EPS arriba = 8 * e2 - C1 * L1 + C0 - 2 * Ti2 + EPS if minimo is None or abajo < minimo[0]: minimo = (abajo, L1) if maximo is None or arriba > maximo[0]: maximo = (arriba, L2) escribir(f"II. 1/30 <= log N <= log 2 ({len(puntos) - 1} subintervals): R(N) >= {minimo[0]:.6f} (near N = {float(minimo[1].exp()):.4f}), R(N) <= {maximo[0]:.6f} (near N = {float(maximo[1].exp()):.4f})", log) return minimo[0] > 0 and maximo[0] < TECHO, minimo[0], maximo[0] def valores_en(L, raiz, A, B, pequeno): """(m_inf, m_sup, r'_inf, r'_sup) en el punto L (N = e^L, raiz = sqrt N) con sumas acumuladas A, B.""" if pequeno: T_inf, T_sup, S_inf, S_sup = series(L, 30) else: T_inf, T_sup, S_inf, S_sup = series_grandes(raiz) base = 8 * raiz - C1 * L + C0 - 2 * A * L + 2 * B deriv = 4 * raiz - C1 return base - 2 * T_sup - EPS, base - 2 * T_inf + EPS, deriv + 2 * S_inf - EPS, deriv + 2 * S_sup + EPS def cota_tramo(La, va, Lb, vb, A): """cota inferior del minimo de la funcion convexa m en [La, Lb] a partir de los valores y pendientes en los extremos.""" D = Lb - La m_a_inf, _, rp_a_inf, _ = va m_b_inf, _, _, rp_b_sup = vb pend_a = rp_a_inf - 2 * A pend_b = rp_b_sup - 2 * A uno = m_a_inf + min(pend_a, Decimal(0)) * D otro = m_b_inf - max(pend_b, Decimal(0)) * D return max(uno, otro) def subdividir(La, Lb, A, B, piezas): """cota inferior en [La, Lb] partiendo en 'piezas' subtramos (N pequeno: series con 30 terminos).""" paso = (Lb - La) / piezas puntos = [La + paso * j for j in range(piezas)] + [Lb] vals = [valores_en(L, (L / 2).exp(), A, B, True) for L in puntos] return min(cota_tramo(puntos[j], vals[j], puntos[j + 1], vals[j + 1], A) for j in range(piezas)) def tramo_III(log): """2 <= N <= X0 por convexidad entre potencias de primo consecutivas.""" t0 = time.time() n, lp, k = potencias_de_primos(X0) escribir(f"III. 2 <= N <= {X0:.0e}: {n.size} prime powers", log) ln_primo = {} A = Decimal(0) B = Decimal(0) minimo = (None, None) minimo_grande = (None, None) maximo = (None, None) subdivididos = 0 previo = None # recorre las potencias de primo en orden, y el tramo [q_i, q_{i+1}] usa las sumas A_i, B_i (hasta q_i incluida) lista = [(int(q), int(e)) for q, e in zip(n, k)] + [(X0, 0)] for q, e in lista: if e == 0: Lq = Decimal(q).ln() raiz = Decimal(q).sqrt() else: if e == 1: p = q else: p = int(round(q ** (1.0 / e))) while p ** e > q: p = p - 1 while (p + 1) ** e <= q: p = p + 1 if p not in ln_primo: ln_primo[p] = Decimal(p).ln() if e == 1 and p > 10 ** 4: # los primos grandes no vuelven a aparecer como base de una potencia: no se guardan lnp = ln_primo.pop(p) else: lnp = ln_primo[p] else: lnp = ln_primo[p] Lq = e * lnp raiz = Decimal(q).sqrt() pequeno = q < 1000 # valores en q con las sumas del tramo anterior (A, B hasta la potencia anterior): extremo derecho de ese tramo v_izq = valores_en(Lq, raiz, A, B, pequeno) if previo is not None: Lp, vp, Ap = previo cota = cota_tramo(Lp, vp, Lq, v_izq, Ap) # si la cota es pobre (tramos largos, N pequeno) se subdivide, primero en 64 y luego en 2048 piezas if cota < Decimal("0.02"): cota = max(cota, subdividir(Lp, Lq, Ap, B, 64)) subdivididos = subdivididos + 1 if cota < Decimal("0.02"): cota = max(cota, subdividir(Lp, Lq, Ap, B, 2048)) if minimo[0] is None or cota < minimo[0]: minimo = (cota, (Lp + Lq) / 2) # minimo aparte para N >= 1000, donde los tramos son cortos y la cota por tangentes es casi exacta if q >= 1000 and (minimo_grande[0] is None or cota < minimo_grande[0]): minimo_grande = (cota, (Lp + Lq) / 2) if maximo[0] is None or v_izq[1] > maximo[0]: maximo = (v_izq[1], Lq) if e == 0: break # se incorpora q a las sumas: A += Lambda(q)/sqrt q, B += Lambda(q) log q/sqrt q termino = lnp / raiz A = A + termino B = B + termino * Lq previo = (Lq, valores_en(Lq, raiz, A, B, pequeno), A) escribir(f" least lower bound for M on an interval: {minimo[0]:.6f} (around N = {float(minimo[1].exp()):.5g}), largest upper bound at a prime power: {maximo[0]:.6f} (N = {float(maximo[1].exp()):.5g}), intervals subdivided: {subdivididos}, {time.time() - t0:.0f} s", log) escribir(f" least lower bound for 1000 <= N <= {X0:.0e}: {minimo_grande[0]:.6f} (around N = {float(minimo_grande[1].exp()):.5g})", log) return minimo[0] > 0 and maximo[0] < TECHO, minimo[0], maximo[0] def main(): os.makedirs(CARPETA, exist_ok=True) log = open(os.path.join(CARPETA, "certificado_sandwich.log"), "w") escribir(f"Decimal arithmetic, {getcontext().prec} digits. c_1 = {C1:.20f}, c_0 = {C0:.20f}", log) ok1 = tramo_I(log) ok2, min2, max2 = tramo_II(log) ok3, min3, max3 = tramo_III(log) escribir("", log) if ok1 and ok2 and ok3: escribir(f"CERTIFIED: 0 < M(N) < {TECHO} for every real N with 1 < N <= {X0:.0e} (and {TECHO} < 4 (log xi)''(1/2) = 0.18484),", log) escribir(f" with M(N) >= {min(min2, min3):.6f} for e^(1/30) = 1.0339 <= N <= {X0:.0e} (below that M = R > 0 by the lemma of I)", log) escribir(f" and M(N) <= {max(max2, max3):.6f} for e^(1/30) <= N <= {X0:.0e} (below that M = R <= 0.102)", log) else: escribir(f"NOT certified: I {ok1}, II {ok2}, III {ok3}", log) log.close() if __name__ == "__main__": main()