#!/usr/bin/env python3
"""
============================================================================
COMPUTATIONAL VERIFICATION SCRIPT FOR "META-AUTOMATIC SEQUENCES"
============================================================================

This single self-contained Python 3 file verifies every computational claim
of the paper "Meta-Automatic Sequences". It has no external dependencies
beyond the Python 3 standard library (tested with Python 3.10 and 3.12).

USAGE
-----

    python3 verify_meta_automatic.py

Optional flags:

    --N             number of sequence terms to generate (default 300000)
    --max-factor    largest factor length to test (default 64)
    --morphism-level   check morphisms through 4^level terms (default 8)

WHAT IS VERIFIED
----------------

  Section 1.  Sequence generation from the defining recurrences (Q, M1, M2).
  Section 2.  Base-4 MSB-first DFAOs of Q, M1, M2.
              Includes Moore minimization to confirm minimal state count
              (3 for Q, 4 for M1, 4 for M2).
  Section 3.  4-uniform morphisms for M1 and M2.
              Fixed point of each morphism, after coding, agrees with the
              generated sequence through 4^level terms.
  Section 4.  Closed digit formula M2(n) = t(q(n)), where t is Thue-Morse
              and q is the bit-masking operator that zeroes binary positions
              m with (m mod 3) = 2.
  Section 5.  Factor complexity values p(1..15) for Q, M1, M2.
  Section 6.  Piecewise factor-complexity formula for Q (Proposition
              Q-complete).
  Section 7.  Piecewise factor-complexity formula for M2 (Proposition
              M2-complete).
  Section 8.  Right-special factor count values for Q and M2 (the sets are
              {2, 3, 4} and {2, 4} respectively).
  Section 9.  Reproduce the Appendix table (Table 1). For each of the 25
              specializations of the base-4 skeleton, compute the size of the
              minimal MSB-first DFAO via signature equivalence on indices.
              The script prints the resulting table alongside the values
              that appear in the paper.

Each section prints a PASS/FAIL summary line.
"""

from __future__ import annotations

import argparse
from typing import Dict, List, Sequence, Tuple

BitSeq = List[int]


# ============================================================================
# SECTION 1.  SEQUENCE GENERATION
# ============================================================================
#
# The three sequences Q, M1, M2 are defined in the paper by the recurrences
# below.  We generate the first N terms of each.
# ----------------------------------------------------------------------------


def generate_Q(N: int) -> BitSeq:
    """Generate the first N terms of the Thue-Morse Quarto sequence Q."""
    a = [0] * N
    if N > 1:
        a[1] = 1
    for m in range(2, N):
        n, r = divmod(m, 4)
        if r == 0:
            a[m] = a[n]
        elif r == 1:
            a[m] = 1 - a[n]
        elif r == 2:
            a[m] = a[2 * n]
        else:
            a[m] = 1 - a[2 * n]
    return a


def generate_M1(N: int) -> BitSeq:
    """Generate the first N terms of M1.

    Defining recurrence:
        M1(0)=0, M1(1)=1.
        For n >= 1:   M1(4n)   = M1(n)
                      M1(4n+1) = 1 - M1(n)
        For n >= 0:   M1(4n+2) = M1(2n + 1 - M1(n))
                      M1(4n+3) = 1 - M1(2n + 1 - M1(n))
    """
    a = [0] * N
    if N > 1:
        a[1] = 1
    for m in range(2, N):
        n, r = divmod(m, 4)
        if r == 0:
            a[m] = a[n]
        elif r == 1:
            a[m] = 1 - a[n]
        elif r == 2:
            a[m] = a[2 * n + 1 - a[n]]
        else:
            a[m] = 1 - a[2 * n + 1 - a[n]]
    return a


def generate_M2(N: int) -> BitSeq:
    """Generate the first N terms of M2.

    Defining recurrence:
        M2(0)=0, M2(1)=1.
        For n >= 1:   M2(4n)   = M2(2n + M2(n))
                      M2(4n+1) = 1 - M2(4n)
        For n >= 0:   M2(4n+2) = M2(2n + 1 - M2(n))
                      M2(4n+3) = 1 - M2(4n+2)
    """
    a = [0] * N
    if N > 1:
        a[1] = 1
    for m in range(2, N):
        n, r = divmod(m, 4)
        if r == 0:
            a[m] = a[2 * n + a[n]]
        elif r == 1:
            a[m] = 1 - a[4 * n]
        elif r == 2:
            a[m] = a[2 * n + 1 - a[n]]
        else:
            a[m] = 1 - a[4 * n + 2]
    return a


