ocas-calc 0.30.0

//! Machine-checkable certificates for integration results (0.28.0).
//!
//! A certificate is a *symbolic* proof that an antiderivative is correct:
//! the difference `D(F) − f` must reduce to zero in a normal form that the
//! engine can evaluate exactly. It replaces "the string contains no
//! `Integral(`" and "sampled numerics looked fine" as the acceptance
//! criterion for an emitted answer; the numeric oracle in
//! `ocas-tests/src/integral_eval.rs` is demoted to a probing and
//! regression tool.
//!
//! Three layers are tried, cheapest first:
//!
//! 1. **Structural** — `normalize(simplify(D(F) − f)) == 0`. Covers the
//!    common case where the engine inverted its own derivative rules.
//! 2. **Elementary field** — the difference is rewritten to
//!    exponent/logarithm form (trigonometric heads through
//!    [`super::trig::trig_to_exp`]) and embedded into the multivariate
//!    rational function field generated by `x`, every residual function
//!    atom, every constant symbol, and every non-integer power left in the
//!    difference. Zero in that field is an exact decision procedure.
//! 3. **Radical field** (0.30.0) — fractional powers `S^(p/q)` and `sqrt`
//!    heads are rewritten to fresh generators `r` carrying the exact
//!    relation `r^q = S` (the relation set is triangular, so reduction to
//!    a normal form is iterated polynomial remainder, no full Gröbner).
//!    This is what certifies answers carrying `sqrt(2)`-style literals
//!    (`sqrt(2)·sqrt(2) = 2` is invisible to the elementary layer) and
//!    radicals of non-constant bases such as `1/sqrt(1−x²)`.
//!
//! # Soundness
//!
//! Layers 2–3 list extra atoms as *independent* generators. Treating a
//! constant such as `sqrt(2)` as a free variable can only produce false
//! negatives (an identity that needs `sqrt(2)·sqrt(2) = 2` is not seen),
//! never false positives: a polynomial identity in free variables stays
//! valid under every substitution. The layers therefore never certify a
//! wrong answer, but they can decline to certify a correct one.
//!
//! Two honest limitations, both documented for callers:
//!
//! - The certificate is computed with this crate's own derivative table
//!   ([`crate::diff`]), so it is relative to that table. The independent
//!   numeric oracle is kept as the cross-check (the CI suite requires
//!   "oracle says mismatch ⇒ no certificate").
//! - Real-valued branches and domains are *not* part of the certificate:
//!   a form-valid antiderivative can still be wrong on a restricted
//!   interval (the half-power sheet factors are the standing example).
//!
//! Diagnostic switch: `OCAS_CERT_TIMING=1` prints the difference and field
//! sizes of every attempt (the same convention as `OCAS_INTEGRATE_TRACE`); it
//! is a debugging aid, not an API.

use ocas_atom::normalize::normalize;
use ocas_atom::walk::collect_funs;
use ocas_atom::{Atom, AtomArena, AtomNode, Symbol};
use ocas_domain::{Domain, RationalDomain};
use ocas_poly::{Lex, SparseMultivariatePolynomial};
use ocas_rewrite::rules::default_rules;
use ocas_rewrite::simplify::simplify;

use super::trig::trig_to_exp;
use super::{contains_integral, fraction_exponent, node_count};
use crate::derivative::diff;
use crate::rules::calculus_rules;
use crate::tower::build::{GenInfo, GenKind};
use crate::tower::convert::atom_to_rational_extended;
use crate::tower::elem::KElem;
use crate::tower::merge;

/// Sparse multivariate polynomial over `ℚ` (the radical layer's normal
/// form carrier; the same representation the tower elements use).
type Sparse = SparseMultivariatePolynomial<RationalDomain, Lex>;

/// Node-count ceiling for the difference `D(F) − f` before any layer runs.
const MAX_CERT_NODES: usize = 20_000;

/// Maximum number of generators the field layers may adjoin.
const MAX_CERT_GENS: usize = 16;

/// Upper bound on the number of additive terms the difference expands to.
///
/// `collect_terms` and the field conversion both expand products over sums,
/// and the expanded form is what blows up on Weierstrass-style answers
/// (measured: a 1 000-node difference whose expansion has `8³`-sized
/// factors took 15 s in the field layer). The estimate is computed from the
/// difference itself — no expansion, no polynomial arithmetic — so it is a
/// deterministic gate rather than a wall-clock one.
const MAX_CERT_EXPANSION: u64 = 4_096;

/// Largest integer exponent whose power may enter the field layer.
const MAX_CERT_EXPONENT: u64 = 64;

/// Work ceiling for the field layer, as `expansion × generators`.
///
/// The field arithmetic normalizes rational functions (multivariate
/// gcd/division over `ℚ`), whose cost grows with both the expanded size and
/// the number of generators: a 2 400-node difference over 10 generators
/// measured 15 s, while the same shape over 5 generators costs ~20 ms. The
/// product is a deterministic proxy for that cost.
///
/// 0.29.0 (Wave G) re-measured the frontier: 400 → 1200 lifted
/// `certified_rate` 29.1% → 31.7% but pushed 3 previously-solved cases over
/// the 10 s per-case budget (timeouts 12 → 15) and cost +14% wall clock —
/// the same coverage-for-time trade 0.28.0 recorded. The genuine fix is the
/// 0.32.0 complexity work (modular certification), not a looser budget.
const MAX_CERT_FIELD_WORK: u64 = 400;

