//! Independent exact checker for finite irreducible covering-set witnesses.
//!
//! This checker deliberately shares no code with the search programs under
//! `workspace/`.  It accepts only instances whose common period is small enough
//! to inspect exhaustively.  Refusing an oversized or over-budget instance is a
//! verification failure, never a PASS.

use rayon::prelude::*;
use std::collections::BTreeSet;
use std::env;
use std::fs;
use std::process;

const MAX_K: usize = 128;
const MAX_PERIOD: u64 = 20_000_000;
const DEFAULT_MAX_SEARCH_NODES: u64 = 50_000_000;

#[derive(Clone, Copy, Debug)]
struct Congruence {
    modulus: u64,
    residue: u64,
}

#[derive(Debug)]
enum SearchOutcome {
    Cover(Vec<u64>),
    NoCover,
    BudgetExceeded,
}

fn gcd(mut a: u64, mut b: u64) -> u64 {
    while b != 0 {
        let r = a % b;
        a = b;
        b = r;
    }
    a
}

fn gcd_u128(mut a: u128, mut b: u128) -> u128 {
    while b != 0 {
        let r = a % b;
        a = b;
        b = r;
    }
    a
}

/// Compute the exact reciprocal sum as a reduced u128 fraction.  `None` means
/// only that this auxiliary representation overflowed.
fn reciprocal_sum_fraction(moduli: &[u64]) -> Option<(u128, u128)> {
    let mut numerator = 0_u128;
    let mut denominator = 1_u128;
    for &n64 in moduli {
        let n = u128::from(n64);
        let common = gcd_u128(denominator, n);
        let n_over_common = n / common;
        let denominator_over_common = denominator / common;
        numerator = numerator
            .checked_mul(n_over_common)?
            .checked_add(denominator_over_common)?;
        denominator = denominator.checked_mul(n_over_common)?;
        let reduction = gcd_u128(numerator, denominator);
        numerator /= reduction;
        denominator /= reduction;
    }
    Some((numerator, denominator))
}

/// Return whether the exact reciprocal sum is below one.  `None` performs no
/// pruning, so auxiliary overflow can never create a false non-cover answer.
fn reciprocal_sum_lt_one(moduli: &[u64]) -> Option<bool> {
    reciprocal_sum_fraction(moduli).map(|(numerator, denominator)| numerator < denominator)
}

fn checked_lcm(a: u64, b: u64) -> Result<u64, String> {
    let q = a / gcd(a, b);
    q.checked_mul(b)
        .ok_or_else(|| format!("LCM overflow while combining {a} and {b}"))
}

fn common_period(moduli: &[u64]) -> Result<u64, String> {
    let mut period = 1_u64;
    for &n in moduli {
        period = checked_lcm(period, n)?;
        if period > MAX_PERIOD {
            return Err(format!(
                "common period {period} exceeds checker limit {MAX_PERIOD}"
            ));
        }
    }
    Ok(period)
}

fn parse_candidate(text: &str) -> Result<Vec<Congruence>, String> {
    let mut words = text.split_whitespace();
    let k_word = words
        .next()
        .ok_or_else(|| "empty candidate: expected k".to_owned())?;
    let k: usize = k_word
        .parse()
        .map_err(|_| format!("invalid size k: {k_word:?}"))?;
    if k == 0 {
        return Err("k must be positive".to_owned());
    }
    if k > MAX_K {
        return Err(format!("k={k} exceeds checker limit {MAX_K}"));
    }

    let mut system: Vec<Congruence> = Vec::with_capacity(k);
    for i in 0..k {
        let n_word = words
            .next()
            .ok_or_else(|| format!("missing modulus on row {}", i + 1))?;
        let a_word = words
            .next()
            .ok_or_else(|| format!("missing residue on row {}", i + 1))?;
        let modulus: u64 = n_word
            .parse()
            .map_err(|_| format!("invalid modulus on row {}: {n_word:?}", i + 1))?;
        let residue: u64 = a_word
            .parse()
            .map_err(|_| format!("invalid residue on row {}: {a_word:?}", i + 1))?;
        if modulus <= 1 {
            return Err(format!("modulus on row {} must exceed 1", i + 1));
        }
        if residue >= modulus {
            return Err(format!(
                "residue {residue} on row {} is not canonical modulo {modulus}",
                i + 1
            ));
        }
        if let Some(previous) = system.last() {
            if modulus <= previous.modulus {
                return Err(format!(
                    "moduli must be strictly increasing: {} then {modulus}",
                    previous.modulus
                ));
            }
        }
        system.push(Congruence { modulus, residue });
    }
    if let Some(extra) = words.next() {
        return Err(format!("unexpected trailing token: {extra:?}"));
    }
    Ok(system)
}

fn verify_displayed_cover(system: &[Congruence], period: u64) -> Result<(), String> {
    for x in 0..period {
        if !system.iter().any(|c| x % c.modulus == c.residue) {
            return Err(format!(
                "the displayed congruences do not cover residue {x} modulo {period}"
            ));
        }
    }
    Ok(())
}