# ============================================================================
# SECTION 2.  DFAO TABLES AND DFAO VERIFICATION
# ============================================================================
#
# We encode each base-4 MSB-first DFAO from the paper as a transition table
# and an output map.  We then check, on the first N generated terms, that
# evaluating the DFAO gives the same sequence, and that Moore minimization
# returns the announced number of states (3, 4, 4).
# ----------------------------------------------------------------------------


Q_TABLE: List[List[int]] = [
    [0, 2, 1, 2],   # state 0, on digits 0,1,2,3
    [1, 2, 2, 1],
    [2, 1, 1, 2],
]
Q_OUTPUT: List[int] = [0, 0, 1]

M1_TABLE: List[List[int]] = [
    [0, 1, 1, 2],
    [1, 2, 3, 0],
    [2, 3, 0, 1],
    [3, 0, 2, 3],
]
M1_OUTPUT: List[int] = [0, 1, 0, 1]

M2_TABLE: List[List[int]] = [
    [0, 1, 3, 2],
    [2, 3, 1, 0],
    [3, 2, 0, 1],
    [1, 0, 2, 3],
]
M2_OUTPUT: List[int] = [0, 1, 0, 1]


def base_digits_msd(n: int, base: int) -> List[int]:
    """Return the most-significant-digit-first base-b representation of n."""
    if n < 0:
        raise ValueError("n must be non-negative")
    if n == 0:
        return [0]
    digs: List[int] = []
    while n:
        digs.append(n % base)
        n //= base
    return digs[::-1]


def eval_dfao(n: int, table: Sequence[Sequence[int]], output: Sequence[int],
              base: int = 4) -> int:
    """Evaluate the DFAO at index n in MSB-first base-b mode."""
    state = 0
    for d in base_digits_msd(n, base):
        state = table[state][d]
    return output[state]


def minimize_moore(table: Sequence[Sequence[int]],
                   output: Sequence[int]) -> Tuple[int, List[int]]:
    """Moore minimization. Return (num_classes, equivalence_class_of_each_state)."""
    nstates = len(table)
    cls: List[int] = []
    partition: Dict[Tuple[int, ...], int] = {}
    for o in output:
        key = (o,)
        if key not in partition:
            partition[key] = len(partition)
        cls.append(partition[key])

    changed = True
    while changed:
        changed = False
        new_partition: Dict[Tuple[int, Tuple[int, ...]], int] = {}
        new_cls = [0] * nstates
        for q in range(nstates):
            key = (output[q], tuple(cls[table[q][d]] for d in range(len(table[q]))))
            if key not in new_partition:
                new_partition[key] = len(new_partition)
            new_cls[q] = new_partition[key]
        if new_cls != cls:
            cls = new_cls
            changed = True
    return len(set(cls)), cls


def check_dfao(name: str, seq: BitSeq, table: Sequence[Sequence[int]],
               output: Sequence[int]) -> None:
    bad = [n for n, value in enumerate(seq) if value != eval_dfao(n, table, output, 4)]
    min_states, _ = minimize_moore(table, output)
    status = "PASS" if not bad else "FAIL"
    print(f"  [{status}] {name}: DFAO matches sequence on {len(seq)} terms; "
          f"minimized states = {min_states}")
    if bad:
        print(f"         first bad indices: {bad[:10]}")


# ============================================================================
# SECTION 3.  4-UNIFORM MORPHISMS
# ============================================================================
#
# Each of M1 and M2 is obtained as a coding of the fixed point of an explicit
# 4-uniform morphism starting from the letter "0".  We iterate the morphism
# `levels` times and check that the coded fixed point agrees with the
# generated sequence on the resulting prefix of length 4^levels.
# ----------------------------------------------------------------------------


M1_MORPHISM: Dict[str, str] = {
    "0": "0112",
    "1": "1230",
    "2": "2301",
    "3": "3023",
}

M2_MORPHISM: Dict[str, str] = {
    "0": "0132",
    "1": "2310",
    "2": "3201",
    "3": "1023",
}

# Both M1 and M2 use the same coding: 0,2 -> 0 and 1,3 -> 1.
MORPHIC_OUTPUT: Dict[str, str] = {"0": "0", "1": "1", "2": "0", "3": "1"}