/// Which layer produced a certificate.
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub enum CertMethod {
    /// `normalize(simplify(D(F) − f))` was exactly zero.
    Structural,
    /// `D(F) − f` vanished in the rational function field generated by the
    /// atoms of the expression.
    ElementaryField,
    /// `D(F) − f` vanished after adjoining radical generators (reserved for
    /// the 0.28.0 W3 layer).
    RadicalField,
}

impl CertMethod {
    /// Stable lowercase name, used by the corpus harness report.
    ///
    /// # Example
    ///
    /// ```
    /// use ocas_calc::integral::certify::CertMethod;
    /// assert_eq!(CertMethod::Structural.as_str(), "structural");
    /// ```
    pub fn as_str(self) -> &'static str {
        match self {
            CertMethod::Structural => "structural",
            CertMethod::ElementaryField => "field",
            CertMethod::RadicalField => "radical",
        }
    }
}

/// A machine-checkable proof that `F` is an antiderivative of `f`.
///
/// The witness kept in the certificate is the reduced difference, which is
/// always the numeric constant zero for a successful certificate.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct Certificate<'a> {
    /// Layer that established the certificate.
    pub method: CertMethod,
    /// The reduced normal form of `D(F) − f` (zero).
    pub difference: Atom<'a>,
}

/// Why [`certify`] declined to certify a result.
///
/// `NonZero` is a **triage hint**, not proof of an error: it means the
/// difference was nonzero in the free-generator normal form, which can also
/// happen when a genuine identity needs a relation the layer does not model
/// (e.g. `I² = −1`, `exp(u)·exp(−u) = 1`, `sqrt(2)·sqrt(2) = 2`). Only the
/// numeric oracle (and human review) may turn it into a wrong-answer claim.
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub enum CertDecline {
    /// The difference could not be embedded into the implemented field
    /// (unknown relation, deeper composite, unmodelled head).
    NotInField,
    /// A deterministic budget (node count, generator count) was exhausted.
    Budget,
    /// The difference had a nonzero normal form in the free-generator field
    /// (triage hint — see the enum docs).
    NonZero,
}

impl CertDecline {
    /// Stable lowercase name, used by the corpus harness report.
    ///
    /// # Example
    ///
    /// ```
    /// use ocas_calc::integral::certify::CertDecline;
    /// assert_eq!(CertDecline::NotInField.as_str(), "notinfield");
    /// ```
    pub fn as_str(self) -> &'static str {
        match self {
            CertDecline::NotInField => "notinfield",
            CertDecline::Budget => "budget",
            CertDecline::NonZero => "nonzero",
        }
    }
}

/// Certify that `antiderivative` differentiates back to `integrand`.
///
/// `integrand` must be the *normalized* integrand that was handed to the
/// integration pipeline, so that both sides live in the same atom
/// representation.
///
/// # Example
///
/// ```
/// use ocas_atom::{AtomArena, Symbol};
/// use ocas_calc::integral::certify::{CertMethod, certify};
/// use ocas_core::arena::Arena;
///
/// let arena = Arena::new();
/// let ctx = AtomArena::new(&arena);
/// let x = ctx.var("x");
/// // ∫ 2·x dx = x²
/// let big_f = ctx.pow(x, ctx.num(2));
/// let integrand = ctx.mul(&[ctx.num(2), x]);
/// let cert = certify(&ctx, integrand, big_f, Symbol::new("x")).expect("certificate");
/// assert_eq!(cert.method, CertMethod::Structural);
/// ```
pub fn certify<'a>(
    ctx: &'a AtomArena<'a>,
    integrand: Atom<'a>,
    antiderivative: Atom<'a>,
    var: Symbol,
) -> Result<Certificate<'a>, CertDecline> {
    certify_with_limits(
        ctx,
        integrand,
        antiderivative,
        var,
        MAX_CERT_NODES,
        MAX_CERT_GENS,
    )
}