/// Add the residue class `residue (mod modulus)` to the coverage counters.
/// Returns the number of points that changed from uncovered to covered.
fn add_class(covered: &mut [u16], modulus: usize, residue: usize) -> usize {
    let mut newly_covered = 0;
    let mut x = residue;
    while x < covered.len() {
        if covered[x] == 0 {
            newly_covered += 1;
        }
        covered[x] = covered[x]
            .checked_add(1)
            .expect("coverage counter cannot overflow when k <= 128");
        x += modulus;
    }
    newly_covered
}

fn remove_class(covered: &mut [u16], modulus: usize, residue: usize) {
    let mut x = residue;
    while x < covered.len() {
        assert!(covered[x] > 0);
        covered[x] -= 1;
        x += modulus;
    }
}

struct DfsState<'a> {
    moduli: &'a [usize],
    period: usize,
    covered: Vec<u16>,
    assigned: Vec<bool>,
    residues: Vec<Option<u64>>,
    nodes: u64,
    max_nodes: u64,
}

impl DfsState<'_> {
    fn dfs(&mut self, uncovered: usize) -> Option<bool> {
        if self.nodes >= self.max_nodes {
            return None;
        }
        self.nodes += 1;
        if uncovered == 0 {
            return Some(true);
        }

        // Even with no overlaps, the remaining classes cannot cover more than
        // this many currently uncovered points.
        let capacity: usize = self
            .moduli
            .iter()
            .enumerate()
            .filter(|(i, _)| !self.assigned[*i])
            .map(|(_, &n)| self.period / n)
            .sum();
        if capacity < uncovered {
            return Some(false);
        }

        // Any extension must cover this point with one of the unassigned
        // moduli.  Branching on that responsible modulus is exhaustive and is
        // much smaller than branching over every possible residue separately.
        let x = self
            .covered
            .iter()
            .position(|&count| count == 0)
            .expect("uncovered count and counters disagree");

        let mut branches: Vec<(usize, usize)> = self
            .moduli
            .iter()
            .enumerate()
            .filter(|(i, _)| !self.assigned[*i])
            .map(|(i, &n)| {
                let residue = x % n;
                let gain = (residue..self.period)
                    .step_by(n)
                    .filter(|&y| self.covered[y] == 0)
                    .count();
                (i, gain)
            })
            .collect();
        branches.sort_unstable_by(|a, b| b.1.cmp(&a.1));

        for (i, _) in branches {
            let n = self.moduli[i];
            let residue = x % n;
            self.assigned[i] = true;
            self.residues[i] = Some(residue as u64);
            let gain = add_class(&mut self.covered, n, residue);
            match self.dfs(uncovered - gain) {
                Some(true) => return Some(true),
                None => {
                    remove_class(&mut self.covered, n, residue);
                    self.assigned[i] = false;
                    self.residues[i] = None;
                    return None;
                }
                Some(false) => {}
            }
            remove_class(&mut self.covered, n, residue);
            self.assigned[i] = false;
            self.residues[i] = None;
        }
        Some(false)
    }
}

/// Decide exactly whether one residue per listed modulus can cover all integers.
/// Translation invariance lets us fix the first listed residue to zero.
fn find_cover(moduli_u64: &[u64], max_nodes: u64) -> Result<(SearchOutcome, u64), String> {
    if moduli_u64.is_empty() {
        return Ok((SearchOutcome::NoCover, 0));
    }
    // The union bound is exact as a non-cover certificate and does not require
    // constructing a potentially enormous common period.
    if reciprocal_sum_lt_one(moduli_u64) == Some(true) {
        return Ok((SearchOutcome::NoCover, 0));
    }
    let period_u64 = common_period(moduli_u64)?;
    let period = usize::try_from(period_u64)
        .map_err(|_| format!("period {period_u64} does not fit usize"))?;
    let moduli: Vec<usize> = moduli_u64
        .iter()
        .map(|&n| usize::try_from(n).map_err(|_| format!("modulus {n} does not fit usize")))
        .collect::<Result<_, _>>()?;

    // The density/union bound is already a complete non-cover certificate in
    // many cases.
    let total_capacity: u128 = moduli_u64.iter().map(|&n| u128::from(period_u64 / n)).sum();
    if total_capacity < u128::from(period_u64) {
        return Ok((SearchOutcome::NoCover, 0));
    }

    let mut state = DfsState {
        moduli: &moduli,
        period,
        covered: vec![0; period],
        assigned: vec![false; moduli.len()],
        residues: vec![None; moduli.len()],
        nodes: 0,
        max_nodes,
    };

    // If a cover has first residue a, translating every class by -a gives a
    // cover whose first residue is zero.  Thus this symmetry break loses none.
    state.assigned[0] = true;
    state.residues[0] = Some(0);
    let first_gain = add_class(&mut state.covered, moduli[0], 0);
    let decision = state.dfs(period - first_gain);
    let outcome = match decision {
        Some(true) => SearchOutcome::Cover(state.residues.iter().map(|a| a.unwrap_or(0)).collect()),
        Some(false) => SearchOutcome::NoCover,
        None => SearchOutcome::BudgetExceeded,
    };
    Ok((outcome, state.nodes))
}

