"""THE THIRD HOFSTADTER — walking the $10,000 sequence to test THE
ADDRESSING DIVIDE.

That bridge said: G(n)=n-G(G(n-1)) is TAME (a nested value, no index-
offset -> a floor of a line) while Q(n)=Q(n-Q(n-1))+Q(n-Q(n-2)) is WILD
(values used as index OFFSETS n-Q(...), a self-chosen address book,
totality unproved). Hofstadter's third recurrence puts BOTH moves in one
line, and Conway offered $10,000 to understand it:

    a(1)=a(2)=1,  a(n) = a( a(n-1) ) + a( n - a(n-1) ).

Term 1 looks up at position p = a(n-1) (a value used DIRECTLY as an index,
G's move). Term 2 looks up at position n - p (a value used as an OFFSET
from n, Q's move). The two positions are COMPLEMENTARY: p + (n-p) = n. This
tool asks, by contact, what that combination buys — and whether the Divide's
clean "direct=tame / offset=wild" reading survives it.

Measurements, each carved with its class:
  (1) TOTALITY, by construction (COMPUTED): 1 <= a(n) <= n for all n, so both
      indices p and n-p always land in [1, n-1]. Unlike Q, no abyss.
  (2) THE POWER-OF-TWO CLOCK (COMPUTED; classical exact): a(2^k) = 2^(k-1)
      exactly, so a(n)/n = 1/2 dead-on at every power of two — a clock that
      does NOT die (contrast Q's Quoting Clock, which broke at k=9).
  (3) THE SELF-CORRECTION (COMPUTED; limit a(n)/n->1/2 is PROVED, Conway/
      Mallows): the complementary lookups drag the ratio back to 1/2. We
      measure sup|a(n)/n - 1/2|, WHERE it peaks (between powers of two), and
      Conway's own question: the last n with |a(n)/n - 1/2| >= 1/20.
"""
import sys


def conway(N):
    a = [0, 1, 1]                       # a[1]=a[2]=1 (index 0 unused)
    ap = a.append
    for n in range(3, N + 1):
        an1 = a[n - 1]
        ap(a[an1] + a[n - an1])
    return a


def main():
    N = 1 << 22                          # ~4.19M
    a = conway(N)

    # (1) totality + index-safety, asserted live
    bad = [n for n in range(1, N + 1) if not (1 <= a[n] <= n)]
    print("(1) TOTALITY  1 <= a(n) <= n  for n=1..%d :  %s"
          % (N, "HOLDS (both lookups always in range)" if not bad else "FAILS %r" % bad[:5]))

    # (2) the power-of-two clock: a(2^k) == 2^(k-1)
    clock = []
    k = 1
    while (1 << k) <= N:
        clock.append((k, a[1 << k], 1 << (k - 1)))
        k += 1
    ok = all(v == e for _, v, e in clock)
    print("(2) CLOCK  a(2^k) == 2^(k-1) :  %s  (checked k=1..%d)"
          % ("EXACT, every rung" if ok else "BREAKS", clock[-1][0]))
    print("    sample: a(2^20)=%d  vs 2^19=%d" % (a[1 << 20], 1 << 19))

    # (3) self-correction. The GLOBAL sup is a trivial transient (a(1)/1=1,
    # dev 0.5); the meaningful quantity is the ASYMPTOTIC sup, so we measure
    # it past the young sequence (n >= 2^10) and, separately, Conway's 1/20 exit.
    ASYM = 1 << 10
    max_dev, at = 0.0, ASYM
    last_20 = 0                          # last n with |a(n)/n - 1/2| >= 0.05
    for n in range(1, N + 1):
        ad = abs(a[n] / n - 0.5)
        if n >= ASYM and ad > max_dev:
            max_dev, at = ad, n
        if ad >= 0.05:
            last_20 = n
    lo = 1 << (at.bit_length() - 1)
    print("(3) SELF-CORRECTION toward 1/2:")
    print("    (global sup is the trivial transient 0.5 at n=1; the real measure is asymptotic)")
    print("    asymptotic sup |a(n)/n - 1/2| over [2^10, %d] :  %.6f  at n=%d" % (N, max_dev, at))
    print("        (that peak sits %.2f of the way through its octave 2^%d..2^%d)"
          % ((at - lo) / lo, at.bit_length() - 1, at.bit_length()))
    print("    CONWAY'S $10,000 VALUE (COMPUTED): a(n)/n stays within 1/20 of 1/2 for ALL n > %d" % last_20)
    print("        (|a(%d)/%d - 1/2| = %.5f was the last violation; Mallows 1991 = 1489)"
          % (last_20, last_20, abs(a[last_20] / last_20 - 0.5)))

    # the deviation is a fractal that RESETS at each power of two (clock (2)):
    # report the peak dev inside each octave to show it does not grow.
    print("    peak |a(n)/n-1/2| inside each octave [2^k, 2^(k+1)) — bounded, not growing:")
    k = 4
    while (1 << (k + 1)) <= N:
        seg = range(1 << k, 1 << (k + 1))
        pk = max(abs(a[n] / n - 0.5) for n in seg)
        if k in (4, 8, 12, 16, 20) or (1 << (k + 2)) > N:
            print("      octave 2^%2d..2^%2d : peak %.6f" % (k, k + 1, pk))
        k += 1


if __name__ == "__main__":
    sys.setrecursionlimit(10000)
    main()