#!/usr/bin/env python3
"""
HDGL / Z[phi] v24
Lucas-Frobenius Congruence Collision Hunter

Primitive closure:
    phi^2 = phi + 1

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

Lucas:
    L_n = 2 F_(n-1) + F_n

Axis A:
    R_A(n) = (L_n - 1) mod n

Frobenius companion:
    For prime p != 5,

        F_p == +1 (mod p)  if p == +/-1 (mod 5)
        F_p == -1 (mod p)  if p == +/-2 (mod 5)

The 5 arises from the discriminant:
    Delta = 1 + 4 = 5

No floating point.
No stored phi.
No decimal conversion of giant integers.
"""

# ------------------------------------------------------------
# FAST DOUBLING
# ------------------------------------------------------------

def fibonacci_pair(n):
    """
    Return (F_n, F_(n+1)) exactly.

    Uses binary fast doubling:
        F_(2k)   = F_k (2F_(k+1) - F_k)
        F_(2k+1) = F_k^2 + F_(k+1)^2
    """
    if n == 0:
        return 0, 1

    a, b = fibonacci_pair(n >> 1)

    c = a * ((b << 1) - a)
    d = a * a + b * b

    if n & 1:
        return d, c + d

    return c, d


# ------------------------------------------------------------
# ITERATIVE FAST DOUBLING
# Avoids recursion depth entirely.
# ------------------------------------------------------------

def fibonacci_pair_iter(n):
    """
    Exact iterative fast doubling.
    Returns (F_n, F_(n+1)).
    """
    a = 0
    b = 1

    for bit in bin(n)[2:]:
        c = a * ((b << 1) - a)
        d = a * a + b * b

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

    return a, b


# ------------------------------------------------------------
# LUCAS
# ------------------------------------------------------------

def lucas_fast(n):
    """
    L_n = F_(n-1) + F_(n+1)
        = 2F_(n-1) + F_n
    """
    if n == 0:
        return 2

    f_n_minus_1, f_n = fibonacci_pair_iter(n - 1)

    return (f_n_minus_1 << 1) + f_n


def lucas_from_pair(n):
    """
    Compute L_n from the pair (F_n, F_(n+1)).
    """
    f_n, f_np1 = fibonacci_pair_iter(n)

    # L_n = F_(n-1) + F_(n+1)
    # F_(n-1) = F_(n+1) - F_n
    #
    # Therefore:
    # L_n = 2F_(n+1) - F_n
    return (f_np1 << 1) - f_n


# ------------------------------------------------------------
# AXIS A
# ------------------------------------------------------------

def axis_a_residue(n):
    """
    R_A(n) = (L_n - 1) mod n
    """
    if n <= 0:
        raise ValueError("n must be positive")

    L = lucas_fast(n)

    return (L - 1) % n


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


# ------------------------------------------------------------
# DISCRIMINANT / FROBENIUS AXIS
# ------------------------------------------------------------

def expected_fibonacci_residue(n):
    """
    Prime Frobenius expectation for F_n mod n.

    The quadratic closure is:

        x^2 - x - 1 = 0

    whose discriminant is:

        Delta = 1 + 4 = 5

    For prime p != 5:

        F_p == +1 mod p    if p == +/-1 mod 5
        F_p == -1 mod p    if p == +/-2 mod 5
    """
    r = n % 5

    if r in (1, 4):
        return 1

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

    # n == 0 mod 5
    return None


def fibonacci_residue(n):
    f_n, _ = fibonacci_pair_iter(n)
    return f_n % n


def axis_b_passes(n):
    """
    Test the Frobenius/Fibonacci congruence.

    n == 5 is handled separately.
    """
    if n == 5:
        return True

    expected = expected_fibonacci_residue(n)

    if expected is None:
        return False

    return fibonacci_residue(n) == expected


# ------------------------------------------------------------
# COMBINED FILTER
# ------------------------------------------------------------

def combined_pass(n):
    """
    Axis A + Axis B.
    """
    return axis_a_passes(n) and axis_b_passes(n)


# ------------------------------------------------------------
# SIMPLE EXACT PRIMALITY
# Used ONLY as ground truth.
# It is NOT part of the HDGL filters.
# ------------------------------------------------------------

def is_prime(n):
    if n < 2:
        return False

    if n == 2:
        return True

    if n % 2 == 0:
        return False

    d = 3

    while d * d <= n:
        if n % d == 0:
            return False
        d += 2

    return True


# ------------------------------------------------------------
# BASIC 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)


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}      {actual:8d}      {expected}")


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

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

for n in TRAIN_PRIMES:
    p = is_prime(n)
    a = axis_a_residue(n)
    b = axis_b_passes(n)
    c = combined_pass(n)

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

print()

for n in ODD_COMPOSITES:
    p = is_prime(n)
    a = axis_a_residue(n)
    b = axis_b_passes(n)
    c = combined_pass(n)

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


# ------------------------------------------------------------
# COLLISION HUNTER
# ------------------------------------------------------------

LIMIT = 10_000

axis_a_collisions = []
combined_collisions = []

prime_a_failures = []
prime_combined_failures = []


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

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

    prime = is_prime(n)

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

    if prime:
        if not a:
            prime_a_failures.append(n)

        if not c:
            prime_combined_failures.append(n)

    else:
        if a:
            axis_a_collisions.append(n)

        if c:
            combined_collisions.append(n)


# ------------------------------------------------------------
# REPORT
# ------------------------------------------------------------

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

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

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

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

print()
print("Composite survivors of BOTH axes:")
print(combined_collisions)

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

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


# ------------------------------------------------------------
# FIRST COLLISIONS
# ------------------------------------------------------------

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

for n in axis_a_collisions[:50]:
    print(
        f"n={n:5d}  "
        f"Lres={axis_a_residue(n):5d}  "
        f"Fres={fibonacci_residue(n):5d}  "
        f"B={'PASS' if axis_b_passes(n) else 'FAIL'}"
    )


print()
print("=" * 72)
print("FIRST COMBINED COLLISIONS")
print("=" * 72)

for n in combined_collisions[:50]:
    print(
        f"n={n:5d}  "
        f"Lres={axis_a_residue(n):5d}  "
        f"Fres={fibonacci_residue(n):5d}"
    )


# ------------------------------------------------------------
# INTERNAL CONSISTENCY
# ------------------------------------------------------------

print()
print("=" * 72)
print("INTERNAL IDENTITY CHECKS")
print("=" * 72)

identity_failures = []

for n in range(1, 1000):

    f_nm1, f_n = fibonacci_pair_iter(n - 1)

    L1 = (f_nm1 << 1) + f_n
    L2 = lucas_from_pair(n)

    if L1 != L2:
        identity_failures.append(
            ("Lucas", n, L1, L2)
        )

    # Fundamental Lucas/Fibonacci identity:
    #
    # L_n^2 - 5 F_n^2 = 4(-1)^n
    #
    f_n, _ = fibonacci_pair_iter(n)
    L_n = lucas_fast(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.")


# ------------------------------------------------------------
# FINAL INTERPRETATION
# ------------------------------------------------------------

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

print("""
phi^2 = phi + 1

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

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

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

L_n = 2 F_(n-1) + F_n

Axis A:
    (L_n - 1) mod n == 0

Axis B:
    F_n mod n follows the discriminant-5 Frobenius branch

The scan does NOT assume that either congruence proves primality.

It searches for the composite collision set.
""")