#!/usr/bin/env python3
# =============================================================================
# HDGL / Z[phi] — VERSION 23
# =============================================================================
#
# EXACT ORBIT
# TRANSLATION SYMMETRY
# LUCAS COORDINATE
# LARGE-PRIME TEST
# CANDIDATE INVARIANT SEARCH
#
# VERSION 23 FIXES:
#
#   * No float arithmetic.
#   * No decimal conversion of giant integers.
#   * No Python 4300-digit string-conversion failure.
#   * Exact Z[phi] coefficient-pair arithmetic.
#   * Exact fast-doubling Fibonacci.
#   * Exact Lucas coordinate.
#   * Exact norm.
#   * Exact translation identities.
#   * Large n = 100003 works without converting 20,000+ digit
#     integers to decimal strings.
#   * Candidate invariant search.
#   * Held-out validation.
#   * Explicit counterexample showing that the zero-sum identity
#     is an oddness invariant rather than, by itself, a primality test.
#
# CORE:
#
#     phi^2 = phi + 1
#
#     X_n = phi^n
#         = F_(n-1) + F_n phi
#
#     N(a+b phi) = a^2 + ab - b^2
#
#     N(phi^n) = (-1)^n
#
# TRANSLATIONS:
#
#     Delta_0(X)  = X + 1
#     Delta_-0(X) = X - 1
#
#     N(X+1) + N(X-1)
#         = 2(N(X)+1)
#         = 2(1+(-1)^n)
#
#     N(X+1) - N(X-1)
#         = 2 L_n
#
# THEREFORE:
#
#     N(X+1) = 1 + (-1)^n + L_n
#     N(X-1) = 1 + (-1)^n - L_n
#
# IMPORTANT:
#
#     The zero translation sum detects parity.
#     It does NOT by itself prove primality.
#
# VERSION 23 searches for an additional structural invariant rather
# than silently equating oddness with primality.
# =============================================================================

import time
from math import gcd


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

BIG_N = 100003

# Representative small values.
SMALL_CASES = [
    0, 1, 2, 3, 4, 5, 6, 7, 8, 9,
    10, 11, 12, 13, 15, 17, 19, 21, 23,
    25, 27, 29, 31
]

# Training primes.
TRAIN_PRIMES = [
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
    31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
    73, 79, 83, 89, 97
]

# Training composites.
TRAIN_COMPOSITES = [
    4, 6, 8, 9, 10, 12, 14, 15, 16, 18,
    20, 21, 22, 24, 25, 26, 27, 28, 30, 32,
    33, 34, 35, 36, 38, 39, 40, 42, 44, 45,
    46, 48, 49, 50, 51, 52, 54, 55, 56, 57,
    58, 60, 62, 63, 64, 65, 66, 68, 69, 70,
    72, 74, 75, 76, 77, 78, 80, 81, 82, 84,
    85, 86, 87, 88, 90, 91, 92, 93, 94, 95,
    96, 98, 99, 100
]

# Held-out primes.
VALIDATE_PRIMES = [
    101, 103, 107, 109, 113, 127, 131, 137,
    139, 149, 151, 157, 163, 167, 173, 179,
    181, 191, 193, 197, 199
]

# Held-out composites.
VALIDATE_COMPOSITES = [
    102, 104, 105, 106, 108, 110, 111, 112,
    114, 115, 116, 117, 118, 119, 120, 121,
    122, 123, 124, 125, 126, 128, 129, 130,
    132, 133, 134, 135, 136, 138, 140, 141,
    142, 143, 144, 145, 146, 147, 148, 150
]