fn verify_irreducible(system: &[Congruence], max_nodes: u64) -> Result<u64, String> {
    let mut total_nodes = 0_u64;
    for deleted in 0..system.len() {
        let moduli: Vec<u64> = system
            .iter()
            .enumerate()
            .filter(|(i, _)| *i != deleted)
            .map(|(_, c)| c.modulus)
            .collect();
        let remaining_budget = max_nodes.saturating_sub(total_nodes);
        if remaining_budget == 0 {
            return Err(format!(
                "irreducibility search exhausted the global budget of {max_nodes} nodes"
            ));
        }
        let (outcome, nodes) = find_cover(&moduli, remaining_budget)?;
        total_nodes = total_nodes
            .checked_add(nodes)
            .ok_or_else(|| "search-node counter overflow".to_owned())?;
        match outcome {
            SearchOutcome::NoCover => {}
            SearchOutcome::BudgetExceeded => {
                return Err(format!(
                    "irreducibility undecided after {total_nodes} search nodes; no PASS"
                ));
            }
            SearchOutcome::Cover(residues) => {
                let witness = moduli
                    .iter()
                    .zip(residues.iter())
                    .map(|(n, a)| format!("{a} (mod {n})"))
                    .collect::<Vec<_>>()
                    .join(", ");
                return Err(format!(
                    "set is reducible: after deleting modulus {}, the remaining set covers via {witness}",
                    system[deleted].modulus
                ));
            }
        }
    }
    Ok(total_nodes)
}

/// Decide irreducibility of a bare modulus set, allowing a fresh residue
/// choice for every deletion.  Any budget exhaustion is an error, never a
/// negative answer.
fn decide_irreducible_moduli(moduli: &[u64], max_nodes: u64) -> Result<(bool, u64), String> {
    let (full, mut total_nodes) = find_cover(moduli, max_nodes)?;
    match full {
        SearchOutcome::NoCover => return Ok((false, total_nodes)),
        SearchOutcome::BudgetExceeded => {
            return Err(format!(
                "full-cover search for {moduli:?} exceeded {max_nodes} nodes"
            ));
        }
        SearchOutcome::Cover(_) => {}
    }
    for deleted in 0..moduli.len() {
        let sub: Vec<u64> = moduli
            .iter()
            .enumerate()
            .filter(|(i, _)| *i != deleted)
            .map(|(_, &n)| n)
            .collect();
        let (outcome, nodes) = find_cover(&sub, max_nodes)?;
        total_nodes = total_nodes
            .checked_add(nodes)
            .ok_or_else(|| "census node counter overflow".to_owned())?;
        match outcome {
            SearchOutcome::NoCover => {}
            SearchOutcome::Cover(_) => return Ok((false, total_nodes)),
            SearchOutcome::BudgetExceeded => {
                return Err(format!(
                    "deletion search for {moduli:?} exceeded {max_nodes} nodes"
                ));
            }
        }
    }
    Ok((true, total_nodes))
}

fn combinations_from(values: &[u64], k: usize) -> Vec<Vec<u64>> {
    fn rec(
        values: &[u64],
        next: usize,
        left: usize,
        current: &mut Vec<u64>,
        output: &mut Vec<Vec<u64>>,
    ) {
        if left == 0 {
            output.push(current.clone());
            return;
        }
        if values.len().saturating_sub(next) < left {
            return;
        }
        for i in next..=values.len() - left {
            current.push(values[i]);
            rec(values, i + 1, left - 1, current, output);
            current.pop();
        }
    }
    let mut output = Vec::new();
    rec(values, 0, k, &mut Vec::with_capacity(k), &mut output);
    output
}

fn candidates_via_common_multiples(k: usize, bound: u64) -> Vec<Vec<u64>> {
    let mut output = BTreeSet::new();
    for common in 1..=bound {
        let divisors: Vec<u64> = (2..=common).filter(|n| common % n == 0).collect();
        for set in combinations_from(&divisors, k) {
            output.insert(set);
        }
    }
    output.into_iter().collect()
}

fn mycielski_weight(mut n: u64) -> u64 {
    let mut weight = 0_u64;
    let mut p = 2_u64;
    while p <= n / p {
        let mut exponent = 0_u64;
        while n % p == 0 {
            n /= p;
            exponent += 1;
        }
        weight += exponent * (p - 1);
        p += 1;
    }
    if n > 1 {
        weight += n - 1;
    }
    weight
}