/// [`certify`] with explicit layer limits; used by the unit tests to exercise
/// the budget paths without building a 20 000-node expression.
pub(crate) fn certify_with_limits<'a>(
    ctx: &'a AtomArena<'a>,
    integrand: Atom<'a>,
    antiderivative: Atom<'a>,
    var: Symbol,
    max_nodes: usize,
    max_gens: usize,
) -> Result<Certificate<'a>, CertDecline> {
    // A result that still carries `Integral(…)` is not an antiderivative.
    if contains_integral(antiderivative) {
        return Err(CertDecline::NotInField);
    }

    let derived = diff(ctx, antiderivative, var);
    // `D(F) − f`, built as `D(F) + (−1)·f` so the arena can hash-cons it.
    let minus_one = ctx.num(-1);
    let reduced = ctx.add(&[derived, ctx.mul(&[minus_one, integrand])]);
    let reduced = normalize(ctx, reduced);
    let default_rules = default_rules(ctx, &crate::pattern_alloc::VecAlloc);
    let calc_rules = calculus_rules(ctx, &crate::pattern_alloc::VecAlloc);
    let reduced = simplify(ctx, reduced, &default_rules, 20);
    let reduced = simplify(ctx, reduced, &calc_rules, 10);
    let reduced = normalize(ctx, reduced);
    if std::env::var_os("OCAS_CERT_TIMING").is_some() {
        eprintln!(
            "[cert] difference nodes={} (integrand {} / antiderivative {})",
            node_count(reduced),
            node_count(integrand),
            node_count(antiderivative)
        );
    }

    if node_count(reduced) > max_nodes {
        return Err(CertDecline::Budget);
    }

    // Layer 1 — structural zero. `collect_terms` supplies the like-term
    // collection (and bounded expansion) the general simplifier lacks, so
    // `2·x − 2·x`, `x − 2·(2⁻¹)·x` and their polynomial relatives collapse.
    let zero = ctx.num(0);
    if is_zero(reduced) {
        return Ok(Certificate {
            method: CertMethod::Structural,
            difference: zero,
        });
    }
    // Everything below expands; skip the expanding layers when the estimate
    // overflows, but still give the radical layer its chance (fractional
    // powers are exactly what the estimate cannot size).
    let expansion = expansion_size(reduced);
    if expansion <= MAX_CERT_EXPANSION {
        let collected = crate::ode::util::collect_terms(ctx, reduced);
        if node_count(collected) > max_nodes {
            return Err(CertDecline::Budget);
        }
        if is_zero(collected) {
            return Ok(Certificate {
                method: CertMethod::Structural,
                difference: zero,
            });
        }
    }

    // Layer 2 — elementary function field over the atoms of the difference.
    // Trigonometric heads are rewritten to exponentials first: they carry
    // identities (`sin² + cos² = 1`) the free-generator field cannot see.
    //
    // The *uncollected* difference is used from here on: `collect_terms`
    // expands products over sums, and the expanded form makes the field
    // arithmetic orders of magnitude more expensive (measured: 33 s versus
    // 22 ms on `(d + e·x)⁸/(a² + 2abx + b²x²)³`).
    let rewritten = trig_to_exp(ctx, reduced).unwrap_or(reduced);
    let rewritten = reduce_imaginary_unit(ctx, rewritten);
    let rewritten = merge_dependent_atoms(ctx, rewritten, var);
    let rewritten = normalize(ctx, rewritten);

    let mut gens: Vec<Atom<'a>> = Vec::with_capacity(8);
    gens.push(ctx.var(var.as_str()));
    let mut overflow = false;
    let complete = collect_generators(rewritten, &mut gens, max_gens, &mut overflow);
    if overflow {
        return Err(CertDecline::Budget);
    }

    if std::env::var_os("OCAS_CERT_TIMING").is_some() {
        eprintln!(
            "[cert] field input nodes={} expansion={} gens={}",
            node_count(rewritten),
            expansion,
            gens.len()
        );
    }
    let layer2 = if !complete {
        Err(CertDecline::NotInField)
    } else if expansion.saturating_mul(gens.len() as u64) > MAX_CERT_FIELD_WORK {
        Err(CertDecline::Budget)
    } else {
        match atom_to_rational_extended(rewritten, &gens, gens.len()) {
            Some(rf) if rf.numerator.n_terms() == 0 => Ok(Certificate {
                method: CertMethod::ElementaryField,
                difference: zero,
            }),
            Some(_) => Err(CertDecline::NonZero),
            None => Err(CertDecline::NotInField),
        }
    };
    match layer2 {
        Ok(_) => layer2,
        // A free-generator normal form can miss exact radical relations
        // (`sqrt(2)·sqrt(2) = 2`, `r² = 1−x²`), and the expansion estimate
        // declines fractional powers outright — both failure modes fall
        // through to the radical layer, which has its own caps. When the
        // radical layer cannot decide either, the elementary layer's
        // verdict stands.
        Err(e) => radical_layer(ctx, rewritten, var, max_gens).or(Err(e)),
    }
}

// ---------------------------------------------------------------------
//  Layer 3 — radical field (0.30.0)
// ---------------------------------------------------------------------

/// Cap on the denominator `q` of a fractional exponent `p/q`.
const MAX_CERT_RADICAL_Q: i64 = 8;
/// Cap on `|p|` for a fractional exponent `p/q`.
const MAX_CERT_RADICAL_P: i64 = 16;
/// Cap on the number of radical generators adjoined by the radical layer
/// (in addition to [`MAX_CERT_GENS`] on the whole generator set).
const MAX_CERT_RADICALS: usize = 8;
/// Term cap for one reduction step (sparse numerator/denominator size while
/// reducing modulo the radical relations).
const MAX_CERT_REDUCTION_TERMS: usize = 8192;
/// Node cap on the rewritten difference entering the radical layer (the
/// expansion estimate cannot see fractional powers, so the layer gates on
/// the tree size directly).
const MAX_CERT_RADICAL_NODES: usize = 2_000;
/// Cap on the numerator/denominator size (terms) entering the triangular
/// reduction — the entry-side counterpart of [`MAX_CERT_REDUCTION_TERMS`]
/// (0.30.0 W4a, measured on a 2816-term denominator case).
const MAX_CERT_RADICAL_POLY_TERMS: usize = 2_048;