# =============================================================================
# EXACT Z[phi]
#
# Every element is represented as:
#
#     (a,b) == a + b phi
#
# with:
#
#     phi^2 = phi + 1
#
# No floating-point phi exists anywhere in the calculation.
# =============================================================================

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


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 phi)(c+d phi)

    phi^2 = phi + 1

    Therefore:

        ac + bd + (ad+bc+bd)phi
    """

    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 zpow(x, n):
    """
    Exact integer-exponentiation in Z[phi].
    """

    result = ONE
    base = x

    while n > 0:
        if n & 1:
            result = zmul(result, base)

        base = zmul(base, base)
        n >>= 1

    return result


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

        N(a+b phi) = a^2 + ab - b^2
    """

    a, b = x

    return (
        a * a
        + a * b
        - b * b
    )


def zfmt(x):
    """
    Human-readable exact representation.
    """

    a, b = x

    if a == 0 and b == 0:
        return "0"

    if b == 0:
        return str(a)

    if a == 0:
        if b == 1:
            return "phi"

        if b == -1:
            return "-phi"

        return f"{b}phi"

    if b == 1:
        return f"{a}+phi"

    if b == -1:
        return f"{a}-phi"

    if b > 0:
        return f"{a}+{b}phi"

    return f"{a}{b}phi"


# =============================================================================
# SAFE INTEGER MAGNITUDE
#
# NEVER converts a giant integer to a decimal string.
#
# Python 3.11+ intentionally limits:
#
#     str(huge_integer)
#
# to 4300 digits by default.
#
# This routine computes the decimal digit count without that conversion.
#
# No floating point is used.
# =============================================================================

def integer_digits(n):
    """
    Exact decimal digit count without str(n).

    The initial estimate uses the rational lower-resolution approximation
    to log10(2):

        log10(2) ~= 30103 / 100000

    The final result is corrected using exact powers of ten.
    """

    n = abs(n)

    if n == 0:
        return 1

    # Initial estimate.
    k = (n.bit_length() * 30103) // 100000

    if k < 1:
        k = 1

    # Exact correction.
    p = 10 ** (k - 1)

    while n < p:
        k -= 1
        p //= 10

    while n >= p * 10:
        p *= 10
        k += 1

    return k


def integer_bit_length(n):
    return abs(n).bit_length()


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

def delta_0(x):
    """
    Primitive positive translation:

        X -> X + 1
    """

    return zadd(x, ONE)


def delta_minus_0(x):
    """
    Primitive negative translation:

        X -> X - 1
    """

    return zsub(x, ONE)


def fixed_closure(x):
    """
    Quadratic closure:

        X^2 - X - 1
    """

    return zsub(
        zsub(
            zsquare(x),
            x
        ),
        ONE
    )


# =============================================================================
# FIBONACCI FAST DOUBLING
#
# Returns:
#
#     (F_n, F_(n+1))
#
# exactly.
# =============================================================================

def fibonacci_pair(n):
    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)


def fibonacci(n):
    if n < 0:
        raise ValueError("negative Fibonacci index")

    return fibonacci_pair(n)[0]


# =============================================================================
# LUCAS
#
#     L_n = F_(n-1) + F_(n+1)
#
# with:
#
#     L_0 = 2
# =============================================================================

def lucas(n):
    if n == 0:
        return 2

    fn, fn1 = fibonacci_pair(n)

    fn_minus_1 = fn1 - fn

    return fn_minus_1 + fn1


# =============================================================================
# PHI ORBIT
#
#     phi^n = F_(n-1) + F_n phi
#
# n=0:
#
#     phi^0 = 1
# =============================================================================

def phi_orbit(n):
    if n < 0:
        raise ValueError("negative orbit index")

    if n == 0:
        return ONE

    fn, fn1 = fibonacci_pair(n)

    fn_minus_1 = fn1 - fn

    return (
        fn_minus_1,
        fn
    )


# =============================================================================
# ORBIT DATA
# =============================================================================

