#!/usr/bin/env python3
"""
======================================================================
HDGL / Z[phi] v25
COMPLETE FROBENIUS MATRIX CLOSURE / COLLISION HUNTER
======================================================================

Primitive closure:

    phi^2 = phi + 1

Orbit:

    phi^n = F_(n-1) + F_n phi

Norm:

    N(phi^n) = (-1)^n

Lucas:

    L_n = 2 F_(n-1) + F_n

Axis A:

    L_n == 1 (mod n)

Axis B:

    F_n == Frobenius branch (mod n)

Axis C:

    COMPLETE FIBONACCI MATRIX CLOSURE (mod n)

The recurrence matrix is:

             [ 1  1 ]
    A   =    [ 1  0 ]

and exactly:

             [ F_(n+1)  F_n     ]
    A^n =    [ F_n      F_(n-1) ]

For prime p != 5:

    s = (5/p)

where s is generated by the discriminant of:

    x^2 - x - 1 = 0

For s = +1:

    A^p == A (mod p)

For s = -1:

    A^p == I - A (mod p)

Equivalently:

             [ (1+s)/2    s       ]
    A^p ==  [ s           (1-s)/2 ]  (mod p)

For p = 5:

    A^5 == 3I (mod 5)

This script:

    1. Computes Fibonacci pairs by binary fast doubling.
    2. Computes Lucas numbers exactly.
    3. Tests Axis A.
    4. Tests Axis B.
    5. Tests complete matrix closure Axis C.
    6. Scans all odd n <= LIMIT.
    7. Separates primes from composites using a sieve.
    8. Reports every collision at each stage.
    9. Tests the previously observed:
           705
           2465
           2737
           3745
           4181
           5777
           6721
   10. Verifies all internal Fibonacci/Lucas/matrix identities.

No floating point.
No stored phi.
No decimal conversion of giant integers.
No external constants except those structurally generated by
the quadratic closure and its discriminant.
======================================================================
"""

# ======================================================================
# CONFIGURATION
# ======================================================================

LIMIT = 10_000


# ======================================================================
# FAST DOUBLING
# ======================================================================

def fibonacci_pair(n):
    """
    Exact iterative binary fast doubling.

    Returns:

        (F_n, F_(n+1))
    """

    if n < 0:
        raise ValueError("n must be non-negative")

    a = 0
    b = 1

    for bit in bin(n)[2:]:
        # F_(2k)
        c = a * ((b << 1) - a)

        # F_(2k+1)
        d = a * a + b * b

        if bit == "0":
            a, b = c, d
        else:
            a, b = d, c + d

    return a, b


# ======================================================================
# FIBONACCI
# ======================================================================

def fibonacci(n):
    return fibonacci_pair(n)[0]


# ======================================================================
# LUCAS
# ======================================================================

def lucas_fast(n):
    """
    L_n = F_(n-1) + F_(n+1)

        = 2F_(n-1) + F_n

        = 2F_(n+1) - F_n
    """

    if n == 0:
        return 2

    f_nm1, f_n = fibonacci_pair(n - 1)

    return (f_nm1 << 1) + f_n


def lucas_from_pair(n):
    """
    L_n from (F_n, F_(n+1)).
    """

    f_n, f_np1 = fibonacci_pair(n)

    return (f_np1 << 1) - f_n


# ======================================================================
# FIBONACCI MATRIX
# ======================================================================

def fibonacci_matrix(n):
    """
    Return the exact matrix:

             [ F_(n+1)  F_n     ]
    A^n  =   [ F_n      F_(n-1) ]

    represented as:

        (m00, m01, m10, m11)
    """

    if n == 0:
        return (1, 0, 0, 1)

    f_n, f_np1 = fibonacci_pair(n)

    f_nm1 = f_np1 - f_n

    return (
        f_np1,
        f_n,
        f_n,
        f_nm1
    )


# ======================================================================
# MATRIX MODULO n
# ======================================================================

def fibonacci_matrix_mod(n):
    """
    A^n modulo n.
    """

    m00, m01, m10, m11 = fibonacci_matrix(n)

    return (
        m00 % n,
        m01 % n,
        m10 % n,
        m11 % n
    )


# ======================================================================
# MATRIX MULTIPLICATION
# ======================================================================

def matrix_mul(A, B, mod=None):
    """
    2x2 matrix multiplication.

    Matrices represented as:

        (a,b,c,d)

    corresponding to:

        [a b]
        [c d]
    """

    a, b, c, d = A
    e, f, g, h = B

    r00 = a * e + b * g
    r01 = a * f + b * h
    r10 = c * e + d * g
    r11 = c * f + d * h

    if mod is not None:
        return (
            r00 % mod,
            r01 % mod,
            r10 % mod,
            r11 % mod
        )

    return (
        r00,
        r01,
        r10,
        r11
    )