/// One radical relation `r^q = base` (`base` already rewritten, so it only
/// mentions generators created *earlier* — the relation set is triangular).
struct RadicalRelation {
    /// Index of the generator variable `r` in the field's generator list.
    var_index: usize,
    /// The exponent `q` of the relation.
    q: u32,
    /// The relation's right-hand side as a sparse polynomial (denominator
    /// free; a rational right-hand side is declined).
    base: Sparse,
}

/// Rewrite every radical site of `expr` to a power of a fresh generator and
/// collect the `(base, q)` groups. Returns `None` when a site exceeds the
/// deterministic caps.
fn rewrite_radicals<'a>(
    ctx: &'a AtomArena<'a>,
    expr: Atom<'a>,
    groups: &mut Vec<(Atom<'a>, i64, Atom<'a>)>,
) -> Option<Atom<'a>> {
    match expr.node() {
        AtomNode::Num(_) | AtomNode::Var(_) => Some(expr),
        AtomNode::Add(args) => {
            let rebuilt: Vec<Atom> = args
                .iter()
                .map(|a| rewrite_radicals(ctx, *a, groups))
                .collect::<Option<_>>()?;
            Some(ctx.add(&rebuilt))
        }
        AtomNode::Mul(args) => {
            let rebuilt: Vec<Atom> = args
                .iter()
                .map(|a| rewrite_radicals(ctx, *a, groups))
                .collect::<Option<_>>()?;
            Some(ctx.mul(&rebuilt))
        }
        AtomNode::Fun(name, args) => {
            if name.as_str() == "sqrt" && args.len() == 1 {
                let base = rewrite_radicals(ctx, args[0], groups)?;
                let (r, _) = radical_group(ctx, base, 2, groups)?;
                return Some(r);
            }
            let rebuilt: Vec<Atom> = args
                .iter()
                .map(|a| rewrite_radicals(ctx, *a, groups))
                .collect::<Option<_>>()?;
            Some(ctx.fun(name.as_str(), &rebuilt))
        }
        AtomNode::Pow(base, exp) => {
            if let Some((p, q)) = fraction_exponent(*exp) {
                if !(2..=MAX_CERT_RADICAL_Q).contains(&q)
                    || p.unsigned_abs() > MAX_CERT_RADICAL_P as u64
                {
                    return None;
                }
                let b = rewrite_radicals(ctx, *base, groups)?;
                let (r, _) = radical_group(ctx, b, q, groups)?;
                return Some(match p {
                    1 => r,
                    _ => ctx.pow(r, ctx.num(p)),
                });
            }
            let b = rewrite_radicals(ctx, *base, groups)?;
            let e = rewrite_radicals(ctx, *exp, groups)?;
            Some(ctx.pow(b, e))
        }
    }
}

/// Find or create the generator for the `(base, q)` group. The generator is
/// a variable whose name (`$cert_r{i}`) the expression parser cannot
/// produce, so it can never collide with user input.
fn radical_group<'a>(
    ctx: &'a AtomArena<'a>,
    base: Atom<'a>,
    q: i64,
    groups: &mut Vec<(Atom<'a>, i64, Atom<'a>)>,
) -> Option<(Atom<'a>, usize)> {
    if let Some(idx) = groups.iter().position(|(b, gq, _)| *b == base && *gq == q) {
        return Some((groups[idx].2, idx));
    }
    if groups.len() >= MAX_CERT_RADICALS {
        return None;
    }
    let r = ctx.var(&format!("$cert_r{}", groups.len()));
    groups.push((base, q, r));
    Some((r, groups.len() - 1))
}

/// Reduce a sparse polynomial modulo the radical relations. Relations are
/// processed in reverse creation order: reducing `rᵢᵉ` (e ≥ q) replaces it
/// by `rᵢ^(e−q)·Sᵢ`, and `Sᵢ` only mentions *earlier* generators
/// (triangularity by construction), so one pass per relation, last-created
/// first, fully normalizes.
///
/// Returns `None` when the term cap trips.
fn reduce_mod_radicals(p: &Sparse, relations: &[RadicalRelation]) -> Option<Sparse> {
    let mut cur = p.clone();
    for rel in relations.iter().rev() {
        let q = rel.q as usize;
        let mut out = cur.zero();
        for (exp, coeff) in cur.terms_ref() {
            let e = exp[rel.var_index];
            if e < q {
                out = out.add(&Sparse::from_terms(
                    RationalDomain,
                    cur.n_vars(),
                    vec![(exp.to_vec(), coeff.clone())],
                ));
                continue;
            }
            let (a, b) = (e / q, e % q);
            let mut exp_red = exp.to_vec();
            exp_red[rel.var_index] = b;
            let mut term =
                Sparse::from_terms(RationalDomain, cur.n_vars(), vec![(exp_red, coeff.clone())]);
            for _ in 0..a {
                term = term.mul(&rel.base);
                if term.n_terms() > MAX_CERT_REDUCTION_TERMS {
                    return None;
                }
            }
            out = out.add(&term);
            if out.n_terms() > MAX_CERT_REDUCTION_TERMS {
                return None;
            }
        }
        cur = out;
    }
    Some(cur)
}