fn candidates_via_simpson(k: usize, bound: u64) -> Result<Vec<Vec<u64>>, String> {
    let mut output = Vec::new();
    for common in 1..=bound {
        // Every accepted set is enumerated once, at its actual LCM.  Simpson's
        // exact weight inequality rules out the other common multiples before
        // any divisor-subset materialization.
        if mycielski_weight(common) >= k as u64 {
            continue;
        }
        let divisors: Vec<u64> = (2..=common).filter(|n| common % n == 0).collect();
        for set in combinations_from(&divisors, k) {
            let actual_lcm = set.iter().try_fold(1_u64, |acc, &n| checked_lcm(acc, n))?;
            if actual_lcm == common {
                output.push(set);
            }
        }
    }
    output.sort();
    if output.windows(2).any(|w| w[0] == w[1]) {
        return Err("internal duplicate in exact-LCM candidate enumeration".to_owned());
    }
    Ok(output)
}

fn combinations(lo: u64, hi: u64, k: usize) -> Vec<Vec<u64>> {
    fn rec(next: u64, hi: u64, left: usize, current: &mut Vec<u64>, output: &mut Vec<Vec<u64>>) {
        if left == 0 {
            output.push(current.clone());
            return;
        }
        if next > hi || (hi - next + 1) < left as u64 {
            return;
        }
        let last = hi - left as u64 + 1;
        for n in next..=last {
            current.push(n);
            rec(n + 1, hi, left - 1, current, output);
            current.pop();
        }
    }
    let mut output = Vec::new();
    if k == 0 {
        output.push(Vec::new());
    } else if lo <= hi {
        rec(lo, hi, k, &mut Vec::with_capacity(k), &mut output);
    }
    output
}

/// Independently verify the complete census made finite by the formally
/// proved bound n_k <= 2^(k-1).  This checker recomputes every cover decision;
/// it shares no source with the generator under workspace/.
fn verify_census(max_k: usize, max_nodes: u64) -> Result<(), String> {
    if max_k > 6 {
        return Err("certified census mode currently supports only k<=6".to_owned());
    }
    let expected_k5 = vec![vec![2, 3, 4, 6, 12]];
    let expected_k6 = vec![
        vec![2, 3, 4, 6, 8, 24],
        vec![2, 3, 4, 8, 12, 24],
        vec![2, 3, 6, 8, 12, 24],
        vec![2, 4, 6, 8, 12, 24],
    ];
    let mut total_candidates = 0_usize;
    let mut total_nodes = 0_u64;
    for k in 1_usize..=max_k {
        let bound = 1_u64 << (k - 1);
        let candidates = combinations(2, bound, k);
        total_candidates += candidates.len();
        let mut answers = Vec::new();
        for moduli in candidates {
            let (is_irreducible, nodes) = decide_irreducible_moduli(&moduli, max_nodes)?;
            total_nodes = total_nodes
                .checked_add(nodes)
                .ok_or_else(|| "global census node counter overflow".to_owned())?;
            if is_irreducible {
                answers.push(moduli);
            }
        }
        let expected: &[Vec<u64>] = match k {
            5 => &expected_k5,
            6 => &expected_k6,
            _ => &[],
        };
        if answers != expected {
            return Err(format!(
                "small census mismatch at k={k}: computed {answers:?}, expected {expected:?}"
            ));
        }
        println!(
            "census k={k}: {} irreducible sets: {answers:?}",
            answers.len()
        );
        if k == 6 {
            let min_largest = answers.iter().filter_map(|s| s.last()).min().copied();
            let max_largest = answers.iter().filter_map(|s| s.last()).max().copied();
            let mut best = (0_u128, 1_u128);
            for set in &answers {
                let value = reciprocal_sum_fraction(set)
                    .ok_or_else(|| "reciprocal-sum overflow in k=6 answers".to_owned())?;
                let left = value
                    .0
                    .checked_mul(best.1)
                    .ok_or_else(|| "fraction comparison overflow".to_owned())?;
                let right = best
                    .0
                    .checked_mul(value.1)
                    .ok_or_else(|| "fraction comparison overflow".to_owned())?;
                if left > right {
                    best = value;
                }
            }
            if (min_largest, max_largest, best) != (Some(24), Some(24), (17, 12)) {
                return Err(format!(
                    "k=6 extrema mismatch: min={min_largest:?}, max={max_largest:?}, reciprocal={best:?}"
                ));
            }
            println!(
                "k=6 extrema: minimum largest modulus=24; maximum largest modulus=24; \
                 maximum reciprocal sum=17/12"
            );
        }
    }
    println!(
        "PASS: exhaustive census through k={max_k} checked {total_candidates} candidate sets \
         ({total_nodes} exhaustive DFS nodes)"
    );
    Ok(())
}