# ======================================================================
# MATRIX POWER
# ======================================================================

def matrix_pow(A, n, mod=None):
    """
    Exact binary matrix exponentiation.

    Included as an independent verification of the
    Fibonacci-coordinate construction.
    """

    result = (1, 0, 0, 1)
    base = A

    while n:
        if n & 1:
            result = matrix_mul(result, base, mod)

        n >>= 1

        if n:
            base = matrix_mul(base, base, mod)

    return result


# ======================================================================
# PRIMITIVE MATRIX
# ======================================================================

A_MATRIX = (
    1, 1,
    1, 0
)

I_MATRIX = (
    1, 0,
    0, 1
)


# ======================================================================
# DISCRIMINANT BRANCH
# ======================================================================

def frobenius_sign(n):
    """
    Branch generated by the discriminant:

        Delta = 1 + 4 = 5

    For primes p != 5:

        p mod 5 = 1,4  -> +1
        p mod 5 = 2,3  -> -1

    For n divisible by 5:

        None

    This function is used only to construct the expected
    prime Frobenius branch.
    """

    r = n % 5

    if r in (1, 4):
        return 1

    if r in (2, 3):
        return -1

    return None


# ======================================================================
# AXIS A
# ======================================================================

def axis_a_residue(n):
    """
    Axis A:

        R_A(n) = (L_n - 1) mod n
    """

    return (lucas_fast(n) - 1) % n


def axis_a_passes(n):
    return axis_a_residue(n) == 0


# ======================================================================
# AXIS B
# ======================================================================

def expected_fibonacci_residue(n):
    """
    Prime Frobenius expectation:

        F_p == +1 mod p
        or
        F_p == -1 mod p

    depending on the discriminant-5 branch.

    p = 5 is handled separately.
    """

    s = frobenius_sign(n)

    if s is None:
        return None

    if s == 1:
        return 1

    return n - 1


def fibonacci_residue(n):
    f_n, _ = fibonacci_pair(n)

    return f_n % n


def axis_b_passes(n):
    """
    Axis B:

        F_n == expected Frobenius residue (mod n)

    p = 5 is accepted as the exceptional discriminant case.
    """

    if n == 5:
        return True

    expected = expected_fibonacci_residue(n)

    if expected is None:
        return False

    return fibonacci_residue(n) == expected


# ======================================================================
# AXIS C
# COMPLETE MATRIX FROBENIUS CLOSURE
# ======================================================================

def expected_frobenius_matrix(n):
    """
    Construct the complete matrix expected for a prime.

    For s = +1:

             [1  1]
        A == [1  0]

    For s = -1:

             [ 0 -1]
        I-A =[ -1 1]

    modulo n.

    For n = 5:

        A^5 == 3I (mod 5)
    """

    if n == 5:
        return (
            3 % n, 0,
            0, 3 % n
        )

    s = frobenius_sign(n)

    if s is None:
        return None

    if s == 1:
        return (
            1 % n,
            1 % n,
            1 % n,
            0
        )

    return (
        0,
        (-1) % n,
        (-1) % n,
        1 % n
    )


def axis_c_matrix(n):
    """
    Actual A^n modulo n.
    """

    return fibonacci_matrix_mod(n)


def axis_c_passes(n):
    """
    Complete matrix closure test.
    """

    expected = expected_frobenius_matrix(n)

    if expected is None:
        return False

    actual = axis_c_matrix(n)

    return actual == expected


# ======================================================================
# COMPLETE FILTER
# ======================================================================

def combined_pass(n):
    """
    All three axes.
    """

    return (
        axis_a_passes(n)
        and axis_b_passes(n)
        and axis_c_passes(n)
    )


# ======================================================================
# PRIME SIEVE
# ======================================================================