/// Layer 3: adjoin the radical relations `rᵢ^qᵢ = Sᵢ` and decide zero after
/// triangular reduction. Sound: every relation is an exact identity, so the
/// reduction only turns genuine zeros into zero (false negatives remain
/// possible when a relation is missing, e.g. `S^(1/2)` versus `S^(1/3)`
/// sharing no common root).
fn radical_layer<'a>(
    ctx: &'a AtomArena<'a>,
    expr: Atom<'a>,
    var: Symbol,
    max_gens: usize,
) -> Result<Certificate<'a>, CertDecline> {
    if node_count(expr) > MAX_CERT_RADICAL_NODES {
        return Err(CertDecline::Budget);
    }
    let mut groups: Vec<(Atom<'a>, i64, Atom<'a>)> = Vec::new();
    let Some(rewritten) = rewrite_radicals(ctx, expr, &mut groups) else {
        if std::env::var_os("OCAS_CERT_TIMING").is_some() {
            eprintln!("[cert] radical rewrite declined: {expr}");
        }
        return Err(CertDecline::Budget);
    };
    if std::env::var_os("OCAS_CERT_TIMING").is_some() {
        eprintln!(
            "[cert] radical layer input: {rewritten} (groups={})",
            groups.len()
        );
    }
    if groups.is_empty() {
        return Err(CertDecline::NotInField);
    }

    // Generator list: x, then the radical generators, then every other atom.
    let mut gens: Vec<Atom<'a>> = Vec::with_capacity(8);
    gens.push(ctx.var(var.as_str()));
    let mut overflow = false;
    for (_, _, r) in &groups {
        if !push_generator(&mut gens, *r, max_gens, &mut overflow) {
            return Err(CertDecline::Budget);
        }
    }
    if !collect_generators(rewritten, &mut gens, max_gens, &mut overflow) || overflow {
        return Err(if overflow {
            CertDecline::Budget
        } else {
            CertDecline::NotInField
        });
    }

    // Relations: each base must be a polynomial in the generators (a
    // rational right-hand side is outside this fragment).
    let mut relations: Vec<RadicalRelation> = Vec::with_capacity(groups.len());
    for (base, q, r) in &groups {
        let rf =
            atom_to_rational_extended(*base, &gens, gens.len()).ok_or(CertDecline::NotInField)?;
        if rf.denominator.n_terms() != 1
            || !RationalDomain.is_one(&rf.denominator.coeff(&vec![0; gens.len()]))
        {
            return Err(CertDecline::NotInField);
        }
        let var_index = gens
            .iter()
            .position(|g| g == r)
            .expect("generator registered");
        relations.push(RadicalRelation {
            var_index,
            q: *q as u32,
            base: rf.numerator,
        });
    }

    let rf =
        atom_to_rational_extended(rewritten, &gens, gens.len()).ok_or(CertDecline::NotInField)?;
    // Entry-size cap for the triangular reduction (0.30.0 W4a): the
    // per-step term cap (`MAX_CERT_REDUCTION_TERMS`) still lets a
    // multi-thousand-term numerator grind for ~14 s through 5 relations
    // (measured on a Möbius multi-radical answer — over the corpus
    // harness's per-case budget). Declining as `Budget` trades certificate
    // coverage for determinism, the same convention as the other caps.
    if rf.numerator.n_terms() > MAX_CERT_RADICAL_POLY_TERMS
        || rf.denominator.n_terms() > MAX_CERT_RADICAL_POLY_TERMS
    {
        return Err(CertDecline::Budget);
    }
    let num = reduce_mod_radicals(&rf.numerator, &relations).ok_or(CertDecline::Budget)?;
    if std::env::var_os("OCAS_CERT_TIMING").is_some() {
        eprintln!(
            "[cert] radical layer: gens={} relations={} num_terms={} den_terms={}",
            gens.len(),
            relations.len(),
            num.n_terms(),
            rf.denominator.n_terms()
        );
    }
    if num.is_zero() {
        return Ok(Certificate {
            method: CertMethod::RadicalField,
            difference: ctx.num(0),
        });
    }
    let den = reduce_mod_radicals(&rf.denominator, &relations).ok_or(CertDecline::Budget)?;
    if den.is_zero() {
        // The denominator vanishes in the quotient ring: no conclusion.
        return Err(CertDecline::NotInField);
    }
    Err(CertDecline::NonZero)
}

/// Whether `atom` is the literal zero.
///
/// `Atom`'s derived `PartialEq` compares node *contents* recursively, so a
/// head check keeps the zero test O(1) on large differences.
fn is_zero(atom: Atom<'_>) -> bool {
    matches!(atom.node(), AtomNode::Num(0))
}

