"""cert_qrh_range.py: the sieve-range margin of the paper's Lemma 20 (under the 7/8 hypothesis) on L = log N in [L_lo, L_hi]. Its margin_iv is imported by second/cert_qrhB3.py and second/cert_qrhB3p.py, which call it with K = M - 1, so that the term bad below is 0, and then subtract their own bad count. Per even N with L = log N, parameters (M, lam, t, K): z = N^{1/M}, H = lam N^{1/8}/L^3 (sieve level QD <= H, Q <= 1155), D = H/1155, s = log D/log z. main = (1/L)(U M/L)(1 - 1.002 M/(z-1))(1 - 1/(2 log^2 z)) F(s), F(s) = f(s) - eps C2 e^2 h(s), eps = 4.77e-3, C2 = 154 [Johnston Thm 2.1, L4.2, L4.3] rem = 0.45 (H/N) [ (cE/L)(1 + 32/(7L)) N^{7/8} + 3 N^{1/2} + 2 log H/log 3 + 1.3841 L/log L ] [L4.4, Proposition 18] one = 1/N (the element N - p = 1) bad = (1 + 1.2762 M/L)(M/L) t^{-K} exp( t (log(M-1) + D(L(M-1)/M) + D(L/M) + e^{-L/M}) + t^2 e^{-L/M}/(2(1 - t e^{-L/M})(1 - e^{-L/M})) ) [a Rankin-type bad count, not used by the paper] need main - rem - one - bad > 0. On an interval [L1, L2] with fixed parameters every term is bounded by its worst value: s, z, F, the two (1 - ...) factors at L1; the factor M U/L^2 at L2; in rem, cE at L2 (cE increases with L and with log H) and the powers of L at L1 (H/N <= lam e^{-7 L1/8}/L1^3); bad at L1 (each factor decreases with L). mode 'float' searches parameters; mode 'iv' certifies with mpmath interval arithmetic. """ import math, sys, json from fractions import Fraction GAMMA = 0.5772156649015329 A_, B_ = 0.298, 4.358 EPS, C2 = 4.77e-3, 154 U_LOW = 0.7413082243 # ---------------- float version (search) ---------------- def f_lin(s): return 2 * math.exp(GAMMA) * math.log(s - 1) / s def h_fn(s): return math.exp(-s) # h(s) <= e^{-s} on [2, 4] (h = e^{-s} on [2,3], 3e^{-s}/s <= e^{-s} on [3,4]) def Fs(s): return f_lin(s) - EPS * C2 * math.e ** 2 * h_fn(s) C2PIE = math.log(2 * math.pi * math.e) def Nplus1(lq): return (lq - C2PIE) / math.pi + A_ * lq + B_ def Z1inc(lq, lT): iT = math.exp(-lT) return (7 * Nplus1(lq) + ((lq - C2PIE) * lT + lT ** 2 / 2) / math.pi + (A_ * lq + B_) * (1 - iT) + A_ * (1 - (lT + 1) * iT) + (4 / math.pi) * (lq + lT - C2PIE + 1)) def Z2(lq): return 6.12 * Nplus1(lq) + (2 / math.pi) * (lq - C2PIE + 1) + A_ * lq + A_ / 2 + B_ def l1pe(x): return x + math.log1p(math.exp(-x)) if x > 30 else math.log1p(math.exp(x)) def cE(L, lH): T1 = (16 / 7) * (1 + math.exp(-L / 16) / (L / 2)) T2 = (L + 1) * (0.65 * math.exp(-L / 4) + (math.log((L + 1) / math.log(2)) + 1) * math.exp(-7 * L / 16)) T3 = 0.5 + 1.5 * math.exp(-7 * L / 16) T4 = (1 + math.exp(-L / 16)) ** (7 / 8) * (Z1inc(lH, l1pe(L / 8)) + 2 * (A_ * lH + A_ / 2 + B_) * math.exp(-l1pe(L / 16)) + 2 * A_ / math.e) T5 = (4.67 * Z2(lH) + L + 1 + 3.1) * math.exp(-7 * L / 16) return T1 + T2 + T3 + T4 + T5 def Dlog(l): return 1 / (10 * l * l) + 4 / (15 * l ** 3) def margin_float(L1, L2, M, lam, t, K): lH1 = L1 / 8 + math.log(lam) - 3 * math.log(L1) lH2 = L2 / 8 + math.log(lam) - 3 * math.log(L2) lz = L1 / M s = (lH1 - math.log(1155)) / lz if s <= 2.0 or s > 4: return None z = math.exp(min(lz, 700.0)) if z < 694 or t >= z: return None F = Fs(s) if F <= 0: return None main = U_LOW * M / L2 ** 2 * (1 - 1.002 * M / (z - 1)) * (1 - 1 / (2 * lz ** 2)) * F # rem / N: (H/N) N^{7/8} = lam/L^3 exactly; bounded by lam/L1^3 rem = 0.45 * lam / L1 ** 3 * (cE(L2, max(lH1, lH2)) / L1 * (1 + 32 / (7 * L1)) + 3 * math.exp(-3 * L1 / 8) + (2 * max(lH1, lH2) / math.log(3) + 1.3841 * L2 / math.log(L2)) * math.exp(-7 * L1 / 8)) one = math.exp(-L1) if K >= M - 1: return main - rem - one, main, rem, 0.0, s # z-rough n < N = z^M have Omega(n) <= M - 1: no bad n bad = (1 + 1.2762 * M / L1) * (M / L1) * math.exp(t * (math.log(M - 1) + Dlog(L1 * (M - 1) / M) + Dlog(lz) + 1 / z) + t * t / (2 * (1 - t / z) * (z - 1)) - K * math.log(t)) return main - rem - one - bad, main, rem, bad, s def search(L1, L2, K): best = None for M in range(19, 70): for lam in (0.1, 1, 10, 100, 1000): tt0 = K / math.log(M - 1) for x in (0.85, 0.9, 0.95, 1.0, 1.05, 1.1): t = round(tt0 * x, 2) if t < 1: continue r = margin_float(L1, L2, M, lam, t, K) if r and r[0] > 0: rel = r[0] / r[1] if best is None or rel > best[0]: best = (rel, M, lam, t) return best # ---------------- interval version (certificate) ---------------- def margin_iv(L1, L2, M, lam, t, K): from mpmath import iv iv.prec = 120 L1 = iv.mpf(L1); L2 = iv.mpf(L2) lam = iv.mpf(lam); tt = iv.mpf(t.numerator) / t.denominator if isinstance(t, Fraction) else iv.mpf(t) ln = iv.log g = iv.euler a, b = iv.mpf('0.298'), iv.mpf('4.358') c2pie = ln(2 * iv.pi * iv.e) def Np1(lq): return (lq - c2pie) / iv.pi + a * lq + b def Z1i(lq, lT): iT = iv.exp(-lT) return (7 * Np1(lq) + ((lq - c2pie) * lT + lT ** 2 / 2) / iv.pi + (a * lq + b) * (1 - iT) + a * (1 - (lT + 1) * iT) + (4 / iv.pi) * (lq + lT - c2pie + 1)) def Z2i(lq): return iv.mpf('6.12') * Np1(lq) + (2 / iv.pi) * (lq - c2pie + 1) + a * lq + a / 2 + b def l1pe_i(x): return x + ln(1 + iv.exp(-x)) def cE_i(Ld, Li, lH): # sup over L in [Ld, Li]: terms decreasing in L are evaluated at Ld, the increasing Z_1 part at Li (T2's factor L + 1 at Li) T1 = (iv.mpf(16) / 7) * (1 + iv.exp(-Ld / 16) / (Ld / 2)) T2 = (Li + 1) * (iv.mpf('0.65') * iv.exp(-Ld / 4) + (ln((Li + 1) / ln(iv.mpf(2))) + 1) * iv.exp(-7 * Ld / 16)) T3 = iv.mpf('0.5') + iv.mpf('1.5') * iv.exp(-7 * Ld / 16) T4 = (1 + iv.exp(-Ld / 16)) ** (iv.mpf(7) / 8) * (Z1i(lH, l1pe_i(Li / 8)) + 2 * (a * lH + a / 2 + b) * iv.exp(-l1pe_i(Ld / 16)) + 2 * a / iv.e) T5 = (iv.mpf('4.67') * Z2i(lH) + Li + 1 + iv.mpf('3.1')) * iv.exp(-7 * Ld / 16) return T1 + T2 + T3 + T4 + T5 lH1 = L1 / 8 + ln(lam) - 3 * ln(L1) lH2 = L2 / 8 + ln(lam) - 3 * ln(L2) lHmax = iv.mpf([min(lH1.a, lH2.a), max(lH1.b, lH2.b)]) lz = L1 / M s = (lH1 - ln(iv.mpf(1155))) / lz # explicit certificate gates if not (isinstance(M, int) and M >= 2): return None, 'gate: integer M >= 2' if not (isinstance(K, int) and K >= 1): return None, 'gate: integer K >= 1' if not (L1.a >= 300 and L2.a >= L1.b): return None, 'gate: L1 >= 300, L2 >= L1' if not (lam.a > 0 and lam.b <= 100): return None, 'gate: 0 < lam <= 100' if not (lH1.a >= ln(iv.mpf(25))): return None, 'gate: H(L1) >= 25' if not (lH2.b < L2.a and lH1.b < L1.a): return None, 'gate: H < N' if not (L1 * (M - 1) / M).a >= ln(iv.mpf(10372)): return None, 'gate: N/z >= 10372' s2 = (lH2 - ln(iv.mpf(1155))) / (L2 / M) if not (s.a > 2 and s.b <= 4 and s2.b <= 4): return None, 'gate: 2 < s1 <= s2 <= 4 (s increases with L)' z = iv.exp(lz) if not (z.a >= 694 and tt.b < z.a): return None, 'z too small' F = 2 * iv.exp(g) * ln(s - 1) / s - iv.mpf('4.77e-3') * 154 * iv.e ** 2 * iv.exp(-s) # h(s) <= e^{-s} on [2,4] if not (F.a > 0): return None, 'gate: F(s1) > 0' if not ((1 - iv.mpf('1.002') * M / (z - 1)).a > 0 and (1 - 1 / (2 * lz ** 2)).a > 0): return None, 'gate: V factors positive' U = 2 * iv.exp(-g) * iv.mpf('0.6601618158') main = U * M / L2 ** 2 * (1 - iv.mpf('1.002') * M / (z - 1)) * (1 - 1 / (2 * lz ** 2)) * F rem = iv.mpf('0.45') * lam / L1 ** 3 * (cE_i(L1, L2, lHmax) / L1 * (1 + iv.mpf(32) / (7 * L1)) + 3 * iv.exp(-3 * L1 / 8) + (2 * lHmax / ln(iv.mpf(3)) + iv.mpf('1.3841') * L2 / ln(L2)) * iv.exp(-7 * L1 / 8)) one = iv.exp(-L1) def Dl(l): return 1 / (10 * l * l) + 4 / (15 * l ** 3) if K >= M - 1: # z-rough n < N = z^M have Omega(n) <= M - 1, so there are no bad n m = main - rem - one return m, (main, rem, 0, s, F) bad = (1 + iv.mpf('1.2762') * M / L1) * (M / L1) * iv.exp(tt * (ln(iv.mpf(M - 1)) + Dl(L1 * (M - 1) / M) + Dl(lz) + 1 / z) + tt * tt / (2 * (1 - tt / z) * (z - 1)) - K * ln(tt)) m = main - rem - one - bad return m, (main, rem, bad, s, F) if __name__ == '__main__': mode = sys.argv[1] if len(sys.argv) > 1 else 'float' K = int(sys.argv[2]) if len(sys.argv) > 2 else 18 Lstart = float(sys.argv[3]) if len(sys.argv) > 3 else 300.0 Lend = float(sys.argv[4]) if len(sys.argv) > 4 else 51046.0 ratio = float(sys.argv[5]) if len(sys.argv) > 5 else 1.01 out = [] L1 = Lstart fails = 0 while L1 < Lend: L2 = min(Lend, L1 * ratio) b = search(L1, L2, K) if b is None: print(f'FAIL float search on [{L1:.3f}, {L2:.3f}] for K = {K}'); fails += 1 out.append((L1, L2, None)) else: rel, M, lam, t = b row = [L1, L2, M, lam, t, rel] if mode == 'iv': tf = Fraction(t).limit_denominator(1000) m, info = margin_iv(L1, L2, M, lam, tf, K) ok = (m is not None) and m.a > 0 row += [str(tf), bool(ok)] if not ok: print(f'IV FAIL on [{L1:.3f}, {L2:.3f}]', b, info); fails += 1 out.append(row) L1 = L2 # coverage: consecutive intervals share endpoints exactly, the first starts at Lstart and the last ends at Lend assert out[0][0] == Lstart and all(out[i][1] == out[i + 1][0] for i in range(len(out) - 1)) and out[-1][1] == Lend, 'coverage' print(f'K = {K}, L in [{Lstart}, {Lend}], ratio {ratio}: {len(out)} intervals, failures: {fails}; coverage and all gates checked in iv mode' if mode == 'iv' else f'K = {K}, L in [{Lstart}, {Lend}], ratio {ratio}: {len(out)} intervals, failures: {fails}') if out: print('first rows:', out[:3]); print('last rows:', out[-3:]) json.dump({'K': K, 'Lstart': Lstart, 'Lend': Lend, 'ratio': ratio, 'mode': mode, 'rows': out, 'failures': fails}, open(f'cert_qrh_range_K{K}_{mode}.json', 'w'))