#!/usr/bin/env python3
# verify_hankel.py
# Ancillary file for the note "The Hankel determinants of an automatic sequence need not be
# automatic". David Niedbala Giraudin, ORCID 0009-0009-1526-1178.
# Technical assistance: Claude (Anthropic).
#
# Python 3 with numpy. Run:  python3 verify_hankel.py      (about a minute)
# Exact arithmetic modulo 3 only. Series in u = 1/T are arrays of coefficients mod 3.
# The script stops at the first failure and ends with "ALL CHECKS PASSED".

import sys
import numpy as np

N = 12000          # number of coefficients of alpha
HMAX = 300         # Hankel determinants computed up to this order
OMEGA_LEVEL = 16   # largest Omega_n built


def check(cond, msg):
    if not cond:
        print("FAILED: " + msg, flush=True)
        sys.exit(1)
    print("  ok  " + msg, flush=True)


def mul_trunc(a, b, n):
    """product of two series in u, truncated mod u**n, coefficients mod 3"""
    r = np.zeros(n, dtype=np.int64)
    nz = np.nonzero(b[:n])[0]
    for i in nz:
        r[i:] += b[i] * a[:n - i]
    return r % 3


print("Mills-Robbins quartic over F_3 and its Hankel determinants", flush=True)

# 1. alpha = u (1 + alpha^2 + alpha^4), coefficient by coefficient:
#    the coefficient of u^n of the right side only involves alpha_1..alpha_{n-1}.
alpha = np.zeros(N + 1, dtype=np.int64)       # alpha[n] = coefficient of u^n
sq = np.zeros(N + 1, dtype=np.int64)          # alpha^2
fo = np.zeros(N + 1, dtype=np.int64)          # alpha^4
for n in range(1, N + 1):
    m = n - 1
    c = (1 if m == 0 else 0) + sq[m] + fo[m]
    alpha[n] = c % 3
    # update alpha^2 at index n and alpha^4 at index n (they use alpha_1..alpha_{n-1}, alpha_n)
    s = int(np.dot(alpha[1:n], alpha[n - 1:0:-1])) if n >= 2 else 0
    sq[n] = s % 3
    f = int(np.dot(sq[1:n], sq[n - 1:0:-1])) if n >= 2 else 0
    fo[n] = f % 3
# residual of the equation multiplied by u^4: 1 + alpha^2 + alpha^4 - alpha/u, i.e.
# u (1 + alpha^2 + alpha^4) - alpha must vanish mod u^(N+1)
lhs = np.zeros(N + 1, dtype=np.int64)
lhs[1] = 1
lhs[1:] += sq[:N] + fo[:N]
check(np.all((lhs - alpha) % 3 == 0) and alpha[1] == 1,
      "alpha = sum alpha_n u^n satisfies alpha^4 + alpha^2 - T alpha + 1 = 0 to %d coefficients" % N)

# 2. the Buck-Robbins word Omega: entries are (sign, degree), all of the form +-T^(3^e)
def omega(n, memo={}):
    if n in memo:
        return memo[n]
    if n == 0:
        w = []
    elif n == 1:
        w = [(1, 1)]
    else:
        a, b = omega(n - 1), omega(n - 2)
        w = a + [(2, 1)] + [(s, 3 * d) for s, d in b] + [(2, 1)] + a   # cube: (cT^d)^3 = c T^(3d)
    memo[n] = w
    return w


W = omega(OMEGA_LEVEL)
D = [sum(d for _, d in omega(n)) for n in range(OMEGA_LEVEL + 1)]
Mx = [max([d for _, d in omega(n)] or [0]) for n in range(OMEGA_LEVEL + 1)]
check(all(D[n] == (3 ** n - 1) // 2 for n in range(1, OMEGA_LEVEL + 1))
      and all(Mx[n] == 3 ** ((n - 1) // 2) for n in range(1, OMEGA_LEVEL + 1)),
      "Lemma 3.1: D_n = (3^n-1)/2 and M_n = 3^floor((n-1)/2) for n = 1..%d" % OMEGA_LEVEL)
check(len(W) == 1136688 and Mx[OMEGA_LEVEL] == 2187,
      "Omega_16 has 1136688 entries, of degrees 1, 3, 9, ..., 2187")

# 3. the convergents: r_j = Q_j alpha - P_j satisfies r_{j+1} = a_{j+1} r_j + r_{j-1},
#    r_{-1} = -1, r_0 = alpha. Since a_{j+1} = c T^m, multiplication by T^m = u^(-m) is a shift.
#    Claim: v_u(r_j) = d_{j+1} exactly, i.e. P_j/Q_j are the convergents with deg Q_j = d_j.
#    r_j is known modulo u^(N+1-d_j).
d = [0]
for _, m in W:
    d.append(d[-1] + m)
r_prev = np.zeros(N + 1, dtype=np.int64)
r_prev[0] = 2                                   # -1
r_cur = alpha.copy()
checked = 0
j = 0
while True:
    prec = N + 1 - d[j]
    if j + 1 >= len(d) or d[j + 1] >= prec:
        break
    nz = np.nonzero(r_cur[:prec] % 3)[0]
    assert len(nz) and nz[0] == d[j + 1], ("valuation", j)
    c, m = W[j]
    nxt = np.zeros(N + 1, dtype=np.int64)
    nxt[:prec - m] = (c * r_cur[m:prec] + r_prev[:prec - m]) % 3
    r_prev, r_cur = r_cur, nxt
    checked += 1
    j += 1
check(checked > 1000,
      "the first %d partial quotients of alpha are those of Omega (|Q_j alpha - P_j| = 3^-d_{j+1})"
      % checked)
K = checked

# 4. Hankel determinants of S = (alpha_{n+1})_{n>=0} over F_3, n = 0..HMAX
S = alpha[1:2 * HMAX + 2] % 3


def det_mod3(A):
    A = A.copy() % 3
    n = A.shape[0]
    det = 1
    for col in range(n):
        piv = np.nonzero(A[col:, col])[0]
        if len(piv) == 0:
            return 0
        p = col + piv[0]
        if p != col:
            A[[col, p]] = A[[p, col]]
            det = -det
        det = det * A[col, col] % 3
        inv = A[col, col]                     # 1 or 2, its own inverse mod 3
        rows = np.arange(col + 1, n)
        f = (A[rows, col] * inv) % 3
        A[rows] = (A[rows] - np.outer(f, A[col])) % 3
    return det % 3


degset = set(x for x in d if x <= HMAX)
bad = [n for n in range(HMAX + 1)
       if (n == 0 or det_mod3(np.array([[S[i + k] for k in range(n)] for i in range(n)])) != 0)
       != (n in degset)]
check(bad == [], "Fact 2.1: for n = 0..%d, H_n != 0 exactly when n = d_j" % HMAX)

# 5. the ratios deg a_{j+1} / d_j, and the bound of Lemma 3.1
ok = all(W[j][1] <= 2 * d[j] ** 0.5 for j in range(len(W)) if d[j] >= 50)
check(ok, "Lemma 3.1: deg a_{j+1} <= 2 sqrt(d_j) for every entry of Omega_16 with d_j >= 50")
print("  ratio max deg a_{j+1}/d_j over dyadic ranges of d_j (Omega_16):", flush=True)
lo = 1
while lo < d[len(W) - 1]:
    vals = [W[j][1] / d[j] for j in range(1, len(W)) if lo <= d[j] < 2 * lo]
    if vals:
        print("      d_j in [%d, %d): %.6f" % (lo, 2 * lo, max(vals)), flush=True)
    lo *= 2

print("ALL CHECKS PASSED", flush=True)