def orbit_data(n):
    x = phi_orbit(n)

    xp = delta_0(x)
    xm = delta_minus_0(x)

    fc = fixed_closure(x)

    nx = znorm(x)
    nxp = znorm(xp)
    nxm = znorm(xm)
    nfc = znorm(fc)

    s = nxp + nxm
    d = nxp - nxm

    L = lucas(n)

    expected_n = -1 if (n & 1) else 1
    expected_sum = 0 if (n & 1) else 4
    expected_difference = 2 * L

    return {
        "n": n,

        "x": x,
        "xp": xp,
        "xm": xm,
        "fixed": fc,

        "N": nx,
        "N_plus": nxp,
        "N_minus": nxm,
        "N_fixed": nfc,

        "sum": s,
        "difference": d,
        "lucas": L,

        "expected_N": expected_n,
        "expected_sum": expected_sum,
        "expected_difference": expected_difference,

        "norm_identity": nx == expected_n,
        "sum_identity": s == expected_sum,
        "difference_identity": d == expected_difference,
    }


# =============================================================================
# EXACT PRIMALITY
#
# Used ONLY as ground truth while searching for a new invariant.
#
# This is not claimed as the HDGL mechanism.
# =============================================================================

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

    if n == 2 or n == 3:
        return True

    if n % 2 == 0:
        return False

    if n % 3 == 0:
        return False

    d = 5
    step = 2

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

        d += step
        step = 6 - step

    return True


# =============================================================================
# EXACT IDENTITY CHECKS
# =============================================================================

def check_norm_identity(n):
    d = orbit_data(n)
    return d["norm_identity"]


def check_sum_identity(n):
    d = orbit_data(n)
    return d["sum_identity"]


def check_difference_identity(n):
    d = orbit_data(n)
    return d["difference_identity"]


def check_all_identities(n):
    d = orbit_data(n)

    return (
        d["norm_identity"]
        and d["sum_identity"]
        and d["difference_identity"]
    )


# =============================================================================
# NORM MULTIPLICATIVITY
# =============================================================================

def check_norm_multiplication(x, y):
    lhs = znorm(zmul(x, y))
    rhs = znorm(x) * znorm(y)

    return lhs == rhs


# =============================================================================
# STRUCTURAL INVARIANTS
#
# These are all generated from the primitive translations and closure.
# =============================================================================

def invariants(n):
    d = orbit_data(n)

    x = d["x"]
    xp = d["xp"]
    xm = d["xm"]
    fc = d["fixed"]

    nx = d["N"]
    np = d["N_plus"]
    nm = d["N_minus"]
    nf = d["N_fixed"]

    s = np + nm
    diff = np - nm

    L = d["lucas"]

    abs_np = abs(np)
    abs_nm = abs(nm)

    translation_gcd = gcd(
        abs_np,
        abs_nm
    )

    return {
        # Direct norms.
        "N(X)": nx,
        "N(X+1)": np,
        "N(X-1)": nm,

        # Translation geometry.
        "SUM": s,
        "DIFF": diff,
        "ABS_SUM": abs(s),
        "ABS_DIFF": abs(diff),

        # GCD structures.
        "GCD_TRANSLATIONS": translation_gcd,
        "GCD_N_X_PLUS": gcd(abs(nx), abs_np),
        "GCD_N_X_MINUS": gcd(abs(nx), abs_nm),

        # Closure.
        "N(FIXED)": nf,
        "ABS_N(FIXED)": abs(nf),

        # Norm powers.
        "N(X)^2": nx * nx,

        # Exact identity residuals.
        "R_NORM": nx - d["expected_N"],
        "R_SUM": s - d["expected_sum"],
        "R_DIFF": diff - 2 * L,

        # Modular structural quantities.
        "SUM_MOD_N": s % n if n else 0,
        "DIFF_MOD_N": diff % n if n else 0,
        "ABS_DIFF_MOD_N": abs(diff) % n if n else 0,

        "LUCAS": L,
        "LUCAS_MOD_N": L % n if n else 0,

        "N_DIVIDES_DIFF": (
            n != 0 and diff % n == 0
        ),

        "N_DIVIDES_ABS_DIFF": (
            n != 0 and abs(diff) % n == 0
        ),

        "N_DIVIDES_LUCAS": (
            n != 0 and L % n == 0
        ),

        "N_DIVIDES_SUM": (
            n != 0 and s % n == 0
        ),

        # GCD against Lucas.
        "GCD_N_LUCAS": gcd(abs(nx), abs(L)),
        "GCD_TRANSLATIONS_LUCAS": gcd(
            translation_gcd,
            abs(L)
        ),
    }


