#!/usr/bin/env python3 """c3_grh.py: the middle range under GRH (the paper's Lemmas 22 and 23): the condition of Lemma 6 with a GRH good count presieved to q_w = prod_{3 <= p <= w} p. Its row w = 11 is Lemma 6 itself (Lemma 7 at (e^62, 9)). Statement. Assume GRH. Let w >= 3, 4e18 <= Y_low <= X, L_low = log Y_low, L_X = log X, q_w^2 <= Y_low - 2, 0 < delta < 1, delta L_low >= 1, and K >= 1 an integer. Every even N in [Y_low, X] with N > 4e18 is p + eta with p an odd prime, eta >= 3, (eta, 2 q_w) = 1 and Omega(eta) <= K, provided L_low e^{-delta L_low} + (1-delta)^{-1} sum_{k=K+1}^{kmax} [c_th R*_{k-1} + sum_{a=2}^k A_a R*_{k-a}] < g_w (1 - 4 e^{-L_low}) - 0.64 phi(q_w) L_low^2 e^{-L_low/2} * 1.01 - L_low e^{-L_low} - 1e-12, g_w := prod_{3<=p<=w} (p-2)/(p-1), kmax := floor(L_X / log w+) (w+ the next prime: (n, 2 q_w) = 1 forces P^-(n) >= w+), R*_j := [z^j] exp(z S'_w(X) + C_w(z)), S'_w(X) := loglog X + B + D(log X) - sigma_w, sigma_w := sum_{p<=w} 1/p, C_w(z) := sum_{k>=2} P_w(k) z^k/k, P_w(k) := sum_{p>w} p^-k, A_a := sum_{p>w} a log p/p^a, c_th = 1 + 1/36260. Proof sketch. Good: by Bordignon-Starichkova's Lemma 3 (arXiv:2211.08844v1, from Ernvall-Hytonen-Palojarvi, Theorem 1), for 3 <= d <= sqrt(X2), N >= X2 >= 4e18 and (a, d) = 1: |pi(N; d, a) - li(N)/phi(d)| <= 0.640 sqrt(N) log N. Put y=N-2 and X_2=max(4e18,Y_low-2). For even N>4e18 in the stated range, y>=X_2. The hypothesis q_w^2<=Y_low-2 implies 3<=q_w<=sqrt(X_2), so the source's .640 bound applies. The theorem for all N uses (I1) at N=4e18; no presieved coprimality conclusion is claimed there from (I1). Alternatively, the unconditional Lemma 6 covers the entire first piece through e^62. With d = q_w: the classes a mod q_w with (a, q_w) = 1 = (N - a, q_w) number #Acal >= prod (p-2) = g_w phi(q_w), so #{odd p <= y : p mod q_w in Acal} >= (#Acal/phi(q_w)) li(y) - #Acal 0.64 sqrt(y) log y - 1 >= g_w li(y) - phi(q_w) 0.64 sqrt(y) log y - 1 (#Acal <= phi(q_w)). Each gives eta = N - p >= 3 odd with (eta, q_w) = 1. Bad: E := {n <= N : (n, 2 q_w) = 1, Omega(n) >= K+1}: #E <= N^{1-delta} + (N/((1-delta)L)) sum_k [...] by Hardy-Ramanujan on n > N^{1-delta} (log n > (1-delta)L; a = 1 terms via theta(N/n') <= c_th N/n' (Dusart), a >= 2 via A_a), and rho_j <= R*_j (Lemma 4 with the primes <= w removed; (I3) at X >= 10372). Normalise by N/L: L li(y)/N >= 1 - 4/N for y = N - 2 >= 4e18 (li(y) >= (y-2)/log y, log y <= L), and the GRH error phi(q_w) 0.64 L sqrt(y) log y/N <= 0.64 phi(q_w) L^2 e^{-L/2}, decreasing in L (L > 4), so take it at L_low (factor 1.01 covers the 1 - 4/N and L/N slack); L e^{-delta L} <= L_low e^{-delta L_low}. QED (sketch).""" import os, sys, math, time from mpmath import iv, mp iv.dps = 30; mp.dps = 30 def I(x): return iv.mpf(x) def lo(x): return mp.make_mpf(x._mpi_[0]) def hi(x): return mp.make_mpf(x._mpi_[1]) import numpy as np NP = 200000 sv = np.ones(NP + 1, dtype=bool); sv[:2] = False for i in range(2, int(NP ** 0.5) + 1): if sv[i]: sv[i * i::i] = False ALLP = [int(p) for p in np.nonzero(sv)[0]] B = I('0.26149721284764278375542683860869585905156664826120'); B = iv.mpf([lo(B) - mp.mpf('1e-28'), hi(B) + mp.mpf('1e-28')]) CTH = 1 + I(1) / 36260 def D(l): return 1 / (10 * l ** 2) + I(4) / (15 * l ** 3) def cut(k): return 200000 if k <= 6 else (20000 if k <= 20 else 3000) class W: def __init__(s, w): s.w = w; s.PR = [p for p in ALLP if p > w]; s.wplus = s.PR[0] s.sigma = sum((I(1) / p for p in ALLP if p <= w), I(0)) s.q = 1; s.phi = 1; s.g = I(1) for p in ALLP: if 3 <= p <= w: s.q *= p; s.phi *= (p - 1); s.g *= I(p - 2) / (p - 1) s._P = {}; s._A = {} def P(s, k): if k not in s._P: Pc = cut(k); tot = I(0) for p in s.PR: if p > Pc: break tot += I(1) / I(p) ** k s._P[k] = iv.mpf([lo(tot), hi(tot + I(Pc) ** (1 - k) / (k - 1))]) return s._P[k] def A(s, a): if a not in s._A: Pc = cut(a); tot = I(0) for p in s.PR: if p > Pc: break tot += iv.log(I(p)) / I(p) ** a lP = iv.log(I(Pc)); tail = lP / I(Pc) ** a + I(Pc) ** (1 - a) * (lP / (a - 1) + 1 / I(a - 1) ** 2) s._A[a] = a * iv.mpf([lo(tot), hi(tot + tail)]) return s._A[a] def lhs_rhs(w, LX, Llow, K, delta): S = W(w) LX = I(LX); Llow = I(Llow); delta = I(delta) assert lo(delta * Llow) >= 1 assert S.q ** 2 <= lo(iv.exp(Llow)) - 2, 'q_w <= sqrt(Y_low - 2)' q = LX / iv.log(I(S.wplus)); kmax = int(mp.floor(lo(q))); assert int(mp.floor(hi(q))) == kmax Sp = iv.log(LX) + B + D(LX) - S.sigma e = [I(1)] for n in range(1, kmax + 1): acc = Sp * e[n - 1] for k in range(2, n + 1): acc += S.P(k) * e[n - k] e.append(acc / n) tot = I(0) for k in range(K + 1, kmax + 1): inner = CTH * e[k - 1] for a in range(2, k + 1): inner += S.A(a) * e[k - a] tot += inner lhs = Llow * iv.exp(-delta * Llow) + tot / (1 - delta) rhs = S.g * (1 - 4 * iv.exp(-Llow)) - I('0.64') * S.phi * Llow ** 2 * iv.exp(-Llow / 2) * I('1.01') - Llow * iv.exp(-Llow) - I('1e-12') return lhs, rhs, S, kmax if __name__ == '__main__': cases = [(11, '58', 'log4e18', 9), (11, '62', 'log4e18', 9), (11, '65', 'log4e18', 9)] for (w, LX, Llow, K) in [(11, '62', None, 9), (13, '65', '60', 9), (17, '68', '62', 9), (19, '70', '62', 9), (23, '73', '65', 9), (29, '73', '65', 9), (31, '73', '67', 9), (29, '73', '65', 8), (31, '73', '67', 8)]: Ll = iv.log(I('4e18')) if Llow is None else I(Llow) best = None for d in ('0.10', '0.15', '0.20', '0.25', '0.30'): if lo(I(d) * Ll) < 1: continue try: l, r, S, kmax = lhs_rhs(w, LX, Ll, K, d) except AssertionError as ex: print('w = %d, [e^%s, e^%s]: gate fails (%s)' % (w, Ll, LX, ex)); best = None; break if best is None or hi(l) < hi(best[0]): best = (l, r, d, S, kmax) if best is None: continue l, r, d, S, kmax = best print(f"w={w}, q_w={S.q}, log-range=[{Ll},{LX}], K={K}, delta={d}: " f"left-side enclosure {l}; right-side enclosure {r} -> " f"{'CERTIFIED' if hi(l) < lo(r) else 'NOT ESTABLISHED'}") sys.stdout.flush()