# ================================================================
# QCA ORIGIN TEST v21
#
# PURPOSE:
#   Separate the two primitive collapses instead of treating
#   X+1 and X²-X-1 as one residual.
#
#   AXIS A : additive / translation collapse
#            X + 1 = 0
#
#   AXIS B : self-similar / fixed-point collapse
#            X² - X - 1 = 0
#
#   The state is represented exactly in Z[φ]:
#
#            X = a + bφ
#
#   with:
#
#            φ² = φ + 1
#
#   and exact norm:
#
#            N(a+bφ) = a² + ab - b²
#
#   NO FLOATS
#   NO STORED φ
#   NO EXTERNAL CONSTANTS
#   NO "INSIDE/OUTSIDE" CLASSIFICATION
# ================================================================

from fractions import Fraction
from math import gcd


# ================================================================
# EXACT Z[φ] ELEMENT
# ================================================================

def zadd(x, y):
    return (x[0] + y[0], x[1] + y[1])


def zsub(x, y):
    return (x[0] - y[0], x[1] - y[1])


def zneg(x):
    return (-x[0], -x[1])


def zscale(k, x):
    return (k * x[0], k * x[1])


def zmul(x, y):
    """
    (a+bφ)(c+dφ)

    = ac + (ad+bc)φ + bdφ²
    = (ac+bd) + (ad+bc+bd)φ
    """
    a, b = x
    c, d = y

    return (
        a * c + b * d,
        a * d + b * c + b * d
    )


def zsquare(x):
    return zmul(x, x)


def znorm(x):
    """
    Algebraic norm:

        N(a+bφ) = (a+bφ)(a+bψ)
                = a² + ab - b²

    where ψ = 1-φ.
    """
    a, b = x
    return a * a + a * b - b * b


def zzero(x):
    return x[0] == 0 and x[1] == 0


def zabsnorm(x):
    return abs(znorm(x))


def fmt(x):
    a, b = x

    if a == 0 and b == 0:
        return "(0+0φ)"

    if b == 0:
        return f"({a}+0φ)"

    if a == 0:
        return f"(0+{b}φ)"

    return f"({a}+{b}φ)"


# ================================================================
# PRIMITIVE OPERATORS
# ================================================================

ONE = (1, 0)
ZERO = (0, 0)
PHI = (0, 1)

TRINITY_MINUS = (-1, 0)
TRINITY_ZERO = (0, 0)
TRINITY_PLUS = (1, 0)


def primitive(x):
    """
    Δ0(X) = X + 1
    """
    return zadd(x, ONE)


def inverse_primitive(x):
    """
    Δ-0(X) = X - 1
    """
    return zsub(x, ONE)


def fixed_collapse(x):
    """
    Δ1(X) = X² - X - 1
    """
    return zsub(
        zsub(
            zsquare(x),
            x
        ),
        ONE
    )


# ================================================================
# ROOT / COLLAPSE TESTS
# ================================================================

def additive_residual(x):
    """
    Residual from:

        X + 1 = 0
    """
    return primitive(x)


def inverse_residual(x):
    """
    Residual from:

        X - 1 = 0
    """
    return inverse_primitive(x)


def fixed_residual(x):
    """
    Residual from:

        X² - X - 1 = 0
    """
    return fixed_collapse(x)


def additive_origin_exact(x):
    return additive_residual(x) == ZERO


def inverse_origin_exact(x):
    return inverse_residual(x) == ZERO


def fixed_origin_exact(x):
    return fixed_residual(x) == ZERO


# ================================================================
# FIBONACCI / φ ORBIT
#
# φ^0 = 1
# φ^1 = φ
# φ²  = 1+φ
# φ³  = 1+2φ
# φ⁴  = 2+3φ
# ...
#
# φ^n = F(n-1) + F(n)φ
# ================================================================

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

    x = ONE

    for _ in range(n):
        x = zmul(x, PHI)

    return x


# ================================================================
# EXACT FIBONACCI
# ================================================================

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

    a = 0
    b = 1

    for _ in range(n):
        a, b = b, a + b

    return a