def iterate_morphism(morphism: Dict[str, str], start: str, levels: int) -> str:
    w = start
    for _ in range(levels):
        w = "".join(morphism[c] for c in w)
    return w


def check_morphism(name: str, seq: BitSeq, morphism: Dict[str, str],
                   levels: int) -> None:
    state_word = iterate_morphism(morphism, "0", levels)
    coded = "".join(MORPHIC_OUTPUT[c] for c in state_word)
    target = "".join(str(x) for x in seq[: len(coded)])
    ok = coded == target
    status = "PASS" if ok else "FAIL"
    print(f"  [{status}] {name}: morphism agrees with sequence "
          f"through 4^{levels} = {len(coded)} terms")


# ============================================================================
# SECTION 4.  CLOSED DIGIT FORMULA  M2(n) = t(q(n))
# ============================================================================
#
# The paper proves M2(n) = t(q(n)), where t is the Thue-Morse sequence and
# q(n) is obtained from the binary expansion of n by zeroing every bit whose
# position m satisfies (m mod 3) == 2.  Equivalently, in base 4, the
# operation acts on each digit d_k according to (k mod 3).
#
# We verify both forms of the formula numerically.
# ----------------------------------------------------------------------------


def thue_morse(n: int) -> int:
    return n.bit_count() & 1


def q_mask(n: int) -> int:
    """Zero the bits in positions m with (m mod 3) == 2."""
    ans = 0
    bit = 1
    pos = 0
    while bit <= n:
        if (n & bit) and (pos % 3 != 2):
            ans |= bit
        bit <<= 1
        pos += 1
    return ans


def m2_closed_digit_formula(n: int) -> int:
    """Base-4 form of the closed formula from Theorem (M2-tq)."""
    ans = 0
    k = 0
    x = n
    while x:
        d = x & 3
        i = d & 1
        j = (d >> 1) & 1
        if k % 3 == 0:
            ans ^= i ^ j
        elif k % 3 == 1:
            ans ^= j
        else:
            ans ^= i
        x >>= 2
        k += 1
    return ans


def check_m2_formula(seq: BitSeq) -> None:
    bad = []
    for n, value in enumerate(seq):
        if value != thue_morse(q_mask(n)) or value != m2_closed_digit_formula(n):
            bad.append(n)
            if len(bad) >= 10:
                break
    status = "PASS" if not bad else "FAIL"
    print(f"  [{status}] M2(n) = t(q(n)) and base-4 digit formula agree "
          f"with sequence on {len(seq)} terms")
    if bad:
        print(f"         first bad indices: {bad}")


# ============================================================================
# SECTION 5.  INITIAL FACTOR-COMPLEXITY VALUES
# ============================================================================
#
# We compute p(1..15) for each of Q, M1, M2 from a long prefix of the
# sequence.  The values listed in the paper are:
#
#   p_Q (1..15) = (2, 4, 6, 10, 12, 16, 20, 23, 26, 30, 34, 38, 42, 45, 48)
#   p_M1(1..15) = (2, 4, 6, 9,  12, 17, 22, 28, 34, 40, 46, 52, 58, 64, 70)
#   p_M2(1..15) = (2, 4, 6, 10, 12, 14, 16, 18, 20, 24, 28, 32, 36, 40, 44)
# ----------------------------------------------------------------------------


P_Q_REFERENCE  = [2, 4, 6, 10, 12, 16, 20, 23, 26, 30, 34, 38, 42, 45, 48]
P_M1_REFERENCE = [2, 4, 6,  9, 12, 17, 22, 28, 34, 40, 46, 52, 58, 64, 70]
P_M2_REFERENCE = [2, 4, 6, 10, 12, 14, 16, 18, 20, 24, 28, 32, 36, 40, 44]


def factor_complexity(seq: BitSeq, max_len: int) -> List[int]:
    """Return [p(1), p(2), ..., p(max_len)] for the binary word seq."""
    s = "".join(str(x) for x in seq)
    return [len({s[i : i + n] for i in range(0, len(s) - n + 1)})
            for n in range(1, max_len + 1)]