/// Rewrite `I^k` by `I² = −1` (the imaginary unit is a `Var`, so the field
/// layer would otherwise treat it as an independent symbol and miss every
/// trigonometric identity).
fn reduce_imaginary_unit<'a>(ctx: &'a AtomArena<'a>, expr: Atom<'a>) -> Atom<'a> {
    match expr.node() {
        AtomNode::Num(_) | AtomNode::Var(_) => expr,
        AtomNode::Fun(name, args) => {
            let rebuilt: Vec<Atom<'a>> = args
                .iter()
                .map(|a| reduce_imaginary_unit(ctx, *a))
                .collect();
            let rebuilt = ctx.fun(name.as_str(), &rebuilt);
            if rebuilt == expr { expr } else { rebuilt }
        }
        AtomNode::Add(args) => {
            let rebuilt: Vec<Atom<'a>> = args
                .iter()
                .map(|a| reduce_imaginary_unit(ctx, *a))
                .collect();
            let rebuilt = ctx.add(&rebuilt);
            if rebuilt == expr { expr } else { rebuilt }
        }
        AtomNode::Mul(args) => {
            let mut power: i64 = 0;
            let mut rest: Vec<Atom<'a>> = Vec::with_capacity(args.len());
            let mut changed = false;
            for a in args.iter() {
                let a = reduce_imaginary_unit(ctx, *a);
                if let Some(n) = imaginary_exponent(a) {
                    power += n;
                    changed = true;
                } else {
                    rest.push(a);
                }
            }
            if !changed {
                return expr;
            }
            match power.rem_euclid(4) {
                0 => {}
                1 => rest.push(ctx.var("I")),
                2 => rest.push(ctx.num(-1)),
                _ => {
                    rest.push(ctx.num(-1));
                    rest.push(ctx.var("I"));
                }
            }
            ctx.mul(&rest)
        }
        AtomNode::Pow(base, exp) => {
            let base = reduce_imaginary_unit(ctx, *base);
            let exp = reduce_imaginary_unit(ctx, *exp);
            if let AtomNode::Num(n) = exp.node()
                && is_imaginary_unit(base)
            {
                return imaginary_power(ctx, *n);
            }
            let rebuilt = ctx.pow(base, exp);
            if rebuilt == expr { expr } else { rebuilt }
        }
    }
}

fn is_imaginary_unit(atom: Atom<'_>) -> bool {
    matches!(atom.node(), AtomNode::Var(s) if s.as_str() == "I")
}

/// `Some(n)` for `I` and `I^n` factors, so products can accumulate the
/// exponent.
fn imaginary_exponent(atom: Atom<'_>) -> Option<i64> {
    if is_imaginary_unit(atom) {
        return Some(1);
    }
    match atom.node() {
        AtomNode::Pow(base, exp) if is_imaginary_unit(*base) => match exp.node() {
            AtomNode::Num(n) => Some(*n),
            _ => None,
        },
        _ => None,
    }
}

/// `I^n` for any integer `n`, using `I² = −1`.
fn imaginary_power<'a>(ctx: &'a AtomArena<'a>, n: i64) -> Atom<'a> {
    match n.rem_euclid(4) {
        0 => ctx.num(1),
        1 => ctx.var("I"),
        2 => ctx.num(-1),
        _ => ctx.mul(&[ctx.num(-1), ctx.var("I")]),
    }
}

/// Maximum dependency-merge rounds before giving up on the field layer.
const MAX_CERT_MERGE_ROUNDS: usize = 8;

/// Express dependent `exp`/`log` atoms in terms of the others, so that the
/// free-generator field can see their relations.
///
/// Reuses the tower's dependency merge ([`crate::tower::merge`]): without it
/// a difference such as `exp(2x) − exp(x)²` (or `exp(x)·exp(−x) − 1`, which
/// every hyperbolic rewrite produces) has a nonzero normal form and would be
/// reported as `NonZero` — a false alarm.
fn merge_dependent_atoms<'a>(ctx: &'a AtomArena<'a>, expr: Atom<'a>, var: Symbol) -> Atom<'a> {
    let mut current = expr;
    for _ in 0..MAX_CERT_MERGE_ROUNDS {
        let mut gens: Vec<GenInfo<'a>> = Vec::new();
        let mut changed = false;
        // Register the *simplest argument first*: the loop only merges a
        // candidate into an already-registered generator, and the integer
        // class `exp(k·v) = exp(v)^k` (merge.rs, 0.29.0) rewrites the
        // *multiple* into the base. Registering `exp(2x)` before `exp(x)`
        // left the ratio `1/2` (fractional — no relation) and the field saw
        // both as independent: the tower proptest's
        // `exp(2x) + exp(x)·exp(2x)` counterexample.
        let mut candidates = collect_funs(current);
        candidates
            .sort_by_key(|(_, app)| app.children().first().map_or(0usize, |a| node_count(*a)));
        for (name, app) in candidates {
            if gens.iter().any(|g| g.atom == app) {
                continue;
            }
            let kind = match name.as_str() {
                "log" => GenKind::Log,
                "exp" => GenKind::Exp,
                _ => continue,
            };
            let args = app.children();
            if args.len() != 1 {
                continue;
            }
            if let Some(merge) = merge::merge_candidate(ctx, kind, args[0], app, &gens, var) {
                current = merge::substitute_atom(ctx, current, app, merge.replacement);
                changed = true;
                break;
            }
            gens.push(GenInfo {
                kind,
                atom: app,
                arg: args[0],
                dt: KElem::zero(0),
            });
        }
        if !changed {
            break;
        }
    }
    current
}