fn parse_expected_sets(text: &str, k: usize) -> Result<Vec<Vec<u64>>, String> {
    let mut sets = Vec::new();
    for (line_no, line) in text.lines().enumerate() {
        if line.trim().is_empty() {
            continue;
        }
        let set: Vec<u64> = line
            .split_whitespace()
            .map(|word| {
                word.parse().map_err(|_| {
                    format!(
                        "invalid expected-set token {word:?} on line {}",
                        line_no + 1
                    )
                })
            })
            .collect::<Result<_, _>>()?;
        if set.len() != k || set.windows(2).any(|w| w[0] >= w[1]) {
            return Err(format!(
                "expected-set line {} is not {k} strictly increasing moduli: {set:?}",
                line_no + 1
            ));
        }
        sets.push(set);
    }
    if sets.windows(2).any(|w| w[0] >= w[1]) {
        return Err("expected sets are not strictly lexicographically sorted".to_owned());
    }
    Ok(sets)
}

fn verify_k7_census(max_nodes: u64) -> Result<(), String> {
    let expected = vec![
        vec![2, 3, 4, 6, 8, 16, 48],
        vec![2, 3, 4, 6, 9, 18, 36],
        vec![2, 3, 4, 6, 16, 24, 48],
        vec![2, 3, 4, 8, 12, 16, 48],
        vec![2, 3, 4, 8, 16, 24, 48],
        vec![2, 3, 4, 9, 12, 18, 36],
        vec![2, 3, 4, 12, 16, 24, 48],
        vec![2, 3, 6, 8, 12, 16, 48],
        vec![2, 3, 6, 9, 12, 18, 36],
        vec![2, 3, 6, 12, 16, 24, 48],
        vec![2, 4, 6, 8, 12, 16, 48],
        vec![2, 4, 6, 8, 16, 24, 48],
        vec![2, 4, 6, 9, 12, 18, 36],
        vec![2, 4, 6, 12, 16, 24, 48],
        vec![2, 4, 8, 12, 16, 24, 48],
    ];
    let candidates = candidates_via_common_multiples(7, 64);
    let candidate_count = candidates.len();
    let mut answers = Vec::new();
    let mut total_nodes = 0_u64;
    for moduli in candidates {
        let (is_irreducible, nodes) = decide_irreducible_moduli(&moduli, max_nodes)?;
        total_nodes = total_nodes
            .checked_add(nodes)
            .ok_or_else(|| "k=7 census node counter overflow".to_owned())?;
        if is_irreducible {
            answers.push(moduli);
        }
    }
    if answers != expected {
        return Err(format!(
            "k=7 census mismatch: computed {answers:?}, expected {expected:?}"
        ));
    }
    let min_largest = answers.iter().filter_map(|s| s.last()).min().copied();
    let max_largest = answers.iter().filter_map(|s| s.last()).max().copied();
    let mut best = (0_u128, 1_u128);
    for set in &answers {
        let value = reciprocal_sum_fraction(set)
            .ok_or_else(|| "reciprocal-sum overflow in k=7 answers".to_owned())?;
        let left = value
            .0
            .checked_mul(best.1)
            .ok_or_else(|| "k=7 fraction comparison overflow".to_owned())?;
        let right = best
            .0
            .checked_mul(value.1)
            .ok_or_else(|| "k=7 fraction comparison overflow".to_owned())?;
        if left > right {
            best = value;
        }
    }
    if (candidate_count, min_largest, max_largest, best) != (378, Some(36), Some(48), (35, 24)) {
        return Err(format!(
            "k=7 aggregate mismatch: candidates={candidate_count}, min={min_largest:?}, \
             max={max_largest:?}, reciprocal={best:?}"
        ));
    }
    println!(
        "census k=7: {} irreducible sets: {answers:?}",
        answers.len()
    );
    println!(
        "k=7 extrema: minimum largest modulus=36; maximum largest modulus=48; \
         maximum reciprocal sum=35/24"
    );
    println!(
        "PASS: exhaustive k=7 common-multiple census checked {candidate_count} candidate sets \
         ({total_nodes} exhaustive DFS nodes)"
    );
    Ok(())
}