def check_initial_complexities(Q: BitSeq, M1: BitSeq, M2: BitSeq) -> None:
    pQ  = factor_complexity(Q,  15)
    pM1 = factor_complexity(M1, 15)
    pM2 = factor_complexity(M2, 15)
    okQ  = pQ  == P_Q_REFERENCE
    okM1 = pM1 == P_M1_REFERENCE
    okM2 = pM2 == P_M2_REFERENCE
    print(f"  [{'PASS' if okQ  else 'FAIL'}] p_Q (1..15) = {pQ}")
    print(f"  [{'PASS' if okM1 else 'FAIL'}] p_M1(1..15) = {pM1}")
    print(f"  [{'PASS' if okM2 else 'FAIL'}] p_M2(1..15) = {pM2}")


# ============================================================================
# SECTION 6.  PIECEWISE FORMULA FOR p_Q   (Proposition Q-complete)
# ============================================================================
#
# For n >= 16, write n = 2^k + i with k = floor(log_2 n) and 0 <= i < 2^k.
# Then
#         13 * 2^(k-2) - 2          if i = 0,
# p_Q(n) = 4n - 3*2^(k-2) - 4       if 1 <= i <= 2^(k-1),
#         3n + 3*2^(k-2) - 3        if 2^(k-1) < i <= 3*2^(k-2),
#         2n + 5*2^(k-1) - 2        if 3*2^(k-2) < i < 2^k.
# ----------------------------------------------------------------------------


def p_Q_formula(n: int) -> int:
    if n < 16:
        raise ValueError("Formula stated for n >= 16")
    k = n.bit_length() - 1
    P = 1 << k
    i = n - P
    if i == 0:
        return 13 * (1 << (k - 2)) - 2
    if 1 <= i <= (1 << (k - 1)):
        return 4 * n - 3 * (1 << (k - 2)) - 4
    if (1 << (k - 1)) < i <= 3 * (1 << (k - 2)):
        return 3 * n + 3 * (1 << (k - 2)) - 3
    return 2 * n + 5 * (1 << (k - 1)) - 2


def check_p_Q_formula(Q: BitSeq, max_factor: int) -> None:
    pQ = factor_complexity(Q, max_factor)
    bad = [n for n in range(16, max_factor + 1) if pQ[n - 1] != p_Q_formula(n)]
    status = "PASS" if not bad else "FAIL"
    print(f"  [{status}] Proposition Q-complete: p_Q(n) matches formula "
          f"for 16 <= n <= {max_factor}")
    if bad:
        print(f"         first bad n: {bad[:10]}")


# ============================================================================
# SECTION 7.  PIECEWISE FORMULA FOR p_M2  (Proposition M2-complete)
# ============================================================================
#
# We have p_M2(1) = 2.  For n >= 2, write k = floor(log_2 n), P = 2^k, and
# r = k mod 3.  The piecewise formula then depends on r and on whether n is
# in [P, P+1] or in a later sub-interval.
# ----------------------------------------------------------------------------


def p_M2_formula(n: int) -> int:
    if n < 1:
        raise ValueError("n must be positive")
    if n == 1:
        return 2
    k = n.bit_length() - 1
    P = 1 << k
    r = k % 3
    if r == 2:
        return 2 * n + (P - 2)
    if r == 0:
        if P <= n <= P + 1:
            return 2 * n + ((P >> 1) - 2)
        if P + 1 <= n <= 7 * (P >> 2) + 1:
            return 4 * n - (3 * (P >> 1) + 4)
        return 2 * n + ((P << 1) - 2)
    # r == 1
    if P <= n <= P + 1:
        return 2 * n + (P - 2)
    if P + 1 <= n <= 3 * (P >> 1) + 1:
        return 4 * n - (P + 4)
    return 2 * n + ((P << 1) - 2)


def check_p_M2_formula(M2: BitSeq, max_factor: int) -> None:
    pM2 = factor_complexity(M2, max_factor)
    bad = [n for n in range(1, max_factor + 1) if pM2[n - 1] != p_M2_formula(n)]
    status = "PASS" if not bad else "FAIL"
    print(f"  [{status}] Proposition M2-complete: p_M2(n) matches formula "
          f"for 1 <= n <= {max_factor}")
    if bad:
        print(f"         first bad n: {bad[:10]}")


# ============================================================================
# SECTION 8.  RIGHT-SPECIAL FACTOR COUNTS
# ============================================================================
#
# A length-n factor u of a is right-special if both u0 and u1 are factors of
# a.  The right-special-factor count RS_a(n) equals p_a(n+1) - p_a(n).
#
# The paper proves:
#   RS_Q (n) in {2, 3, 4} for all n >= 1
#   RS_M2(n) in {2, 4}    for all n >= 1
# ----------------------------------------------------------------------------


