commutative_lrs_review.py

#!/usr/bin/env python3
"""Independent finite audit for the CPA and indexed-LRS review units."""

from __future__ import annotations

import json
import math
from fractions import Fraction
from itertools import product
from pathlib import Path


OUT = Path(__file__).with_name("audit-output.json")


def lcm(values: list[int]) -> int:
    result = 1
    for value in values:
        result = math.lcm(result, value)
    return result


def multiply_polynomials(left: list[int], right: list[int]) -> list[int]:
    result = [0] * (len(left) + len(right) - 1)
    for i, a in enumerate(left):
        for j, b in enumerate(right):
            result[i + j] += a * b
    return result


def factorial_audit() -> dict[str, int]:
    transition_checks = 0
    derivative_checks = 0
    x, y = 1, 1
    for n in range(25):
        assert x == math.factorial(n)
        assert y == n + 1
        transition_checks += 2
        x, y = x * y, y + 1

    # The same length-only computation is a valid unary control system and a
    # two-letter control system with identical transitions on both letters.
    for alphabet_size in (1, 2):
        for n in range(13):
            for r in range(13):
                rising = math.prod(range(n + 1, n + r + 1))
                assert math.factorial(n + r) == math.factorial(n) * rising
                derivative_checks += alphabet_size

    # P_r(y)=y(y+1)...(y+r-1) is monic of degree r.  Therefore the first
    # R+1 rising factorials have unit triangular coefficient matrix and span
    # every polynomial of degree at most R.
    triangular_checks = 0
    polynomial = [1]
    for r in range(13):
        if r:
            polynomial = multiply_polynomials(polynomial, [r - 1, 1])
        assert len(polynomial) == r + 1
        assert polynomial[-1] == 1
        triangular_checks += len(polynomial)

    cusp_exponents = {
        2 * a + 3 * b for a in range(20) for b in range(20)
    }
    assert 0 in cusp_exponents
    assert 1 not in cusp_exponents
    for exponent in range(2, 40):
        assert exponent in cusp_exponents

    return {
        "factorial_transition_equalities": transition_checks,
        "unary_and_binary_derivative_equalities": derivative_checks,
        "rising_factorial_coefficient_checks": triangular_checks,
        "cusp_exponents_checked": 40,
    }


def recurrence_terms(coefficients: tuple[Fraction, ...], initial: tuple[Fraction, ...], limit: int) -> list[Fraction]:
    d = len(coefficients)
    terms = list(initial)
    while len(terms) <= limit:
        k = len(terms) - d
        terms.append(sum(coefficients[i] * terms[k + i] for i in range(d)))
    return terms


def matmul(left: list[list[int]], right: list[list[int]]) -> list[list[int]]:
    rows, inner, cols = len(left), len(right), len(right[0])
    assert len(left[0]) == inner
    return [
        [sum(left[i][k] * right[k][j] for k in range(inner)) for j in range(cols)]
        for i in range(rows)
    ]


def matpow(matrix: list[list[int]], exponent: int) -> list[list[int]]:
    size = len(matrix)
    result = [[int(i == j) for j in range(size)] for i in range(size)]
    base = matrix
    while exponent:
        if exponent & 1:
            result = matmul(result, base)
        base = matmul(base, base)
        exponent //= 2
    return result


def is_prime(value: int) -> bool:
    if value < 2:
        return False
    for divisor in range(2, math.isqrt(value) + 1):
        if value % divisor == 0:
            return False
    return True


def lrs_audit() -> dict[str, int]:
    values = (Fraction(-1), Fraction(-1, 2), Fraction(0), Fraction(1, 2), Fraction(1))
    normalization_cases = 0
    term_checks = 0
    block_checks = 0
    prime_certificates = 0

    for d in (1, 2, 3):
        # Exhaust d<=2 and use a deterministic 1-in-13 slice at d=3 so the
        # certificate remains fast while covering genuinely third-order data.
        for index, data in enumerate(product(values, repeat=2 * d)):
            if d == 3 and index % 13:
                continue
            coefficients = tuple(data[:d])
            initial = tuple(data[d:])
            q = lcm([value.denominator for value in coefficients])
            e = lcm([value.denominator for value in initial])
            terms = recurrence_terms(coefficients, initial, 10)
            scaled = [Fraction(e * (q**k)) * terms[k] for k in range(11)]
            assert all(value.denominator == 1 for value in scaled)
            integer_terms = [int(value) for value in scaled]
            integer_coefficients = [int(coefficients[i] * (q ** (d - i))) for i in range(d)]
            for k in range(11 - d):
                assert integer_terms[k + d] == sum(
                    integer_coefficients[i] * integer_terms[k + i]
                    for i in range(d)
                )
                assert (terms[k] == 0) == (integer_terms[k] == 0)
                term_checks += 2

            # Column companion representation g_n=e_1^T A^n beta and the
            # displayed block matrix identity (M^(n+2))_(1,d+2)=g_n.
            companion = [[0 for _ in range(d)] for _ in range(d)]
            for i in range(d - 1):
                companion[i][i + 1] = 1
            companion[-1] = integer_coefficients
            beta = [[integer_terms[i]] for i in range(d)]
            for n in range(8):
                value = matmul([[1] + [0] * (d - 1)], matmul(matpow(companion, n), beta))[0][0]
                assert value == integer_terms[n]

                size = d + 2
                block = [[0 for _ in range(size)] for _ in range(size)]
                for j in range(d):
                    block[0][1 + j] = int(j == 0)
                    block[1 + j][size - 1] = beta[j][0]
                    for k in range(d):
                        block[1 + j][1 + k] = companion[j][k]
                assert matpow(block, n + 2)[0][size - 1] == integer_terms[n]
                block_checks += 2

            nonzero = next((abs(value) for value in integer_terms if value), None)
            if nonzero is not None:
                witness = next(p for p in range(2, 1000) if is_prime(p) and nonzero % p)
                assert nonzero % witness
                prime_certificates += 1
            normalization_cases += 1

    # Hostile denominator example: multiplying the whole sequence by one
    # common denominator does not make the recurrence coefficient integral,
    # whereas the index-dependent Q^k scaling does.
    hostile = recurrence_terms((Fraction(1, 2),), (Fraction(1),), 12)
    assert hostile[4] == Fraction(1, 16)
    assert [int((2**k) * hostile[k]) for k in range(13)] == [1] * 13

    return {
        "normalized_recurrences": normalization_cases,
        "normalization_and_zero_equalities": term_checks,
        "companion_and_block_matrix_equalities": block_checks,
        "escaping_prime_certificates": prime_certificates,
        "hostile_denominator_cases": 1,
    }


def main() -> None:
    output = {
        "protocol": "fprd-result-review.v1",
        "review_units": [
            "commutative-polynomial-automata-series",
            "indexed-zeroness-linear-recurrence-sequences",
        ],
        "commutative_polynomial_automata": factorial_audit(),
        "indexed_lrs": lrs_audit(),
        "status": "pass",
    }
    OUT.write_text(json.dumps(output, indent=2, sort_keys=True) + "\n", encoding="utf-8")
    print(json.dumps(output, indent=2, sort_keys=True))


if __name__ == "__main__":
    main()