# =============================================================================
# CANDIDATE SEARCH
# =============================================================================

def invariant_names():
    return list(invariants(7).keys())


def candidate_zero_separates(name):
    prime_states = set()
    composite_states = set()

    for n in TRAIN_PRIMES:
        prime_states.add(
            invariants(n)[name] == 0
        )

    for n in TRAIN_COMPOSITES:
        composite_states.add(
            invariants(n)[name] == 0
        )

    if len(prime_states) != 1:
        return False

    if len(composite_states) != 1:
        return False

    return (
        next(iter(prime_states))
        != next(iter(composite_states))
    )


def candidate_divisibility_separates(name):
    prime_states = set()
    composite_states = set()

    for n in TRAIN_PRIMES:
        v = invariants(n)[name]
        prime_states.add(v % n == 0)

    for n in TRAIN_COMPOSITES:
        v = invariants(n)[name]
        composite_states.add(v % n == 0)

    if len(prime_states) != 1:
        return False

    if len(composite_states) != 1:
        return False

    return (
        next(iter(prime_states))
        != next(iter(composite_states))
    )


def search_candidates():
    names = invariant_names()

    zero_candidates = []
    divisibility_candidates = []

    for name in names:
        if candidate_zero_separates(name):
            zero_candidates.append(name)

        if candidate_divisibility_separates(name):
            divisibility_candidates.append(name)

    print()
    print("=" * 78)
    print("STRUCTURAL CANDIDATE SEARCH")
    print("=" * 78)

    print()
    print("ZERO-SEPARATING CANDIDATES:")

    if zero_candidates:
        for name in zero_candidates:
            print(f"  {name}")
    else:
        print("  NONE")

    print()
    print("DIVISIBILITY-SEPARATING CANDIDATES:")

    if divisibility_candidates:
        for name in divisibility_candidates:
            print(f"  n | {name}")
    else:
        print("  NONE")

    return (
        zero_candidates,
        divisibility_candidates
    )


# =============================================================================
# HELD-OUT VALIDATION
# =============================================================================

def validate_zero_candidate(name):
    prime_ok = all(
        invariants(n)[name] == 0
        for n in VALIDATE_PRIMES
    )

    composite_ok = all(
        invariants(n)[name] != 0
        for n in VALIDATE_COMPOSITES
    )

    return prime_ok and composite_ok


def validate_divisibility_candidate(name):
    prime_ok = all(
        invariants(n)[name] % n == 0
        for n in VALIDATE_PRIMES
    )

    composite_ok = all(
        invariants(n)[name] % n != 0
        for n in VALIDATE_COMPOSITES
    )

    return prime_ok and composite_ok


def validate_candidates(
    zero_candidates,
    divisibility_candidates
):
    print()
    print("=" * 78)
    print("HELD-OUT VALIDATION")
    print("=" * 78)

    print()
    print("ZERO CANDIDATES:")

    if not zero_candidates:
        print("  NONE")

    for name in zero_candidates:
        result = validate_zero_candidate(name)

        print(
            f"  {name:<32}"
            f"{'VALIDATED' if result else 'REJECTED'}"
        )

    print()
    print("DIVISIBILITY CANDIDATES:")

    if not divisibility_candidates:
        print("  NONE")

    for name in divisibility_candidates:
        result = validate_divisibility_candidate(name)

        print(
            f"  n | {name:<28}"
            f"{'VALIDATED' if result else 'REJECTED'}"
        )


# =============================================================================
# EXACT ALGEBRA CONSISTENCY
# =============================================================================