fn verify_k8_census(max_nodes: u64) -> Result<(), String> {
    let expected = parse_expected_sets(include_str!("../fixtures/k8_expected.txt"), 8)?;
    if expected.len() != 65 {
        return Err(format!(
            "expected k=8 fixture has {} sets, not 65",
            expected.len()
        ));
    }
    let candidates = candidates_via_common_multiples(8, 128);
    let candidate_count = candidates.len();
    let mut answers = Vec::new();
    let mut total_nodes = 0_u64;
    for moduli in candidates {
        let (is_irreducible, nodes) = decide_irreducible_moduli(&moduli, max_nodes)?;
        total_nodes = total_nodes
            .checked_add(nodes)
            .ok_or_else(|| "k=8 census node counter overflow".to_owned())?;
        if is_irreducible {
            answers.push(moduli);
        }
    }
    if answers != expected {
        return Err(format!(
            "k=8 census mismatch: computed {} sets, expected {}",
            answers.len(),
            expected.len()
        ));
    }
    let min_largest = answers.iter().filter_map(|s| s.last()).min().copied();
    let max_largest = answers.iter().filter_map(|s| s.last()).max().copied();
    let mut best = (0_u128, 1_u128);
    for set in &answers {
        let value = reciprocal_sum_fraction(set)
            .ok_or_else(|| "reciprocal-sum overflow in k=8 answers".to_owned())?;
        let left = value
            .0
            .checked_mul(best.1)
            .ok_or_else(|| "k=8 fraction comparison overflow".to_owned())?;
        let right = best
            .0
            .checked_mul(value.1)
            .ok_or_else(|| "k=8 fraction comparison overflow".to_owned())?;
        if left > right {
            best = value;
        }
    }
    if (candidate_count, min_largest, max_largest, best) != (7443, Some(36), Some(96), (14, 9)) {
        return Err(format!(
            "k=8 aggregate mismatch: candidates={candidate_count}, min={min_largest:?}, \
             max={max_largest:?}, reciprocal={best:?}"
        ));
    }
    println!(
        "census k=8: {} irreducible sets (exact list in fixtures/k8_expected.txt)",
        answers.len()
    );
    println!(
        "k=8 extrema: minimum largest modulus=36; maximum largest modulus=96; \
         maximum reciprocal sum=14/9"
    );
    println!(
        "PASS: exhaustive k=8 common-multiple census checked {candidate_count} candidate sets \
         ({total_nodes} exhaustive DFS nodes)"
    );
    Ok(())
}

fn verify_k9_census(max_nodes: u64) -> Result<(), String> {
    let expected = parse_expected_sets(include_str!("../fixtures/k9_expected.txt"), 9)?;
    if expected.len() != 318 {
        return Err(format!(
            "expected k=9 fixture has {} sets, not 318",
            expected.len()
        ));
    }
    let candidates = candidates_via_common_multiples(9, 256);
    let candidate_count = candidates.len();
    let decisions: Vec<Result<(bool, u64), String>> = candidates
        .par_iter()
        .map(|moduli| decide_irreducible_moduli(moduli, max_nodes))
        .collect();
    let mut answers = Vec::new();
    let mut total_nodes = 0_u64;
    for (moduli, decision) in candidates.into_iter().zip(decisions) {
        let (is_irreducible, nodes) = decision?;
        total_nodes = total_nodes
            .checked_add(nodes)
            .ok_or_else(|| "k=9 census node counter overflow".to_owned())?;
        if is_irreducible {
            answers.push(moduli);
        }
    }
    if answers != expected {
        return Err(format!(
            "k=9 census mismatch: computed {} sets, expected {}",
            answers.len(),
            expected.len()
        ));
    }
    let min_largest = answers.iter().filter_map(|s| s.last()).min().copied();
    let max_largest = answers.iter().filter_map(|s| s.last()).max().copied();
    let mut best = (0_u128, 1_u128);
    for set in &answers {
        let value = reciprocal_sum_fraction(set)
            .ok_or_else(|| "reciprocal-sum overflow in k=9 answers".to_owned())?;
        let left = value
            .0
            .checked_mul(best.1)
            .ok_or_else(|| "k=9 fraction comparison overflow".to_owned())?;
        let right = best
            .0
            .checked_mul(value.1)
            .ok_or_else(|| "k=9 fraction comparison overflow".to_owned())?;
        if left > right {
            best = value;
        }
    }
    if (candidate_count, min_largest, max_largest, best) != (159_225, Some(30), Some(192), (17, 10))
    {
        return Err(format!(
            "k=9 aggregate mismatch: candidates={candidate_count}, min={min_largest:?}, \
             max={max_largest:?}, reciprocal={best:?}"
        ));
    }
    println!(
        "census k=9: {} irreducible sets (exact list in fixtures/k9_expected.txt)",
        answers.len()
    );
    println!(
        "k=9 extrema: minimum largest modulus=30; maximum largest modulus=192; \
         maximum reciprocal sum=17/10"
    );
    println!(
        "PASS: exhaustive k=9 common-multiple census checked {candidate_count} candidate sets \
         ({total_nodes} exhaustive DFS nodes)"
    );
    Ok(())
}