def prime_sieve(limit):
    """
    Exact Eratosthenes sieve.

    Used ONLY as ground truth for validation.

    It is not part of the HDGL/Frobenius filter.
    """

    prime = bytearray(b"\x01") * (limit + 1)

    if limit >= 0:
        prime[0] = 0

    if limit >= 1:
        prime[1] = 0

    p = 2

    while p * p <= limit:

        if prime[p]:

            start = p * p

            count = ((limit - start) // p) + 1

            prime[start:limit + 1:p] = b"\x00" * count

        p += 1

    return prime


# ======================================================================
# MATRIX FORMAT
# ======================================================================

def format_matrix(M):
    a, b, c, d = M

    return (
        f"[[{a}, {b}], "
        f"[{c}, {d}]]"
    )


# ======================================================================
# BASIC PRIME VALIDATION
# ======================================================================

TRAIN_PRIMES = [
    3, 5, 7, 11, 13,
    17, 19, 23, 29, 31
]

ODD_COMPOSITES = [
    9, 15, 21, 25, 27,
    33, 35, 39, 45, 49,
    51, 55, 57, 63, 65,
    69, 75, 77, 81, 85,
    87, 91, 93, 95, 99
]


print("=" * 72)
print("AXIS A: (L_n - 1) % n")
print("=" * 72)

prime_a = []

for n in TRAIN_PRIMES:

    r = axis_a_residue(n)

    prime_a.append((n, r))

print("PRIMES:")
print(prime_a)

print()
print("ODD COMPOSITES:")

composite_a = []

for n in ODD_COMPOSITES:

    r = axis_a_residue(n)

    composite_a.append((n, r))

print(composite_a)


# ======================================================================
# AXIS B
# ======================================================================

print()
print("=" * 72)
print("AXIS B: F_n FROBENIUS RESIDUE")
print("=" * 72)

print()
print("n        F_n mod n       expected")

for n in TRAIN_PRIMES:

    actual = fibonacci_residue(n)

    expected = expected_fibonacci_residue(n)

    print(
        f"{n:3d}      "
        f"{actual:8d}      "
        f"{str(expected):>8s}"
    )


# ======================================================================
# AXIS C
# ======================================================================

print()
print("=" * 72)
print("AXIS C: COMPLETE FROBENIUS MATRIX")
print("=" * 72)

for n in TRAIN_PRIMES:

    actual = axis_c_matrix(n)

    expected = expected_frobenius_matrix(n)

    print()
    print(f"n = {n}")
    print("actual:   ", format_matrix(actual))
    print("expected: ", format_matrix(expected))
    print("PASS:", actual == expected)


# ======================================================================
# COMBINED TRAINING BLOCK
# ======================================================================

print()
print("=" * 72)
print("COMBINED AXIS A + B + C")
print("=" * 72)

print()
print(
    "n        prime   A-residue   B-pass   C-pass   combined"
)

for n in TRAIN_PRIMES:

    a = axis_a_residue(n)
    b = axis_b_passes(n)
    c = axis_c_passes(n)
    combined = a == 0 and b and c

    print(
        f"{n:3d}      "
        f"{str(True):5s}   "
        f"{a:9d}   "
        f"{str(b):6s}   "
        f"{str(c):6s}   "
        f"{str(combined):8s}"
    )

print()

for n in ODD_COMPOSITES:

    a = axis_a_residue(n)
    b = axis_b_passes(n)
    c = axis_c_passes(n)
    combined = a == 0 and b and c

    print(
        f"{n:3d}      "
        f"{str(False):5s}   "
        f"{a:9d}   "
        f"{str(b):6s}   "
        f"{str(c):6s}   "
        f"{str(combined):8s}"
    )


# ======================================================================
# PREVIOUS COLLISION SET
# ======================================================================

KNOWN_COLLISIONS = [
    705,
    2465,
    2737,
    3745,
    4181,
    5777,
    6721
]

print()
print("=" * 72)
print("PREVIOUS AXIS-A COLLISIONS")
print("=" * 72)

print()
print(
    "n        A-residue   F-residue   "
    "B-pass   actual matrix              C-pass"
)

for n in KNOWN_COLLISIONS:

    a = axis_a_residue(n)

    f = fibonacci_residue(n)

    b = axis_b_passes(n)

    actual = axis_c_matrix(n)

    c = axis_c_passes(n)

    print(
        f"{n:5d}    "
        f"{a:9d}   "
        f"{f:9d}   "
        f"{str(b):6s}   "
        f"{format_matrix(actual):27s} "
        f"{str(c):6s}"
    )


# ======================================================================
# FULL SCAN
# ======================================================================

print()
print("=" * 72)
print(f"SCANNING ODD n: 3 <= n <= {LIMIT}")
print("=" * 72)

prime_table = prime_sieve(LIMIT)

axis_a_collisions = []
axis_b_collisions = []
axis_c_collisions = []
combined_collisions = []

prime_a_failures = []
prime_b_failures = []
prime_c_failures = []
prime_combined_failures = []


for n in range(3, LIMIT + 1, 2):

    prime = bool(prime_table[n])

    a = axis_a_passes(n)
    b = axis_b_passes(n)
    c = axis_c_passes(n)

    combined = a and b and c

    if prime:

        if not a:
            prime_a_failures.append(n)

        if not b:
            prime_b_failures.append(n)

        if not c:
            prime_c_failures.append(n)

        if not combined:
            prime_combined_failures.append(n)

    else:

        if a:
            axis_a_collisions.append(n)

        if b:
            axis_b_collisions.append(n)

        if c:
            axis_c_collisions.append(n)

        if combined:
            combined_collisions.append(n)


# ======================================================================
# RESULTS
# ======================================================================

print()
print("=" * 72)
print("RESULTS")
print("=" * 72)

print()
print("Prime failures on Axis A:")
print(prime_a_failures)

print()
print("Prime failures on Axis B:")
print(prime_b_failures)

print()
print("Prime failures on Axis C:")
print(prime_c_failures)

print()
print("Prime failures on A + B + C:")
print(prime_combined_failures)

print()
print("Composite survivors of Axis A:")
print(axis_a_collisions)

print()
print("Composite survivors of Axis B:")
print(axis_b_collisions)

print()
print("Composite survivors of Axis C:")
print(axis_c_collisions)

print()
print("Composite survivors of ALL THREE:")
print(combined_collisions)

print()
print("Axis A composite survivors:",
      len(axis_a_collisions))

print("Axis B composite survivors:",
      len(axis_b_collisions))

print("Axis C composite survivors:",
      len(axis_c_collisions))

print("Combined composite survivors:",
      len(combined_collisions))


# ======================================================================
# COMPLETE COLLISION REPORT
# ======================================================================

print()
print("=" * 72)
print("COMPLETE A+B+C COLLISION REPORT")
print("=" * 72)

if not combined_collisions:

    print()
    print("NO COMPOSITE SURVIVORS FOUND.")

else:

    print()

    for n in combined_collisions:

        A_actual = axis_c_matrix(n)
        A_expected = expected_frobenius_matrix(n)

        print(
            f"n={n:6d}  "
            f"Lres={axis_a_residue(n):6d}  "
            f"Fres={fibonacci_residue(n):6d}  "
            f"matrix={format_matrix(A_actual)}"
        )

        print(
            f"         expected={format_matrix(A_expected)}"
        )


# ======================================================================
# FIRST AXIS-A COLLISIONS
# ======================================================================

print()
print("=" * 72)
print("FIRST AXIS-A COLLISIONS WITH COMPLETE MATRIX DATA")
print("=" * 72)

for n in axis_a_collisions:

    print()
    print(f"n = {n}")
    print(f"L_n - 1 mod n = {axis_a_residue(n)}")
    print(f"F_n mod n      = {fibonacci_residue(n)}")
    print(
        f"Axis B         = "
        f"{'PASS' if axis_b_passes(n) else 'FAIL'}"
    )

    print(
        "Actual A^n     =",
        format_matrix(axis_c_matrix(n))
    )

    expected = expected_frobenius_matrix(n)

    if expected is not None:

        print(
            "Expected prime =",
            format_matrix(expected)
        )

    print(
        "Axis C         =",
        "PASS" if axis_c_passes(n) else "FAIL"
    )


# ======================================================================
# MATRIX IDENTITY VALIDATION
# ======================================================================

print()
print("=" * 72)
print("MATRIX IDENTITY VALIDATION")
print("=" * 72)

matrix_failures = []

for n in range(0, 500):

    direct = matrix_pow(A_MATRIX, n)

    fibmat = fibonacci_matrix(n)

    if direct != fibmat:

        matrix_failures.append(
            (n, direct, fibmat)
        )


if matrix_failures:

    print("FAILURES:")

    for failure in matrix_failures[:20]:
        print(failure)

else:

    print(
        "ALL MATRIX POWERS MATCH "
        "THE FIBONACCI COORDINATE FORM."
    )


# ======================================================================
# LUCAS IDENTITY VALIDATION
# ======================================================================

print()
print("=" * 72)
print("LUCAS / FIBONACCI IDENTITY VALIDATION")
print("=" * 72)

identity_failures = []

for n in range(1, 1000):

    f_nm1, f_n = fibonacci_pair(n - 1)

    L1 = (f_nm1 << 1) + f_n

    L2 = lucas_fast(n)

    L3 = lucas_from_pair(n)

    if not (L1 == L2 == L3):

        identity_failures.append(
            ("Lucas", n, L1, L2, L3)
        )

    f_n, _ = fibonacci_pair(n)

    L_n = lucas_fast(n)

    # L_n^2 - 5 F_n^2 = 4(-1)^n

    lhs = L_n * L_n - 5 * f_n * f_n

    rhs = 4 if (n & 1) == 0 else -4

    if lhs != rhs:

        identity_failures.append(
            ("Lucas-Fibonacci", n, lhs, rhs)
        )


if identity_failures:

    print("FAILURES:")

    for failure in identity_failures[:20]:
        print(failure)

else:

    print("ALL EXACT IDENTITIES PASSED.")


# ======================================================================
# DETERMINANT IDENTITY
# ======================================================================

print()
print("=" * 72)
print("MATRIX DETERMINANT IDENTITY")
print("=" * 72)

det_failures = []

for n in range(0, 1000):

    m00, m01, m10, m11 = fibonacci_matrix(n)

    determinant = m00 * m11 - m01 * m10

    expected = -1 if (n & 1) else 1

    if determinant != expected:

        det_failures.append(
            (n, determinant, expected)
        )


if det_failures:

    print("FAILURES:")

    for failure in det_failures[:20]:
        print(failure)

else:

    print(
        "ALL DETERMINANTS SATISFY:"
    )

    print(
        "det(A^n) = (-1)^n"
    )


# ======================================================================
# MATRIX TRACE / LUCAS IDENTITY
# ======================================================================

print()
print("=" * 72)
print("MATRIX TRACE = LUCAS VALIDATION")
print("=" * 72)

trace_failures = []

for n in range(0, 1000):

    m00, m01, m10, m11 = fibonacci_matrix(n)

    trace = m00 + m11

    L = lucas_fast(n)

    if trace != L:

        trace_failures.append(
            (n, trace, L)
        )


if trace_failures:

    print("FAILURES:")

    for failure in trace_failures[:20]:
        print(failure)

else:

    print(
        "ALL TRACES SATISFY:"
    )

    print(
        "trace(A^n) = L_n"
    )


# ======================================================================
# PRIME FROBENIUS MATRIX VALIDATION
# ======================================================================

print()
print("=" * 72)
print("PRIME FROBENIUS MATRIX VALIDATION")
print("=" * 72)

frobenius_failures = []

for p in range(2, LIMIT + 1):

    if not prime_table[p]:
        continue

    actual = fibonacci_matrix_mod(p)

    expected = expected_frobenius_matrix(p)

    if actual != expected:

        frobenius_failures.append(
            (p, actual, expected)
        )


if frobenius_failures:

    print("FAILURES:")

    for failure in frobenius_failures[:20]:
        print(failure)

else:

    print(
        f"ALL PRIMES <= {LIMIT} "
        "PASS COMPLETE FROBENIUS MATRIX CLOSURE."
    )


# ======================================================================
# COLLISION STRUCTURE
# ======================================================================

print()
print("=" * 72)
print("COLLISION STRUCTURE")
print("=" * 72)

print()
print(
    "A = Lucas scalar closure"
)

print(
    "B = Fibonacci discriminant branch"
)

print(
    "C = Complete 2x2 orbit closure"
)

print()

print(
    "Axis A survivors:",
    axis_a_collisions
)

print()

print(
    "Axis A + B survivors:",
    [
        n for n in axis_a_collisions
        if axis_b_passes(n)
    ]
)

print()

print(
    "Axis A + B + C survivors:",
    combined_collisions
)


# ======================================================================
# FINAL STRUCTURAL SUMMARY
# ======================================================================

print()
print("=" * 72)
print("HDGL CLOSURE")
print("=" * 72)

print(
"""
Primitive closure:

    phi^2 = phi + 1

Orbit:

    phi^n = F_(n-1) + F_n phi

Norm:

    N(phi^n) = (-1)^n

Translation duality:

    N(phi^n + 1) - N(phi^n - 1) = 2L_n

Lucas:

    L_n = 2F_(n-1) + F_n

Axis A:

    L_n == 1 (mod n)

Axis B:

    F_n == Frobenius branch (mod n)

Complete orbit:

             [ F_(n+1)  F_n     ]
    A^n  =   [ F_n      F_(n-1) ]

Prime Frobenius closure:

    A^p == A       (mod p)    when  p == +/-1 (mod 5)

    A^p == I - A   (mod p)    when  p == +/-2 (mod 5)

Exceptional discriminant branch:

    A^5 == 3I      (mod 5)

Axis C therefore tests the complete orbit state,
rather than only one or two scalar projections.

The collision set is measured explicitly rather than assumed away.
"""
)

print()
print("=" * 72)
print("DONE")
print("=" * 72)