def consistency_tests():
    print()
    print("=" * 78)
    print("EXACT ALGEBRA CONSISTENCY")
    print("=" * 78)

    passed = 0
    total = 0

    def test(name, condition):
        nonlocal passed
        nonlocal total

        total += 1

        if condition:
            passed += 1
            print(f"[PASS] {name}")
        else:
            print(f"[FAIL] {name}")

    # phi^2 = phi + 1
    test(
        "phi^2 = phi + 1",
        zmul(PHI, PHI)
        == zadd(PHI, ONE)
    )

    # N(phi) = -1
    test(
        "N(phi) = -1",
        znorm(PHI) == -1
    )

    # Norm orbit.
    test(
        "N(phi^n)=(-1)^n for n=0..64",
        all(
            check_norm_identity(n)
            for n in range(65)
        )
    )

    # Translation sum.
    test(
        "translation sum identity n=0..64",
        all(
            check_sum_identity(n)
            for n in range(65)
        )
    )

    # Translation difference.
    test(
        "translation difference identity n=0..64",
        all(
            check_difference_identity(n)
            for n in range(65)
        )
    )

    # Norm multiplicativity.
    multiplication_pairs = [
        (0, 0),
        (1, 1),
        (1, 2),
        (2, 3),
        (3, 5),
        (5, 8),
        (7, 11),
        (13, 17),
        (19, 23),
        (31, 37),
    ]

    test(
        "norm multiplicativity",
        all(
            check_norm_multiplication(
                phi_orbit(a),
                phi_orbit(b)
            )
            for a, b in multiplication_pairs
        )
    )

    # Fibonacci embedding.
    embedding_ok = True

    for n in range(1, 65):
        fn, fn1 = fibonacci_pair(n)

        expected = (
            fn1 - fn,
            fn
        )

        if phi_orbit(n) != expected:
            embedding_ok = False
            break

    test(
        "Fibonacci embedding",
        embedding_ok
    )

    # Direct powers versus Fibonacci orbit.
    powers_ok = True

    for n in range(0, 33):
        if zpow(PHI, n) != phi_orbit(n):
            powers_ok = False
            break

    test(
        "direct phi powers = Fibonacci orbit",
        powers_ok
    )

    print()
    print(
        f"RESULT: {passed}/{total} exact tests passed"
    )


# =============================================================================
# SMALL STRUCTURAL TABLE
# =============================================================================

def structural_table():
    values = [
        2, 3, 4, 5, 6, 7, 8, 9,
        10, 11, 12, 13, 15, 17, 19,
        21, 23, 25, 27, 29, 31
    ]

    print()
    print("=" * 78)
    print("STRUCTURAL ORBIT TABLE")
    print("=" * 78)

    print(
        f"{'n':>5} "
        f"{'P/C':>5} "
        f"{'N(X)':>6} "
        f"{'SUM':>8} "
        f"{'DIFF':>14} "
        f"{'L_n':>14} "
        f"{'n|DIFF':>8}"
    )

    print("-" * 78)

    for n in values:
        d = orbit_data(n)

        print(
            f"{n:>5} "
            f"{'P' if is_prime(n) else 'C':>5} "
            f"{d['N']:>6} "
            f"{d['sum']:>8} "
            f"{d['difference']:>14} "
            f"{d['lucas']:>14} "
            f"{str(d['difference'] % n == 0):>8}"
        )


# =============================================================================
# ONE CASE
# =============================================================================