/// Estimate the number of additive terms `expr` expands to, saturating at
/// [`MAX_CERT_EXPANSION`]-plus-one (i.e. "too big").
///
/// Function applications and variables count as one term — they are the
/// field's generators, not sums. A non-integer or oversized exponent cannot
/// be estimated and is reported as "too big", which makes the caller decline
/// the expanding layers.
fn expansion_size(expr: Atom<'_>) -> u64 {
    let cap = MAX_CERT_EXPANSION.saturating_add(1);
    match expr.node() {
        AtomNode::Num(_) | AtomNode::Var(_) | AtomNode::Fun(_, _) => 1,
        AtomNode::Add(args) => args
            .iter()
            .fold(0u64, |acc, a| acc.saturating_add(expansion_size(*a)))
            .min(cap),
        AtomNode::Mul(args) => args
            .iter()
            .fold(1u64, |acc, a| acc.saturating_mul(expansion_size(*a)))
            .min(cap),
        AtomNode::Pow(base, exp) => {
            let AtomNode::Num(n) = exp.node() else {
                return cap;
            };
            let power = n.unsigned_abs();
            if power > MAX_CERT_EXPONENT {
                return cap;
            }
            let base_size = expansion_size(*base);
            let mut acc = 1u64;
            for _ in 0..power {
                acc = acc.saturating_mul(base_size).min(cap);
            }
            acc
        }
    }
}

/// Collect the independent generators of `expr` into `gens`.
///
/// Order does not matter for correctness (generators are matched by whole
/// atom equality), but inner atoms are added after the atoms that contain
/// them so a generator is never shadowed by its children.
///
/// Returns `false` when a node makes the expression non-rational in any
/// possible generator set (a `Derivative(…)` or `Integral(…)` head, or a
/// non-integer power whose base is itself non-constant — the latter is the
/// radical layer's job, not the elementary layer's).
fn collect_generators<'a>(
    expr: Atom<'a>,
    gens: &mut Vec<Atom<'a>>,
    max_gens: usize,
    overflow: &mut bool,
) -> bool {
    match expr.node() {
        AtomNode::Num(_) => true,
        AtomNode::Var(_) => push_generator(gens, expr, max_gens, overflow),
        AtomNode::Fun(name, args) => {
            let n = name.as_str();
            if n == "Integral" || n == "Derivative" {
                return false;
            }
            if !push_generator(gens, expr, max_gens, overflow) {
                return false;
            }
            args.iter()
                .all(|a| collect_generators(*a, gens, max_gens, overflow))
        }
        AtomNode::Add(args) | AtomNode::Mul(args) => args
            .iter()
            .all(|a| collect_generators(*a, gens, max_gens, overflow)),
        AtomNode::Pow(base, exp) => {
            if let AtomNode::Num(_) = exp.node() {
                // An integer power of a rational expression stays rational.
                collect_generators(*base, gens, max_gens, overflow)
                    && collect_generators(*exp, gens, max_gens, overflow)
            } else if is_rational_constant(*base) {
                // `sqrt(2)`, `2^(1/3)`, … — a constant; adjoin it whole so
                // the field sees an independent symbol (sound: only false
                // negatives are possible).
                push_generator(gens, expr, max_gens, overflow)
            } else {
                // A genuine radical of a non-constant base: the elementary
                // layer declines; the radical layer (0.30.0) adjoins the
                // relation `r^q = base` instead.
                false
            }
        }
    }
}

/// Register `atom` as a generator unless it is already present; reports the
/// generator budget through `overflow`.
fn push_generator<'a>(
    gens: &mut Vec<Atom<'a>>,
    atom: Atom<'a>,
    max_gens: usize,
    overflow: &mut bool,
) -> bool {
    if gens.contains(&atom) {
        return true;
    }
    if gens.len() >= max_gens {
        *overflow = true;
        return false;
    }
    gens.push(atom);
    true
}

/// Whether `atom` is a rational constant (contains no variables at all).
fn is_rational_constant(atom: Atom<'_>) -> bool {
    match atom.node() {
        AtomNode::Num(_) => true,
        AtomNode::Var(_) | AtomNode::Fun(_, _) => false,
        AtomNode::Add(args) | AtomNode::Mul(args) => args.iter().all(|a| is_rational_constant(*a)),
        AtomNode::Pow(b, e) => is_rational_constant(*b) && is_rational_constant(*e),
    }
}

#[cfg(test)]
mod tests {
    use ocas_core::arena::Arena;
    use ocas_parse::parse;

    use super::*;