def phi_power_fast(n):
    return (fibonacci(n - 1), fibonacci(n)) if n > 0 else ONE


# ================================================================
# PRIMALITY
#
# Used only as a label for the displayed test set.
# The QCA algebra itself does not depend on primality.
# ================================================================

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


# ================================================================
# EXACT RATIO
#
# No floating point.
# ================================================================

def ratio(num, den):
    if den == 0:
        return None

    return Fraction(num, den)


def ratio_text(r):
    if r is None:
        return "UNDEFINED"

    if r.denominator == 1:
        return str(r.numerator)

    return f"{r.numerator}/{r.denominator}"


# ================================================================
# OPERATOR PROFILE
# ================================================================

def operator_profile(x):
    nx = znorm(x)

    p = primitive(x)
    i = inverse_primitive(x)
    f = fixed_collapse(x)

    np = znorm(p)
    ni = znorm(i)
    nf = znorm(f)

    if nx == 0:
        rp = None
        ri = None
        rf = None
    else:
        rp = ratio(abs(np), abs(nx))
        ri = ratio(abs(ni), abs(nx))
        rf = ratio(abs(nf), abs(nx) * abs(nx))

    return {
        "x": x,
        "nx": nx,

        "primitive": p,
        "primitive_norm": np,
        "primitive_ratio": rp,

        "inverse": i,
        "inverse_norm": ni,
        "inverse_ratio": ri,

        "fixed": f,
        "fixed_norm": nf,
        "fixed_ratio": rf,
    }


# ================================================================
# TRANSLATION DUALITY
#
# Compare:
#
#     N(X+1)
#     N(X-1)
#
# Their difference has a simple exact form.
#
# N(X+1)-N(X-1)=2(a+b)
#
# for X=a+bφ.
# ================================================================

def translation_difference(x):
    return znorm(primitive(x)) - znorm(inverse_primitive(x))


def translation_sum(x):
    return znorm(primitive(x)) + znorm(inverse_primitive(x))


# ================================================================
# FIXED COLLAPSE EXPANSION
#
# For X=a+bφ:
#
# X²-X-1
#
# is computed entirely in Z[φ].
# ================================================================

def fixed_components(x):
    a, b = x

    # X²:
    #
    # a² + 2abφ + b²(φ+1)
    #
    # = (a²+b²) + (2ab+b²)φ

    x2_a = a * a + b * b
    x2_b = 2 * a * b + b * b

    # X²-X-1

    c = x2_a - a - 1
    d = x2_b - b

    return (c, d)


# ================================================================
# ORIGIN CLASSIFICATION
#
# This does NOT call anything "inside" or "outside".
#
# The state is described by which collapse, if any, is exact,
# and by the relative operator magnitudes.
# ================================================================

def classify_origin(x):
    p = primitive(x)
    i = inverse_primitive(x)
    f = fixed_collapse(x)

    exact = []

    if zzero(p):
        exact.append("ADDITIVE(+1)")

    if zzero(i):
        exact.append("ADDITIVE(-1)")

    if zzero(f):
        exact.append("FIXED(φ)")

    if exact:
        return "EXACT:" + ",".join(exact)

    return "NONZERO-RESIDUAL"


# ================================================================
# DOMINANT OPERATOR
#
# Dominance is based on normalized exact ratios.
#
# Primitive / inverse:
#
#     |N(X±1)| / |N(X)|
#
# Fixed:
#
#     |N(X²-X-1)| / |N(X)|²
#
# ================================================================

def dominant_operator(profile):
    x = profile["x"]
    nx = profile["nx"]

    if nx == 0:
        return "N(X)=0"

    candidates = [
        ("primitive", profile["primitive_ratio"]),
        ("inverse", profile["inverse_ratio"]),
        ("fixed", profile["fixed_ratio"]),
    ]

    candidates = [
        item for item in candidates
        if item[1] is not None
    ]

    # Compare exact Fractions.
    candidates.sort(key=lambda item: item[1])

    return candidates[0][0]