def print_case(n):
    d = orbit_data(n)

    print()
    print("-" * 78)
    print(
        f"n = {n}    "
        f"{'PRIME' if is_prime(n) else 'COMPOSITE'}"
    )
    print("-" * 78)

    print(f"X             = {zfmt(d['x'])}")
    print(f"X + 1         = {zfmt(d['xp'])}")
    print(f"X - 1         = {zfmt(d['xm'])}")

    print()
    print(f"N(X)          = {d['N']}")
    print(f"N(X + 1)      = {d['N_plus']}")
    print(f"N(X - 1)      = {d['N_minus']}")

    print()
    print(f"SUM            = {d['sum']}")
    print(f"EXPECTED SUM   = {d['expected_sum']}")

    print()
    print(f"DIFFERENCE     = {d['difference']}")
    print(f"2*L_n          = {d['expected_difference']}")
    print(f"L_n            = {d['lucas']}")

    print()
    print(f"FIXED CLOSURE  = {zfmt(d['fixed'])}")
    print(f"N(FIXED)       = {d['N_fixed']}")

    print()
    print(f"NORM IDENTITY  = {d['norm_identity']}")
    print(f"SUM IDENTITY   = {d['sum_identity']}")
    print(f"DIFF IDENTITY  = {d['difference_identity']}")
    print(
        f"ALL IDENTITIES = "
        f"{d['norm_identity'] and d['sum_identity'] and d['difference_identity']}"
    )


# =============================================================================
# LARGE EXACT TEST
# =============================================================================

def large_test(n):
    print()
    print("=" * 78)
    print(f"LARGE EXACT ORBIT TEST: n = {n}")
    print("=" * 78)

    start = time.perf_counter()

    d = orbit_data(n)

    elapsed = time.perf_counter() - start

    a, b = d["x"]

    print()
    print(f"n                  = {n}")
    print(
        f"ground-truth       = "
        f"{'PRIME' if is_prime(n) else 'COMPOSITE'}"
    )

    print()
    print("EXACT COORDINATE SIZE")
    print(
        f"digits(F_(n-1))    = "
        f"{integer_digits(a)}"
    )
    print(
        f"digits(F_n)        = "
        f"{integer_digits(b)}"
    )

    print(
        f"bits(F_(n-1))      = "
        f"{integer_bit_length(a)}"
    )
    print(
        f"bits(F_n)          = "
        f"{integer_bit_length(b)}"
    )

    print()
    print("EXACT NORM")

    print(
        f"N(X)               = "
        f"{d['N']}"
    )

    print()
    print("TRANSLATION CLOSURE")

    print(
        f"N(X+1)+N(X-1)      = "
        f"{d['sum']}"
    )

    print(
        f"expected sum       = "
        f"{d['expected_sum']}"
    )

    print(
        f"sum identity       = "
        f"{d['sum'] == d['expected_sum']}"
    )

    print()
    print("LUCAS SEPARATION")

    print(
        f"Lucas digits        = "
        f"{integer_digits(d['lucas'])}"
    )

    print(
        f"difference matches  = "
        f"{d['difference'] == 2 * d['lucas']}"
    )

    print(
        f"difference digits   = "
        f"{integer_digits(d['difference'])}"
    )

    print()
    print("TIMING")

    print(
        f"exact orbit time    = "
        f"{elapsed:.6f} sec"
    )

    print()
    print("STRUCTURAL RESULT")
    print()
    print("    N(phi^n)              = (-1)^n")
    print("    N(phi^n+1)+N(phi^n-1) = 2(1+(-1)^n)")
    print("    N(phi^n+1)-N(phi^n-1) = 2L_n")

    if n & 1:
        print()
        print("    odd orbit:")
        print("        translation sum = 0")
    else:
        print()
        print("    even orbit:")
        print("        translation sum = 4")

    print()
    print("NO GIANT INTEGER WAS CONVERTED TO DECIMAL TEXT.")


# =============================================================================
# LARGE NUMBER REPORT
# =============================================================================