fn verify_k10_census(max_nodes: u64) -> Result<(), String> {
    let expected = parse_expected_sets(include_str!("../fixtures/k10_expected.txt"), 10)?;
    if expected.len() != 2_102 {
        return Err(format!(
            "expected k=10 fixture has {} sets, not 2102",
            expected.len()
        ));
    }
    let candidates = candidates_via_simpson(10, 512)?;
    let candidate_count = candidates.len();
    let decisions: Vec<Result<(bool, u64), String>> = candidates
        .par_iter()
        .map(|moduli| decide_irreducible_moduli(moduli, max_nodes))
        .collect();
    let mut answers = Vec::new();
    let mut total_nodes = 0_u64;
    for (moduli, decision) in candidates.into_iter().zip(decisions) {
        let (is_irreducible, nodes) = decision?;
        total_nodes = total_nodes
            .checked_add(nodes)
            .ok_or_else(|| "k=10 census node counter overflow".to_owned())?;
        if is_irreducible {
            answers.push(moduli);
        }
    }
    if answers != expected {
        return Err(format!(
            "k=10 census mismatch: computed {} sets, expected {}",
            answers.len(),
            expected.len()
        ));
    }
    let min_largest = answers.iter().filter_map(|s| s.last()).min().copied();
    let max_largest = answers.iter().filter_map(|s| s.last()).max().copied();
    let mut best = (0_u128, 1_u128);
    for set in &answers {
        let value = reciprocal_sum_fraction(set)
            .ok_or_else(|| "reciprocal-sum overflow in k=10 answers".to_owned())?;
        let left = value
            .0
            .checked_mul(best.1)
            .ok_or_else(|| "k=10 fraction comparison overflow".to_owned())?;
        let right = best
            .0
            .checked_mul(value.1)
            .ok_or_else(|| "k=10 fraction comparison overflow".to_owned())?;
        if left > right {
            best = value;
        }
    }
    if (candidate_count, min_largest, max_largest, best) != (28_457, Some(30), Some(384), (109, 60))
    {
        return Err(format!(
            "k=10 aggregate mismatch: candidates={candidate_count}, min={min_largest:?}, \
             max={max_largest:?}, reciprocal={best:?}"
        ));
    }
    println!(
        "census k=10: {} irreducible sets (exact list in fixtures/k10_expected.txt)",
        answers.len()
    );
    println!(
        "k=10 extrema: minimum largest modulus=30; maximum largest modulus=384; \
         maximum reciprocal sum=109/60"
    );
    println!(
        "PASS: exhaustive k=10 Simpson-filtered common-multiple census checked \
         {candidate_count} candidate sets ({total_nodes} exhaustive DFS nodes)"
    );
    Ok(())
}

fn verify_k11_census(max_nodes: u64) -> Result<(), String> {
    let expected = parse_expected_sets(include_str!("../fixtures/k11_expected.txt"), 11)?;
    if expected.len() != 17_040 {
        return Err(format!(
            "expected k=11 fixture has {} sets, not 17040",
            expected.len()
        ));
    }
    let candidates = candidates_via_simpson(11, 1024)?;
    let candidate_count = candidates.len();
    let decisions: Vec<Result<(bool, u64), String>> = candidates
        .par_iter()
        .map(|moduli| decide_irreducible_moduli(moduli, max_nodes))
        .collect();
    let mut answers = Vec::new();
    let mut total_nodes = 0_u64;
    for (moduli, decision) in candidates.into_iter().zip(decisions) {
        let (is_irreducible, nodes) = decision?;
        total_nodes = total_nodes
            .checked_add(nodes)
            .ok_or_else(|| "k=11 census node counter overflow".to_owned())?;
        if is_irreducible {
            answers.push(moduli);
        }
    }
    if answers != expected {
        return Err(format!(
            "k=11 census mismatch: computed {} sets, expected {}",
            answers.len(),
            expected.len()
        ));
    }
    let min_largest = answers.iter().filter_map(|s| s.last()).min().copied();
    let max_largest = answers.iter().filter_map(|s| s.last()).max().copied();
    let mut best = (0_u128, 1_u128);
    for set in &answers {
        let value = reciprocal_sum_fraction(set)
            .ok_or_else(|| "reciprocal-sum overflow in k=11 answers".to_owned())?;
        let left = value
            .0
            .checked_mul(best.1)
            .ok_or_else(|| "k=11 fraction comparison overflow".to_owned())?;
        let right = best
            .0
            .checked_mul(value.1)
            .ok_or_else(|| "k=11 fraction comparison overflow".to_owned())?;
        if left > right {
            best = value;
        }
    }
    if (candidate_count, min_largest, max_largest, best) != (343_875, Some(36), Some(768), (28, 15))
    {
        return Err(format!(
            "k=11 aggregate mismatch: candidates={candidate_count}, min={min_largest:?}, \
             max={max_largest:?}, reciprocal={best:?}"
        ));
    }
    println!(
        "census k=11: {} irreducible sets (exact list in fixtures/k11_expected.txt)",
        answers.len()
    );
    println!(
        "k=11 extrema: minimum largest modulus=36; maximum largest modulus=768; \
         maximum reciprocal sum=28/15"
    );
    println!(
        "PASS: exhaustive k=11 Simpson-filtered exact-LCM census checked \
         {candidate_count} candidate sets ({total_nodes} exhaustive DFS nodes)"
    );
    Ok(())
}