def right_special_count(seq: BitSeq, n: int) -> int:
    s = "".join(str(x) for x in seq)
    followers: Dict[str, set] = {}
    for i in range(0, len(s) - n):
        f = s[i : i + n]
        followers.setdefault(f, set()).add(s[i + n])
    return sum(1 for vals in followers.values() if len(vals) >= 2)


def check_rs_counts(Q: BitSeq, M2: BitSeq, max_rs: int) -> None:
    rs_q  = {right_special_count(Q,  n) for n in range(1, max_rs + 1)}
    rs_m2 = {right_special_count(M2, n) for n in range(1, max_rs + 1)}
    okQ  = rs_q  <= {2, 3, 4} and rs_q  == {2, 3, 4}
    okM2 = rs_m2 <= {2, 4}    and rs_m2 == {2, 4}
    print(f"  [{'PASS' if okQ  else 'FAIL'}] RS_Q (n) values for n <= {max_rs}: "
          f"{sorted(rs_q)}")
    print(f"  [{'PASS' if okM2 else 'FAIL'}] RS_M2(n) values for n <= {max_rs}: "
          f"{sorted(rs_m2)}")


# ============================================================================
# SECTION 9.  REPRODUCE THE APPENDIX TABLE
# ============================================================================
#
# The Appendix of the paper lists the minimal base-4 MSB-first DFAO size for
# 25 specializations of the balanced skeleton
#       a(4n+0) = F(n),         a(4n+1) = 1 - F(n),
#       a(4n+2) = G(n),         a(4n+3) = 1 - G(n).
# Here F and G each range over five basic choices:
#       a(n),  a(2n),  a(2n+1),  a(2n + a(n)),  a(2n + 1 - a(n)).
# The first three are "non-nested" and the last two contain a nested call.
#
# For each (F, G) we generate a long prefix of the resulting sequence, detect
# ultimate periodicity, and otherwise estimate the minimal MSB-first DFAO
# size by counting the distinct depth-d future-signatures of indices.
#
# A future-signature of index n at depth d is the tuple
#       (a(base^L * n + r)) for 0 <= L <= d and 0 <= r < base^L.
# Two indices yielding equal signatures lie in the same state of the minimal
# MSB-first DFAO. With d=4 and base=4, signatures have length 341 and reliably
# distinguish small state classes.
# ----------------------------------------------------------------------------


# Five basic selectors. The third element of each tuple is the "nested" flag:
# 1 if the selector contains a nested call a(... + a(n)), 0 otherwise.
APPENDIX_SELECTORS: List[Tuple[str, "callable", int]] = [
    ("a(n)",         lambda a, n: n,             0),
    ("a(2n)",        lambda a, n: 2 * n,         0),
    ("a(2n+1)",      lambda a, n: 2 * n + 1,     0),
    ("a(2n+a(n))",   lambda a, n: 2 * n + a[n],  1),
    ("a(2n+1-a(n))", lambda a, n: 2 * n + 1 - a[n], 1),
]


def gen_appendix_sequence(F_idx, G_idx, N: int) -> BitSeq:
    """Generate a prefix of the sequence defined by F (branches 0, 1) and
    G (branches 2, 3) of the base-4 balanced skeleton.
    Convention: a(0)=0, a(1)=1; recurrence rules are applied for m >= 2,
    so the initial values are never overwritten."""
    a = [0] * N
    if N > 1:
        a[1] = 1
    for m in range(2, N):
        n, r = divmod(m, 4)
        if r == 0:
            idx = F_idx(a, n)
            a[m] = a[idx] if 0 <= idx < m else 0
        elif r == 1:
            idx = F_idx(a, n)
            a[m] = 1 - a[idx] if 0 <= idx < m else 0
        elif r == 2:
            idx = G_idx(a, n)
            a[m] = a[idx] if 0 <= idx < m else 0
        else:
            idx = G_idx(a, n)
            a[m] = 1 - a[idx] if 0 <= idx < m else 0
    return a


def msb_signature(seq: BitSeq, n: int, depth: int, base: int = 4) -> tuple:
    """Future-signature of index n at depth d, used to identify MSB-first DFAO states."""
    sig = []
    for L in range(depth + 1):
        step = base ** L
        for r in range(step):
            sig.append(seq[step * n + r])
    return tuple(sig)