def report_large_number(n):
    d = orbit_data(n)

    a, b = d["x"]

    print()
    print("=" * 78)
    print("LARGE-N REPRESENTATION REPORT")
    print("=" * 78)

    print()
    print(f"n               : {n}")
    print(
        f"prime           : "
        f"{is_prime(n)}"
    )

    print()
    print("COORDINATE A = F_(n-1)")
    print(
        f"    decimal digits : "
        f"{integer_digits(a)}"
    )
    print(
        f"    binary bits    : "
        f"{integer_bit_length(a)}"
    )

    print()
    print("COORDINATE B = F_n")
    print(
        f"    decimal digits : "
        f"{integer_digits(b)}"
    )
    print(
        f"    binary bits    : "
        f"{integer_bit_length(b)}"
    )

    print()
    print("LUCAS COORDINATE")
    print(
        f"    decimal digits : "
        f"{integer_digits(d['lucas'])}"
    )
    print(
        f"    binary bits    : "
        f"{integer_bit_length(d['lucas'])}"
    )

    print()
    print("NORM")
    print(
        f"    N(X) = {d['N']}"
    )

    print()
    print("TRANSLATION SUM")
    print(
        f"    N(X+1)+N(X-1) = {d['sum']}"
    )

    print()
    print("TRANSLATION DIFFERENCE")
    print(
        f"    digits = "
        f"{integer_digits(d['difference'])}"
    )

    print()
    print("THE ACTUAL COORDINATE PAIR IS NOT PRINTED.")
    print("THIS AVOIDS GIANT DECIMAL STRING CONVERSION.")


# =============================================================================
# PARITY COUNTEREXAMPLE
#
# 100005 is composite but odd.
#
# Therefore the zero translation sum cannot itself be a prime detector.
# =============================================================================

def demonstrate_parity_not_primality():
    n = 100005

    d = orbit_data(n)

    print()
    print("=" * 78)
    print("PARITY != PRIMALITY")
    print("=" * 78)

    print()
    print(f"n                  = {n}")
    print(
        f"ground-truth       = "
        f"{'PRIME' if is_prime(n) else 'COMPOSITE'}"
    )

    print(
        f"N(phi^n)           = "
        f"{d['N']}"
    )

    print(
        f"N(X+1)+N(X-1)      = "
        f"{d['sum']}"
    )

    print()
    print("RESULT:")
    print()
    print("    zero translation sum")
    print("            =>")
    print("    odd orbit")
    print()
    print("but NOT necessarily:")
    print()
    print("    zero translation sum")
    print("            =>")
    print("    prime exponent")


# =============================================================================
# PRIME / COMPOSITE COMPARISON
# =============================================================================

def compare_prime_composite():
    print()
    print("=" * 78)
    print("PRIME / COMPOSITE STRUCTURAL COMPARISON")
    print("=" * 78)

    print()
    print("PRIMES")
    print()
    print(
        f"{'n':>5} "
        f"{'N(X)':>8} "
        f"{'SUM':>8} "
        f"{'DIFF%n':>12} "
        f"{'L_n%n':>12}"
    )

    print("-" * 55)

    for n in TRAIN_PRIMES[:20]:
        d = orbit_data(n)

        print(
            f"{n:>5} "
            f"{d['N']:>8} "
            f"{d['sum']:>8} "
            f"{d['difference'] % n:>12} "
            f"{d['lucas'] % n:>12}"
        )

    print()
    print("COMPOSITES")
    print()
    print(
        f"{'n':>5} "
        f"{'N(X)':>8} "
        f"{'SUM':>8} "
        f"{'DIFF%n':>12} "
        f"{'L_n%n':>12}"
    )

    print("-" * 55)

    for n in TRAIN_COMPOSITES[:20]:
        d = orbit_data(n)

        print(
            f"{n:>5} "
            f"{d['N']:>8} "
            f"{d['sum']:>8} "
            f"{d['difference'] % n:>12} "
            f"{d['lucas'] % n:>12}"
        )


# =============================================================================
# DERIVED EQUATIONS
# =============================================================================

