"""THE GOLDEN TWIN — Hofstadter's tame sibling, walked by contact.

The country's Q region stands over an abyss: Q(n) = Q(n-Q(n-1)) + Q(n-Q(n-2))
is not even proved total, and every Q landmark (THE SPIKE'S BITE, THE FORBIDDEN
ADDRESS) turns on Q *addressing itself* -- the previous values are used as INDEX
OFFSETS (n - Q(...)), so the recursion reaches into a self-chosen address book.

Hofstadter wrote a milder twin the same year:
    G(0) = 0,   G(n) = n - G(G(n-1)).
Here the previous value is NOT an address -- it is fed straight back INTO G,
nested as a VALUE, never used to index. This tool asks, by contact, what that
one difference buys.

Three measurements, each carved with its evidence class:
  (1) EXACT closed form. Claim: G(n) = floor((n+1)/phi), phi=(1+sqrt5)/2.
      Verified with NO floating point -- pure big-integer isqrt -- so the
      verdict is exact, not "close". Derivation of the integer form is in
      _exact_floor() below (checkable at the trailhead, per law 9).
  (2) The float LIE. floor((n+1)/phi) computed in float64 eventually disagrees
      with the exact form near an integer boundary. We find the first n where
      it lies -- a deliberate rhyme with THE STAIRCASE OF e ("float64 lies from
      term 20"): exactness beats float here too.
  (3) The contrast, live. Q's error e(n) = Q(n) - n/2 is measured over the same
      horizon: G's error is bounded and periodic-ish (it IS a Beatty/Wythoff
      staircase), Q's wanders. Same inventor, one word of difference, two fates.
"""
import math


def build_G(N):
    """G[0..N] by the recursion itself, iteratively (no recursion-depth wall)."""
    G = [0] * (N + 1)
    for n in range(1, N + 1):
        G[n] = n - G[G[n - 1]]
    return G


def _exact_floor(n):
    """floor((n+1)*(sqrt5 - 1)/2) with exact integer arithmetic.

    Let a = n+1. The target is (a*sqrt5 - a)/2 = (sqrt(5a^2) - a)/2.
    5a^2 is never a perfect square for a>=1 (sqrt5 irrational), so
    s = isqrt(5a^2) = floor(sqrt(5a^2)) and sqrt(5a^2) = s + f, 0<f<1.
    floor((s - a + f)/2): with d = s - a an integer and 0<f<1,
      d even -> (d + f)/2 has floor d/2;   d odd -> floor (d-1)/2.
    Both are exactly Python's d//2. Hence the one-liner below is EXACT.
    """
    a = n + 1
    s = math.isqrt(5 * a * a)
    return (s - a) // 2


def main():
    N = 2_000_000
    G = build_G(N)

    # (1) EXACT closed form -- zero mismatches expected.
    bad = [n for n in range(N + 1) if G[n] != _exact_floor(n)]
    print("(1) EXACT closed form  G(n) == floor((n+1)/phi)  [integer isqrt]")
    print("    checked n = 0..%d :  mismatches = %d  %s"
          % (N, len(bad), "(EXACT MATCH)" if not bad else ("first: %r" % bad[:5])))

    # totality-of-a-different-kind: G is manifestly total (each step defined),
    # but let's confirm the self-index G[n-1] never runs off the array -- i.e.
    # G(n) <= n always (so G[G[n-1]] is always a computed cell).
    off = [n for n in range(1, N + 1) if G[n] > n]
    print("    self-index safe (G(n)<=n) for all n:  %s" % ("YES" if not off else "NO %r" % off[:5]))

    # (2) the float lie -- and WHY it comes so late. phi is the worst-
    # approximable number (CF = [1;1,1,1,...]), so (n+1)/phi stays farther from
    # integer boundaries than for any other constant: ||a/phi|| >= 1/(sqrt5*a).
    # float64's absolute error in a/phi is ~ eps*(a/phi) ~ 1.4e-16*a, so the lie
    # cannot appear until 1.4e-16*a ~ 1/(sqrt5*a), i.e. a ~ 5.7e7. We scan (no
    # array needed) for the actual first witness. Inverse rhyme with THE
    # STAIRCASE OF e: e's CF has LARGE terms so float lies at term 20; phi's CF
    # is all 1s so float holds longest. The continued fraction governs the lie.
    phi = (1 + 5 ** 0.5) / 2
    first_lie, cap = None, 150_000_000
    n = 0
    while n <= cap:
        if int((n + 1) / phi) != _exact_floor(n):
            first_lie = n
            break
        n += 1
    if first_lie is None:
        print("(2) float64 floor((n+1)/phi): still no lie through n=%d" % cap)
    else:
        a = first_lie + 1
        print("(2) float64 floor((n+1)/phi) FIRST LIES at n=%d  (~%.2eM)"
              % (first_lie, first_lie / 1e6))
        print("    float says %d, truth is %d;  boundary distance ||a/phi|| = %.3e  (bound 1/(sqrt5*a)=%.3e)"
              % (int(a / phi), _exact_floor(first_lie),
                 abs(a / phi - round(a / phi)), 1 / (5 ** 0.5 * a)))

    # (3) the contrast, live: Q's error vs G's error over the same horizon.
    M = 1_000_000
    Q = [0, 1, 1]  # Q[1]=Q[2]=1
    total_ok = True
    for n in range(3, M + 1):
        v = Q[n - Q[n - 1]] + Q[n - Q[n - 2]]
        Q.append(v)
        if n - Q[n - 1] < 1 or n - Q[n - 2] < 1:
            total_ok = False
    # error signals -- each centered on the line it actually hugs.
    # G rides n/phi (Beatty), so |G(n)-(n+1)/phi| < 1 BY THE FLOOR: pinned.
    # Q rides n/2, but its deviation WANDERS with growing amplitude.
    gdev = [abs(G[n] - (n + 1) / phi) for n in range(M + 1)]
    qerr = [Q[n] - n / 2 for n in range(1, M + 1)]
    print("(3) same recursion family, one word apart:")
    print("    |G(n)-(n+1)/phi| over 0..%d :  max %.6f  (< 1 always -- PINNED to a straight line)"
          % (M, max(gdev)))
    print("    Q(n)-n/2         over 1..%d :  min %+.1f  max %+.1f  (amplitude ~%.0f -- WANDERS off the line)"
          % (M, min(qerr), max(qerr), max(abs(min(qerr)), abs(max(qerr)))))
    print("    Q totality holds over 1..%d :  %s" % (M, "YES" if total_ok else "NO"))

    # the phi-density signature, exact: #{n<=N : G(n)=G(n+1)} / N -> 2 - phi = 1/phi^2
    flats = sum(1 for n in range(N) if G[n] == G[n + 1])
    print("    G's 'flat step' density (G(n)=G(n+1)) = %.6f   vs  2-phi = %.6f"
          % (flats / N, 2 - phi))


if __name__ == "__main__":
    main()