fn usage(program: &str) -> String {
    format!(
        "usage: {program} <candidate.txt> [--max-nodes N]\n\
         or: {program} --verify-small-census [--max-nodes N]\n\
         or: {program} --verify-k6-census [--max-nodes N]\n\
         or: {program} --verify-k7-census [--max-nodes N]\n\
         or: {program} --verify-k8-census [--max-nodes N]\n\
         or: {program} --verify-k9-census [--max-nodes N]\n\
         or: {program} --verify-k10-census [--max-nodes N]\n\
         or: {program} --verify-k11-census [--max-nodes N]\n\
         candidate format: first token k, followed by exactly k pairs 'modulus residue'"
    )
}

fn run() -> Result<(), String> {
    let args: Vec<String> = env::args().collect();
    if args.len() != 2 && args.len() != 4 {
        return Err(usage(&args[0]));
    }
    let mut max_nodes = DEFAULT_MAX_SEARCH_NODES;
    if args.len() == 4 {
        if args[2] != "--max-nodes" {
            return Err(usage(&args[0]));
        }
        max_nodes = args[3]
            .parse()
            .map_err(|_| format!("invalid search-node budget: {:?}", args[3]))?;
        if max_nodes == 0 {
            return Err("search-node budget must be positive".to_owned());
        }
    }
    if args[1] == "--verify-small-census" {
        return verify_census(5, max_nodes);
    }
    if args[1] == "--verify-k6-census" {
        return verify_census(6, max_nodes);
    }
    if args[1] == "--verify-k7-census" {
        return verify_k7_census(max_nodes);
    }
    if args[1] == "--verify-k8-census" {
        return verify_k8_census(max_nodes);
    }
    if args[1] == "--verify-k9-census" {
        return verify_k9_census(max_nodes);
    }
    if args[1] == "--verify-k10-census" {
        return verify_k10_census(max_nodes);
    }
    if args[1] == "--verify-k11-census" {
        return verify_k11_census(max_nodes);
    }

    let text = fs::read_to_string(&args[1])
        .map_err(|e| format!("cannot read candidate {:?}: {e}", args[1]))?;
    let system = parse_candidate(&text)?;
    let moduli: Vec<u64> = system.iter().map(|c| c.modulus).collect();
    let period = common_period(&moduli)?;
    verify_displayed_cover(&system, period)?;
    let nodes = verify_irreducible(&system, max_nodes)?;
    println!(
        "PASS: {} distinct moduli; displayed classes cover period {}; every one-modulus deletion is non-covering ({} exhaustive DFS nodes)",
        system.len(), period, nodes
    );
    Ok(())
}

fn main() {
    if let Err(message) = run() {
        eprintln!("FAIL: {message}");
        process::exit(1);
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn parses_strictly() {
        assert!(parse_candidate("2 2 0 3 1").is_ok());
        assert!(parse_candidate("2 2 0 2 1").is_err());
        assert!(parse_candidate("1 2 2").is_err());
        assert!(parse_candidate("1 2 0 trailing").is_err());
    }

    #[test]
    fn classic_five_classes_cover() {
        let s = parse_candidate("5 2 0 3 0 4 1 6 5 12 7").unwrap();
        verify_displayed_cover(&s, 12).unwrap();
    }

    #[test]
    fn four_moduli_do_not_cover() {
        let (outcome, _) = find_cover(&[2, 3, 4, 6], 1_000_000).unwrap();
        assert!(matches!(outcome, SearchOutcome::NoCover));
    }

    #[test]
    fn combination_counts_are_exact() {
        assert_eq!(combinations(2, 4, 3), vec![vec![2, 3, 4]]);
        assert_eq!(combinations(2, 8, 4).len(), 35);
        assert_eq!(combinations(2, 32, 6).len(), 736_281);
        assert_eq!(candidates_via_common_multiples(7, 64).len(), 378);
        assert_eq!(candidates_via_common_multiples(8, 128).len(), 7_443);
        assert_eq!(candidates_via_common_multiples(9, 256).len(), 159_225);
        assert_eq!(candidates_via_simpson(10, 512).unwrap().len(), 28_457);
        assert!(combinations(2, 1, 1).is_empty());
    }

    #[test]
    fn expected_fixture_parses() {
        let sets = parse_expected_sets(include_str!("../fixtures/k8_expected.txt"), 8).unwrap();
        assert_eq!(sets.len(), 65);
        let sets = parse_expected_sets(include_str!("../fixtures/k9_expected.txt"), 9).unwrap();
        assert_eq!(sets.len(), 318);
        let sets = parse_expected_sets(include_str!("../fixtures/k10_expected.txt"), 10).unwrap();
        assert_eq!(sets.len(), 2_102);
    }

    #[test]
    fn exact_density_pruning() {
        assert_eq!(reciprocal_sum_lt_one(&[3, 4, 5]), Some(true));
        assert_eq!(reciprocal_sum_lt_one(&[2, 3, 6]), Some(false));
        assert_eq!(reciprocal_sum_lt_one(&[2, 3, 7]), Some(true));
        assert_eq!(
            reciprocal_sum_fraction(&[2, 3, 4, 6, 8, 24]),
            Some((17, 12))
        );
    }
}