def min_msb_dfao_size(seq: BitSeq, base: int = 4, depth: int = 4):
    """Estimate the size of the minimal MSB-first DFAO from a prefix of seq."""
    N = len(seq)
    step = base ** depth
    M = (N - step) // step
    if M < 4:
        return None
    sigs = {msb_signature(seq, n, depth, base) for n in range(M + 1)}
    return len(sigs)


def is_ultimately_periodic(seq: BitSeq, max_period: int = 200) -> bool:
    """Detect ultimate periodicity with a moderate prefix and period bound."""
    N = len(seq)
    for prefix in range(min(100, N // 3)):
        for period in range(1, max_period + 1):
            if N - prefix < 4 * period:
                break
            tail = seq[prefix:]
            if all(tail[i] == tail[i + period] for i in range(len(tail) - period)):
                return True
    return False


def reproduce_appendix_table(N_terms: int = 20000, depth: int = 4) -> None:
    """Print the Appendix table of minimal DFAO sizes for the 25 (F, G) pairs."""
    print(f"{'F(n)':<15} | {'G(n)':<15} | {'DFAO':<6} | {'Nested':<6}")
    print("-" * 50)
    for F_name, F_idx, F_nest in APPENDIX_SELECTORS:
        for G_name, G_idx, G_nest in APPENDIX_SELECTORS:
            seq = gen_appendix_sequence(F_idx, G_idx, N_terms)
            nested = F_nest + G_nest
            if is_ultimately_periodic(seq):
                dfao_str = "---"
            else:
                sz = min_msb_dfao_size(seq, depth=depth)
                dfao_str = str(sz) if sz is not None else "?"
            print(f"{F_name:<15} | {G_name:<15} | {dfao_str:<6} | {nested:<6}")


# ============================================================================
# MAIN
# ============================================================================


def main() -> None:
    parser = argparse.ArgumentParser(description=__doc__,
                                     formatter_class=argparse.RawDescriptionHelpFormatter)
    parser.add_argument("--N", type=int, default=300_000,
                        help="number of sequence terms to generate (default 300000)")
    parser.add_argument("--max-factor", type=int, default=64,
                        help="largest factor length to test (default 64)")
    parser.add_argument("--morphism-level", type=int, default=8,
                        help="check morphisms through 4^level terms (default 8)")
    args = parser.parse_args()

    needed = max(args.N, 10 * args.max_factor + 10, 4 ** args.morphism_level)

    print("=" * 72)
    print("Verification script for 'Meta-Automatic Sequences'")
    print("=" * 72)
    print(f"Generating {needed} terms of Q, M1, M2...")
    Q  = generate_Q (needed)
    M1 = generate_M1(needed)
    M2 = generate_M2(needed)

    print()
    print("Section 2: DFAO checks")
    print("-" * 72)
    check_dfao("Q ", Q [: args.N], Q_TABLE,  Q_OUTPUT)
    check_dfao("M1", M1[: args.N], M1_TABLE, M1_OUTPUT)
    check_dfao("M2", M2[: args.N], M2_TABLE, M2_OUTPUT)

    print()
    print("Section 3: Morphism checks")
    print("-" * 72)
    check_morphism("M1", M1, M1_MORPHISM, args.morphism_level)
    check_morphism("M2", M2, M2_MORPHISM, args.morphism_level)

    print()
    print("Section 4: Closed digit formula for M2")
    print("-" * 72)
    check_m2_formula(M2[: args.N])

    print()
    print("Section 5: Initial factor-complexity values")
    print("-" * 72)
    check_initial_complexities(Q, M1, M2)

    print()
    print("Section 6: Piecewise formula for p_Q")
    print("-" * 72)
    check_p_Q_formula(Q, args.max_factor)

    print()
    print("Section 7: Piecewise formula for p_M2")
    print("-" * 72)
    check_p_M2_formula(M2, args.max_factor)

    print()
    print("Section 8: Right-special factor counts")
    print("-" * 72)
    max_rs = min(args.max_factor - 1, 128)
    check_rs_counts(Q, M2, max_rs)

    print()
    print("Section 9: Appendix table (minimal MSB-first DFAO size for the 25")
    print("           specializations of the base-4 skeleton)")
    print("-" * 72)
    reproduce_appendix_table(N_terms=20000, depth=4)

    print()
    print("=" * 72)
    print("End of verification.")
    print("=" * 72)


if __name__ == "__main__":
    main()