# ================================================================
# PRIME ORBIT / PERTURBATION
#
# The primary orbit is φ^n.
#
# The perturbation is not another "answer".
# It is the additive translation of the orbit index:
#
#     φ^n -> φ^n - φ^(n-1)
#
# which is deliberately kept visible rather than collapsed into
# a Boolean classification.
# ================================================================

def orbit_perturbation(n):
    if n <= 0:
        return ZERO

    return zsub(
        phi_power_fast(n),
        phi_power_fast(n - 1)
    )


# ================================================================
# PRINT PROFILE
# ================================================================

def print_operator(name, value, norm_value, r):
    print(f"{name} = {fmt(value)}")
    print(f"N    = {norm_value}")

    if r is None:
        print("R    = UNDEFINED")
    else:
        print(f"R    = {ratio_text(r)}")

    print()


def print_case(n):
    prime_label = "prime" if is_prime(n) else "composite"

    x = phi_power_fast(n)
    p = operator_profile(x)

    print()
    print("=" * 64)
    print(f"{n:3d} {prime_label}")
    print("=" * 64)

    print()
    print("INPUT:")
    print(f"Π = {fmt(x)}")
    print(f"N = {p['nx']}")

    print()
    print("ORIGIN CLASS:")
    print(classify_origin(x))

    print()
    print("AXIS A — ADDITIVE COLLAPSE")
    print()
    print("Equation:")
    print("X+1=0")
    print()
    print("Residual:")
    print(f"Δ0 = {fmt(p['primitive'])}")
    print(f"N  = {p['primitive_norm']}")
    print(f"R  = {ratio_text(p['primitive_ratio'])}")

    print()
    print("DUAL ADDITIVE COLLAPSE")
    print()
    print("Equation:")
    print("X-1=0")
    print()
    print("Residual:")
    print(f"Δ-0 = {fmt(p['inverse'])}")
    print(f"N   = {p['inverse_norm']}")
    print(f"R   = {ratio_text(p['inverse_ratio'])}")

    print()
    print("TRANSLATION DUALITY")
    print()
    print(f"N(X+1)-N(X-1) = {translation_difference(x)}")
    print(f"N(X+1)+N(X-1) = {translation_sum(x)}")

    print()
    print("AXIS B — FIXED-POINT COLLAPSE")
    print()
    print("Equation:")
    print("X²-X-1=0")
    print()
    print("Residual:")
    print(f"Δ1 = {fmt(p['fixed'])}")
    print(f"N  = {p['fixed_norm']}")
    print(f"R  = {ratio_text(p['fixed_ratio'])}")

    print()
    print("TRINITY NEIGHBORHOOD")
    print()
    print(f"X-1 = {fmt(inverse_primitive(x))}")
    print(f"X   = {fmt(x)}")
    print(f"X+1 = {fmt(primitive(x))}")

    print()
    print("PERTURBATION")
    print()

    q = orbit_perturbation(n)

    print(f"Πp = {fmt(q)}")
    print(f"N  = {znorm(q)}")

    print()
    print("DOMINANT NORMALIZED OPERATOR:")
    print(dominant_operator(p))

    print()
    print("INTERPRETATION:")
    print(
        "translation and fixed-point collapse remain separate axes"
    )

    print()


# ================================================================
# SYMBOLIC ORIGIN DISPLAY
# ================================================================

def print_header():
    print()
    print("=" * 64)
    print("QCA ORIGIN TEST v21")
    print("=" * 64)

    print()
    print("SUBSTRATE:")
    print()
    print("X = 0")

    print()
    print()
    print("CO-EMERGENCE:")
    print()
    print("(-1,0,+1)")

    print()
    print()
    print("AXIS A:")
    print()
    print("X+1=0")
    print()
    print("ROOT:")
    print()
    print("X=-1")

    print()
    print()
    print("AXIS B:")
    print()
    print("X²-X-1=0")

    print()
    print("FIXED POINT:")
    print()
    print("Ω = α = φ")

    print()
    print()
    print("IDENTITY:")
    print()
    print("φ²=φ+1")

    print()
    print()
    print("TRINITY:")
    print()
    print("alpha-1=(-1+1φ)")
    print("alpha  =( 0+1φ)")
    print("alpha+1=( 1+1φ)")

    print()
    print()
    print("ALGEBRA:")
    print()
    print("X = a+bφ")
    print()
    print("N(a+bφ)=a²+ab-b²")

    print()
    print()
    print("OPERATORS:")
    print()
    print("Δ0(X)   = X+1")
    print("Δ-0(X)  = X-1")
    print("Δ1(X)   = X²-X-1")

    print()
    print()
    print("NORMALIZATION:")
    print()
    print("|N(X+1)| / |N(X)|")
    print("|N(X-1)| / |N(X)|")
    print("|N(X²-X-1)| / |N(X)|²")

    print()
    print()
    print("=" * 64)


