collatz_reverse_add_audit.py

#!/usr/bin/env python3
"""Independent audit of the Collatz guard and cross-length reverse-add units.

The implementation is dependency-free and does not import the archived packet
checkers.  All arithmetic is exact.
"""

from __future__ import annotations

from collections import Counter, defaultdict
from fractions import Fraction
from itertools import product
import argparse
import json
import math


def bits(k):
    return product((0, 1), repeat=k)


def collatz_block(word):
    a = 0
    ones = 0
    for j, bit in enumerate(word):
        if bit:
            a = 3 * a + (1 << j)
            ones += 1
    return ones, a


def collatz_residue(word):
    k = len(word)
    if not k:
        return 0
    ones, a = collatz_block(word)
    modulus = 1 << k
    return (-a * pow(pow(3, ones, modulus), -1, modulus)) % modulus


def itinerary(n, k):
    out = []
    for _ in range(k):
        out.append(n & 1)
        n = n // 2 if n % 2 == 0 else (3 * n + 1) // 2
    return tuple(out), n


def primitive_necklace(word):
    k = len(word)
    if any(k % d == 0 and word == word[:d] * (k // d) for d in range(1, k)):
        return False
    return word == min(word[i:] + word[:i] for i in range(k))


def selector_values(depth, mask):
    vertices = list(bits(depth))
    names = {v: i for i, v in enumerate(vertices)}
    successor = {
        v: v[1:] + ((mask >> names[v]) & 1,) for v in vertices
    }
    values = {}
    for root in vertices:
        if root in values:
            continue
        path = []
        position = {}
        v = root
        while v not in values and v not in position:
            position[v] = len(path)
            path.append(v)
            v = successor[v]
        if v in position:
            cycle = path[position[v]:]
            word = tuple(x[0] for x in cycle)
            ones, a = collatz_block(word)
            values[cycle[0]] = Fraction(a, (1 << len(word)) - 3**ones)
            for x in cycle[:-1]:
                y = values[x]
                values[successor[x]] = y / 2 if x[0] == 0 else (3 * y + 1) / 2
            x = cycle[-1]
            y = values[x]
            assert (y / 2 if x[0] == 0 else (3 * y + 1) / 2) == values[cycle[0]]
        for x in reversed(path):
            if x in values:
                continue
            y = values[successor[x]]
            values[x] = 2 * y if x[0] == 0 else (2 * y - 1) / 3
    for v, value in values.items():
        assert value.denominator % 2 == 1
        modulus = 1 << depth
        observed = value.numerator * pow(value.denominator, -1, modulus) % modulus
        assert observed == collatz_residue(v)
    return values


def collatz_audit():
    word_count = 0
    seam_count = 0
    swap_count = 0
    necklaces = 0
    integral_necklaces = 0
    for k in range(1, 17):
        seen = set()
        for word in bits(k):
            word_count += 1
            ones, a = collatz_block(word)
            r = collatz_residue(word)
            observed, _ = itinerary(r, k)
            assert observed == word
            seen.add(r)
            if k <= 12:
                for split in range(k + 1):
                    u, v = word[:split], word[split:]
                    su, au = collatz_block(u)
                    sv, av = collatz_block(v)
                    assert a == 3**sv * au + (1 << split) * av
                    seam_count += 1
            if k <= 14:
                modulus = 1 << k
                for i in range(k - 1):
                    if word[i:i + 2] != (0, 1):
                        continue
                    swapped = word[:i] + (1, 0) + word[i + 2:]
                    delta = (collatz_residue(swapped) - r) % modulus
                    assert delta and (delta & -delta).bit_length() - 1 == i
                    swap_count += 1
            if primitive_necklace(word):
                necklaces += 1
                denominator = (1 << k) - 3**ones
                if a % denominator == 0:
                    integral_necklaces += 1
        assert seen == set(range(1 << k))

    sharp_response_pairs = 0
    for h in range(1, 9):
        for k in range(1, 9):
            n, m = 1 << (h + k - 1), 1 << (h + k)
            assert n % (1 << (h + k - 1)) == m % (1 << (h + k - 1))
            _, tn = itinerary(n, k)
            _, tm = itinerary(m, k)
            assert tn % (1 << h) != tm % (1 << h)
            sharp_response_pairs += 1

    assert collatz_block((0, 0, 1)) == (1, 4)
    assert Fraction(4, 8 - 3) == Fraction(4, 5)
    for k in range(1, 17):
        assert collatz_residue((1,) * k) == (1 << k) - 1

    selector_counts = {}
    for depth in range(1, 5):
        counts = Counter()
        for mask in range(1 << (1 << depth)):
            values = selector_values(depth, mask)
            counts["selectors"] += 1
            integral = all(v.denominator == 1 for v in values.values())
            if integral:
                counts["integral"] += 1
                if all(v > 0 for v in values.values()):
                    counts["positive"] += 1
        selector_counts[str(depth)] = dict(counts)
    assert [selector_counts[str(d)].get("integral", 0) for d in range(1, 5)] == [3, 5, 6, 1]
    assert [selector_counts[str(d)].get("positive", 0) for d in range(1, 5)] == [1, 0, 0, 0]

    cycles = set()
    bound = 100_000
    done = set()
    for start in range(1, bound + 1):
        path, position = [], {}
        n = start
        while n <= bound and n not in done and n not in position:
            position[n] = len(path)
            path.append(n)
            n = n // 2 if n % 2 == 0 else (3 * n + 1) // 2
        if n in position:
            cycle = tuple(path[position[n]:])
            cycles.add(min(cycle[i:] + cycle[:i] for i in range(len(cycle))))
        done.update(path)
    assert cycles == {(1, 2)}

    return {
        "parity_words_checked": word_count,
        "concatenation_seams_checked": seam_count,
        "adjacent_swaps_checked": swap_count,
        "sharp_response_pairs": sharp_response_pairs,
        "primitive_necklaces_through_16": necklaces,
        "integral_primitive_necklaces": integral_necklaces,
        "selector_counts": selector_counts,
        "hostile_examples": {
            "integral_vs_2_adic": "w=001 has fixed lift 4/5 and no integral lift",
            "completion_escape": "rho(1^k)=2^k-1 converges 2-adically to -1",
            "bounded_cycle_scan": "the induced graph on 1..100000 contains only the cycle 1,2",
        },
    }


def digit_length(n, b):
    length = 1
    while n >= b:
        n //= b
        length += 1
    return length


def profiles(b, m):
    if m == 1:
        for a in range(1, b):
            yield (2 * a,)
        return
    pairs = m // 2
    prefix_ranges = [range(1, 2 * b - 1)] + [range(2 * b - 1)] * (pairs - 1)
    if m % 2 == 0:
        yield from product(*prefix_ranges)
    else:
        for prefix in product(*prefix_ranges):
            for center in range(0, 2 * b - 1, 2):
                yield prefix + (center,)


def profile_value(b, m, q):
    pairs = m // 2
    total = sum(q[i] * (b**i + b**(m - 1 - i)) for i in range(pairs))
    if m % 2:
        total += q[pairs] * b**pairs
    return total


def reverse_digits(n, b):
    digits = []
    while n:
        digits.append(n % b)
        n //= b
    out = 0
    for digit in digits:
        out = out * b + digit
    return out


def reverse_add_audit():
    total_profiles = 0
    overlap_outputs = 0
    bases = {}
    for b in range(2, 13):
        fibres = defaultdict(list)
        per_length = {}
        for m in range(1, 9):
            chart = {}
            for q in profiles(b, m):
                y = profile_value(b, m, q)
                assert y not in chart
                chart[y] = q
                assert digit_length(y, b) in (m, m + 1)
                fibres[y].append((m, q))
            per_length[str(m)] = len(chart)
            total_profiles += len(chart)
        overlaps = 0
        for y, fibre in fibres.items():
            assert len(fibre) <= 2
            M = digit_length(y, b)
            assert all(m in (M - 1, M) for m, _ in fibre)
            assert len({M - m for m, _ in fibre}) == len(fibre)
            if len(fibre) == 2:
                overlaps += 1
        overlap_outputs += overlaps

        if b == 2:
            n0, n1 = 0b110, 0b1000
        else:
            n0, n1 = (b - 1) * b + 2, b * b + b
        y0, y1 = n0 + reverse_digits(n0, b), n1 + reverse_digits(n1, b)
        assert y0 == y1
        bases[str(b)] = {
            "profiles_by_length": per_length,
            "overlap_outputs": overlaps,
            "sharp_collision": {"sources": [n0, n1], "output": y0},
        }
    assert total_profiles == 1_208_130
    assert overlap_outputs == 310
    return {
        "profiles_checked": total_profiles,
        "overlap_outputs": overlap_outputs,
        "bases": bases,
        "hostile_boundary": "base 2 uses 110_2 and 1000_2, both mapping to 1001_2",
    }


def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("--output")
    args = parser.parse_args()
    report = {
        "schema": "fprd.result-review.collatz-reverse-add.v1",
        "arithmetic": "exact Python integers and fractions; no project imports",
        "collatz": collatz_audit(),
        "reverse_add": reverse_add_audit(),
        "status": "pass",
    }
    rendered = json.dumps(report, indent=2, sort_keys=True) + "\n"
    if args.output:
        with open(args.output, "w", encoding="utf-8") as handle:
            handle.write(rendered)
    print(rendered, end="")


if __name__ == "__main__":
    main()