    fn certify_str<'a>(
        ctx: &'a AtomArena<'a>,
        f: &str,
        big_f: &str,
    ) -> Result<Certificate<'a>, CertDecline> {
        let integrand = normalize(ctx, parse(ctx, f).expect("parse integrand"));
        let antiderivative = parse(ctx, big_f).expect("parse antiderivative");
        certify(ctx, integrand, antiderivative, Symbol::new("x"))
    }

    #[test]
    fn structural_layer_certifies_polynomial_pairs() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        // ∫ 2·x dx = x²
        let cert = certify_str(&ctx, "2*x", "x^2").expect("certificate");
        assert_eq!(cert.method, CertMethod::Structural);
        assert_eq!(cert.difference.to_string(), "0");
    }

    #[test]
    fn structural_layer_collects_like_terms() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        // ∫ x dx = x²/2. The arena keeps `p·q⁻¹` as a product, so the
        // difference is `x − 2·(2⁻¹)·x`; `collect_terms` folds it to zero.
        let cert = certify_str(&ctx, "x", "(1/2)*x^2").expect("certificate");
        assert_eq!(cert.method, CertMethod::Structural);
    }

    #[test]
    fn field_layer_certifies_a_sum_of_logs() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        // ∫ 1/(x·(x+1)) dx = log(x) − log(x+1). The structural layer sees
        // the three separate terms; the field layer puts them over a common
        // denominator and sees zero.
        let cert = certify_str(&ctx, "1/(x*(x+1))", "log(x)-log(x+1)").expect("certificate");
        assert_eq!(cert.method, CertMethod::ElementaryField);
    }

    #[test]
    fn trig_heads_are_rewritten_before_the_field_layer() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        // ∫ cos(x) dx = sin(x) — the derivative table returns cos(x) but the
        // normalization of the difference may need the elementary rewrite.
        let cert = certify_str(&ctx, "cos(x)", "sin(x)").expect("certificate");
        assert!(matches!(
            cert.method,
            CertMethod::Structural | CertMethod::ElementaryField
        ));
    }

    #[test]
    fn wrong_antiderivative_is_declined() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        // D(x²) = 2x ≠ x.
        assert_eq!(
            certify_str(&ctx, "x", "x^2").unwrap_err(),
            CertDecline::NonZero
        );
        // A dropped constant factor.
        assert!(certify_str(&ctx, "2*sin(x)", "sin(x)").is_err());
    }

    #[test]
    fn residue_carrying_answers_are_not_certifiable() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        assert_eq!(
            certify_str(&ctx, "1/(x^2+5)", "Integral(1/(x^2+5), x)").unwrap_err(),
            CertDecline::NotInField
        );
    }

    #[test]
    fn budget_is_reported_separately() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        let integrand = normalize(&ctx, parse(&ctx, "x").unwrap());
        let antiderivative = parse(&ctx, "log(x)-log(x+1)").unwrap();
        // A one-node ceiling forces the budget decline before any layer runs.
        assert_eq!(
            certify_with_limits(&ctx, integrand, antiderivative, Symbol::new("x"), 1, 16)
                .unwrap_err(),
            CertDecline::Budget
        );
    }

    #[test]
    fn generator_budget_is_reported() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        // D(x) = 1 ≠ log(x): the difference `1 − log(x)` keeps a logarithm
        // atom, so the field layer needs a second generator.
        let integrand = normalize(&ctx, parse(&ctx, "log(x)").unwrap());
        let antiderivative = parse(&ctx, "x").unwrap();
        // Zero generators allowed beyond `x`: the generator collection
        // overflows before the field is built.
        assert_eq!(
            certify_with_limits(&ctx, integrand, antiderivative, Symbol::new("x"), 20_000, 1)
                .unwrap_err(),
            CertDecline::Budget
        );
    }

    #[test]
    fn radical_of_a_non_constant_base_certifies_via_radical_layer() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        // ∫ 1/sqrt(1−x²) dx = asin(x): the difference may carry
        // (1−x²)^(−1/2) — a radical of a non-constant base. Before the
        // radical layer this declined; now the relation r² = 1−x² decides it.
        let verdict = certify_str(&ctx, "1/sqrt(1-x^2)", "asin(x)");
        assert!(verdict.is_ok(), "expected a certificate, got: {verdict:?}");
    }

    #[test]
    fn radical_layer_uses_the_constant_radical_relation() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        // ∫ dx/(x²−2) = (1/√8)·log((2x−√8)/(2x+√8)): the difference needs
        // √8·√8 = 8, which the free-generator field cannot see — the answer
        // can only be certified with the relation r² = 8 adjoined.
        let cert = certify_str(
            &ctx,
            "1/(x^2-2)",
            "8^(-1/2)*log((2*x-8^(1/2))/(2*x+8^(1/2)))",
        )
        .expect("certificate");
        assert_eq!(cert.method, CertMethod::RadicalField);
    }

    #[test]
    fn radical_layer_merges_sqrt_head_and_fractional_power() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        // D(−sqrt(1−x²)) = x/sqrt(1−x²): the derivative table emits the
        // `sqrt` head while the integrand carries the fractional power —
        // both rewrite to the same (base, 2) generator.
        let cert = certify_str(&ctx, "x*(1-x^2)^(-1/2)", "-sqrt(1-x^2)").expect("certificate");
        assert!(matches!(
            cert.method,
            CertMethod::RadicalField | CertMethod::ElementaryField | CertMethod::Structural
        ));
    }

    #[test]
    fn merge_orders_integer_multiples_after_their_base() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        // Regression (tower proptest counterexample, 0.30.0 W5b review):
        // `∫ (exp(2x) + exp(x)·exp(2x)) dx = exp(2x)/2 + exp(x)³/3`. The
        // difference needs `exp(2x) = exp(x)²`; registering `exp(2x)` first
        // (source order) left the ratio 1/2 and the field saw two
        // independent generators, so the certificate declined.
        let cert = certify_str(
            &ctx,
            "exp(2*x) + exp(x)*exp(2*x)",
            "exp(2*x)/2 + exp(x)^3/3",
        )
        .expect("certificate");
        assert!(matches!(
            cert.method,
            CertMethod::Structural | CertMethod::ElementaryField
        ));
    }

    #[test]
    fn radical_layer_rejects_wrong_answers() {
        let arena = Arena::new();
        let ctx = AtomArena::new(&arena);
        // A wrong coefficient inside the log argument is not a relation
        // artifact: no layer may certify it.
        assert!(
            certify_str(
                &ctx,
                "1/(x^2-2)",
                "8^(-1/2)*log((3*x-8^(1/2))/(2*x+8^(1/2)))"
            )
            .is_err()
        );
    }
}