# ================================================================
# CONSISTENCY TESTS
# ================================================================

def run_consistency_tests():
    print()
    print("=" * 64)
    print("EXACT CONSISTENCY TESTS")
    print("=" * 64)

    # φ² = φ+1
    lhs = zsquare(PHI)
    rhs = zadd(PHI, ONE)

    print()
    print("1. φ² = φ+1")
    print(f"   φ²     = {fmt(lhs)}")
    print(f"   φ+1    = {fmt(rhs)}")
    print(f"   PASS   = {lhs == rhs}")

    # Norm multiplicativity
    a = (3, 5)
    b = (2, 7)

    product = zmul(a, b)

    lhs_norm = znorm(product)
    rhs_norm = znorm(a) * znorm(b)

    print()
    print("2. N(XY)=N(X)N(Y)")
    print(f"   N(XY)    = {lhs_norm}")
    print(f"   N(X)N(Y) = {rhs_norm}")
    print(f"   PASS     = {lhs_norm == rhs_norm}")

    # φ norm
    print()
    print("3. N(φ)")
    print(f"   N(φ) = {znorm(PHI)}")
    print("   PASS =", znorm(PHI) == -1)

    # Fibonacci orbit
    print()
    print("4. Fibonacci embedding")

    passed = True

    for n in range(1, 16):
        x1 = phi_power(n)
        x2 = phi_power_fast(n)

        ok = x1 == x2
        passed = passed and ok

        print(
            f"   n={n:2d} "
            f"{fmt(x1):>24s} "
            f"PASS={ok}"
        )

    print()
    print("   OVERALL =", passed)

    # Norm alternation
    print()
    print("5. Norm orbit")

    passed = True

    for n in range(1, 16):
        x = phi_power_fast(n)
        actual = znorm(x)
        expected = -1 if n % 2 else 1

        ok = actual == expected
        passed = passed and ok

        print(
            f"   n={n:2d} "
            f"N(φ^n)={actual:3d} "
            f"expected={expected:3d} "
            f"PASS={ok}"
        )

    print()
    print("   OVERALL =", passed)


# ================================================================
# SUMMARY
# ================================================================

def print_summary(values):
    print()
    print("=" * 64)
    print("ORIGIN SUMMARY")
    print("=" * 64)

    print()
    print(
        f"{'n':>3} "
        f"{'N(X)':>12} "
        f"{'R+':>12} "
        f"{'R-':>12} "
        f"{'Rφ':>18} "
        f"{'DOMINANT':>12}"
    )

    print("-" * 64)

    for n in values:
        x = phi_power_fast(n)
        p = operator_profile(x)

        print(
            f"{n:3d} "
            f"{p['nx']:12d} "
            f"{ratio_text(p['primitive_ratio']):>12} "
            f"{ratio_text(p['inverse_ratio']):>12} "
            f"{ratio_text(p['fixed_ratio']):>18} "
            f"{dominant_operator(p):>12}"
        )

    print()
    print("=" * 64)


# ================================================================
# MAIN
# ================================================================

def main():

    print_header()

    run_consistency_tests()

    values = [
        2,
        3,
        4,
        5,
        6,
        7,
        8,
        9,
        11,
        13,
        17,
        19,
        23,
        29,
        31,
    ]

    for n in values:
        print_case(n)

    print_summary(values)

    print()
    print("=" * 64)
    print("END")
    print("=" * 64)
    print()


if __name__ == "__main__":
    main()