def print_derived_equations():
    print()
    print("=" * 78)
    print("DERIVED EQUATIONS")
    print("=" * 78)

    print()
    print("PRIMITIVE")
    print("    X = 0")

    print()
    print("TRANSLATIONS")
    print("    Delta_0(X)  = X + 1")
    print("    Delta_-0(X) = X - 1")

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

    print()
    print("ORBIT")
    print("    X_n = phi^n")
    print("        = F_(n-1) + F_n phi")

    print()
    print("NORM")
    print("    N(a+b phi) = a^2 + ab - b^2")

    print()
    print("ORBIT NORM")
    print("    N(X_n) = (-1)^n")

    print()
    print("TRANSLATION SUM")
    print("    N(X_n+1) + N(X_n-1)")
    print("        = 2(N(X_n)+1)")
    print("        = 2(1+(-1)^n)")

    print()
    print("TRANSLATION DIFFERENCE")
    print("    N(X_n+1) - N(X_n-1)")
    print("        = 2L_n")

    print()
    print("SEPARATED STATES")
    print("    N(X_n+1) = 1 + (-1)^n + L_n")
    print("    N(X_n-1) = 1 + (-1)^n - L_n")


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

def main():
    print()
    print("=" * 78)
    print("HDGL / Z[phi] — VERSION 23")
    print("EXACT LARGE-INTEGER ORBIT ARCHITECTURE")
    print("=" * 78)

    # -------------------------------------------------------------------------
    # Algebra.
    # -------------------------------------------------------------------------

    print_derived_equations()

    # -------------------------------------------------------------------------
    # Internal consistency.
    # -------------------------------------------------------------------------

    consistency_tests()

    # -------------------------------------------------------------------------
    # Small structural table.
    # -------------------------------------------------------------------------

    structural_table()

    # -------------------------------------------------------------------------
    # Prime/composite observation.
    # -------------------------------------------------------------------------

    compare_prime_composite()

    # -------------------------------------------------------------------------
    # Demonstrate that zero-sum = oddness, not automatically primality.
    # -------------------------------------------------------------------------

    demonstrate_parity_not_primality()

    # -------------------------------------------------------------------------
    # Search for additional structural candidates.
    # -------------------------------------------------------------------------

    zero_candidates, divisibility_candidates = search_candidates()

    # -------------------------------------------------------------------------
    # Validate only on held-out data.
    # -------------------------------------------------------------------------

    validate_candidates(
        zero_candidates,
        divisibility_candidates
    )

    # -------------------------------------------------------------------------
    # Exact representative cases.
    # -------------------------------------------------------------------------

    for n in [
        2, 3, 5, 7, 11,
        13, 17, 19, 23,
        29, 31
    ]:
        print_case(n)

    # -------------------------------------------------------------------------
    # The requested extremely large prime.
    # -------------------------------------------------------------------------

    large_test(BIG_N)

    # -------------------------------------------------------------------------
    # Large-number magnitude report.
    # -------------------------------------------------------------------------

    report_large_number(BIG_N)

    # -------------------------------------------------------------------------
    # Final state.
    # -------------------------------------------------------------------------

    print()
    print("=" * 78)
    print("VERSION 23 CONCLUSION")
    print("=" * 78)

    print()
    print("EXACT ORBIT:")
    print("    X_n = phi^n = F_(n-1) + F_n phi")

    print()
    print("EXACT NORM:")
    print("    N(X_n) = (-1)^n")

    print()
    print("EXACT TRANSLATION BALANCE:")
    print("    N(X_n+1)+N(X_n-1)")
    print("        = 2(1+(-1)^n)")

    print()
    print("EXACT LUCAS SEPARATION:")
    print("    N(X_n+1)-N(X_n-1)")
    print("        = 2L_n")

    print()
    print("LARGE-N RESULT:")
    print(f"    n = {BIG_N}")
    print(f"    prime = {is_prime(BIG_N)}")
    print(f"    N(phi^n) = {orbit_data(BIG_N)['N']}")
    print(
        f"    translation sum = "
        f"{orbit_data(BIG_N)['sum']}"
    )
    print(
        f"    identities = "
        f"{check_all_identities(BIG_N)}"
    )

    print()
    print("NO FLOATING-POINT PHI.")
    print("NO DECIMAL PHI.")
    print("NO GIANT INTEGER STRING CONVERSION.")
    print("NO STORED LARGE DECIMAL CONSTANT.")
    print()
    print("DONE.")


if __name__ == "__main__":
    main()