diff --git a/problemreductions-expr/tests/fixtures/sympy_oracle.json b/problemreductions-expr/tests/fixtures/sympy_oracle.json index 999eaf1a4..27feff24d 100644 --- a/problemreductions-expr/tests/fixtures/sympy_oracle.json +++ b/problemreductions-expr/tests/fixtures/sympy_oracle.json @@ -849,6 +849,111 @@ "right": "exp(n)", "ratio_limit": "0", "relation": "right_dominates" + }, + { + "name": "polynomial_exponent", + "left": "2^(n^2)", + "right": "3^n", + "ratio_limit": "oo", + "relation": "left_dominates" + }, + { + "name": "sqrt_exponent_over_polynomial", + "left": "2^sqrt(n)", + "right": "n^100", + "ratio_limit": "oo", + "relation": "left_dominates" + }, + { + "name": "exponential_over_sqrt_exponent", + "left": "1.0001^n", + "right": "2^sqrt(n)", + "ratio_limit": "oo", + "relation": "left_dominates" + }, + { + "name": "quasi_polynomial", + "left": "2^(log(n)^2)", + "right": "n^100", + "ratio_limit": "oo", + "relation": "left_dominates" + }, + { + "name": "self_power_over_binary_power", + "left": "n^n", + "right": "2^(n * log(n))", + "ratio_limit": "oo", + "relation": "left_dominates" + }, + { + "name": "variable_base_rate", + "left": "n^(2 * n)", + "right": "n^n", + "ratio_limit": "oo", + "relation": "left_dominates" + }, + { + "name": "mixed_base_product", + "left": "2^(2 * n) * 3^n", + "right": "2^n * 4^n", + "ratio_limit": "oo", + "relation": "left_dominates" + }, + { + "name": "close_exponential_rates", + "left": "2^(0.773 * n)", + "right": "1.1996^n", + "ratio_limit": "oo", + "relation": "left_dominates" + }, + { + "name": "exponent_cancellation", + "left": "2^((n + 1)^2 - n^2)", + "right": "2^(2 * n)", + "ratio_limit": "2", + "relation": "equivalent" + }, + { + "name": "decaying_factor", + "left": "2^n * 0.5^n", + "right": "1", + "ratio_limit": "1", + "relation": "equivalent" + }, + { + "name": "reverse_variable_base", + "left": "2^n", + "right": "n^n", + "ratio_limit": "0", + "relation": "right_dominates" + }, + { + "name": "stirling_factorial", + "left": "factorial(n)", + "right": "n^n * sqrt(n) * exp(-n)", + "ratio_limit": "sqrt(2)*sqrt(pi)", + "relation": "equivalent" + }, + { + "name": "factorial_below_self_power", + "left": "n^n", + "right": "factorial(n)", + "ratio_limit": "oo", + "relation": "left_dominates" + }, + { + "name": "factorial_above_base_e_power", + "left": "factorial(n)", + "right": "(n / 3)^n", + "ratio_limit": "oo", + "relation": "left_dominates" + }, + { + "name": "scaled_factorial", + "left": "factorial(2 * n)", + "right": "(2 * n)^(2 * n) * sqrt(n) * exp(-2 * n)", + "ratio_limit": "2*sqrt(pi)", + "relation": "equivalent" } ], "factorial_domain_cases": [ diff --git a/scripts/generate_symbolic_expr_fixture.py b/scripts/generate_symbolic_expr_fixture.py index 082136b82..52bf21f45 100644 --- a/scripts/generate_symbolic_expr_fixture.py +++ b/scripts/generate_symbolic_expr_fixture.py @@ -148,6 +148,21 @@ ("exponential_poly_tie_break", "2^n * n", "2^n"), ("reverse_polynomial_degree", "n", "n^2"), ("reverse_exponential", "n^100", "exp(n)"), + ("polynomial_exponent", "2^(n^2)", "3^n"), + ("sqrt_exponent_over_polynomial", "2^sqrt(n)", "n^100"), + ("exponential_over_sqrt_exponent", "1.0001^n", "2^sqrt(n)"), + ("quasi_polynomial", "2^(log(n)^2)", "n^100"), + ("self_power_over_binary_power", "n^n", "2^(n * log(n))"), + ("variable_base_rate", "n^(2 * n)", "n^n"), + ("mixed_base_product", "2^(2 * n) * 3^n", "2^n * 4^n"), + ("close_exponential_rates", "2^(0.773 * n)", "1.1996^n"), + ("exponent_cancellation", "2^((n + 1)^2 - n^2)", "2^(2 * n)"), + ("decaying_factor", "2^n * 0.5^n", "1"), + ("reverse_variable_base", "2^n", "n^n"), + ("stirling_factorial", "factorial(n)", "n^n * sqrt(n) * exp(-n)"), + ("factorial_below_self_power", "n^n", "factorial(n)"), + ("factorial_above_base_e_power", "factorial(n)", "(n / 3)^n"), + ("scaled_factorial", "factorial(2 * n)", "(2 * n)^(2 * n) * sqrt(n) * exp(-2 * n)"), ] diff --git a/src/growth.rs b/src/growth.rs index 10cb97aae..ed6206ee9 100644 --- a/src/growth.rs +++ b/src/growth.rs @@ -1,5 +1,5 @@ //! Symbolic growth domain: a dedicated asymptotic normal form for reduction -//! parameter expressions. +//! parameter expressions and complexity bounds. //! //! Where full monomial canonicalization answers Big-O questions by expanding an //! [`Expr`] to monomial normal form, with exponential cost in nesting depth, the @@ -10,24 +10,39 @@ //! //! # Representation //! -//! One internal growth term is a monomial +//! One internal growth term is stored in log space, //! //! ```text -//! ∏_v ∏_f base[f]^(coefficient[f] · v) -//! · ∏_v v^(poly[v]) · ∏_v (log v)^(logs[v]) +//! exp( Σ_atom weight[atom] · atom ), atom ∈ { ∏_v v^p · (ln v)^l , ln ln v } +//! weight = Σ_b c_b · ln b (exact rational c_b, symbolic base b) //! ``` //! -//! and a [`Growth`] is either a known antichain of pairwise-incomparable dominant -//! terms (each summand of an asymptotic sum), or an unknown result with explicit -//! reasons for content we cannot represent symbolically. +//! so `n^2` is the atom `ln n` with weight `2·ln e`, `2^(n·m)` is the atom +//! `n·m` with weight `ln 2`, `k^n` is the atom `n·ln k` with weight `ln e`, and +//! `(log n)^3` is the atom `ln ln n` with weight `3·ln e`. A [`Growth`] is +//! either a known antichain of pairwise-incomparable dominant terms (each +//! summand of an asymptotic sum), or an unknown result with explicit reasons +//! for content we cannot represent symbolically. //! //! # Semantic foundation (the trust contract) //! //! Every expression admitted to the domain is assumed **nonnegative** and -//! **weakly monotone** (nondecreasing in each variable) on `vars ≥ 2`. Under -//! these axioms Howell's multivariate-O inconsistencies vanish and -//! `f + g ≍ max(f, g)` up to a constant factor, which licenses -//! `add = antichain union + prune`. All bounds produced are **upper** bounds. +//! **weakly monotone** (nondecreasing in each variable) on `vars ≥ 2`, and +//! asymptotic statements quantify over all variables growing together +//! (`∃N: ∀ vars ≥ N`). Under these axioms `f + g ≍ max(f, g)` up to a constant +//! factor, which licenses `add = antichain union + prune`. Every known growth is +//! an **upper** bound on the expression's magnitude, including subexpressions +//! that are negative (a subtracted summand, `log` of a decaying term); +//! [`Growth::relation`] reports +//! [`ParameterRelation::Exact`] when it is also a lower bound (`f = Θ(growth)`) +//! and [`ParameterRelation::UpperBound`] when any step below widened it. +//! +//! Exactness is preserved by sums and products of exact operands, constant +//! powers, exactly expanded exponents with exactly bounded bases, logarithms +//! whose argument grows without bound, and sums whose widened summands are +//! strictly dominated by an exact summand (so `n − 1` and `n² − n` stay exact). +//! `factorial(x)` is exact by Stirling, `x! ≍ x^(x+1/2)·e^(−x)`, when `x` +//! expands exactly to one monomial. //! //! Widening (always toward a valid upper bound): //! - Subtraction is normalized to addition of a negative term, and @@ -35,26 +50,42 @@ //! This also covers the //! `sqrt((a − b)^2)` absolute-value idiom (`|a − b| ≤ a + b`). //! - Constants and constant multipliers/divisors are dropped on entry. -//! - Exponentials with a **linear** exponent (`c^x`, `c^(r·x)`, `exp(x)`) are -//! first-class via symbolic base/coefficient factors. The original base is -//! authoritative: it is never normalized through a floating-point logarithm -//! and never reconstructed by rounding. Nonlinear exponents (`2^(n·k)`, -//! `2^sqrt(n)`), `factorial(·)`, and negative polynomial exponents widen to -//! an unknown result, which preserves its reasons through every operation. +//! - Exponents are expanded exactly into polynomials over `v` and `ln v`; +//! fractional powers of sums and logarithms of sums widen to upper bounds, +//! which are only ever combined monotonically. Exponential bases are +//! authoritative symbolic data: they are never normalized through a +//! floating-point logarithm and never reconstructed by rounding. +//! - A variable base `B` with exact expansion `Σ c_i μ_i` is bounded by +//! `C · M`, where `C = Σ c_i` and `M` takes the largest power of every +//! variable, so `B^E ≤ C^E · M^E`. Any other `factorial(x)` widens to `x^x`. +//! - A logarithm of an argument with only an upper bound is unknown: bounding +//! its magnitude also requires a lower bound on its positive argument. +//! - `ln ln v` inside an exponent widens to `ln v`. //! - The explicit approximation boundary treats [`Expr::log`] as the natural //! logarithm, but all fixed //! logarithm bases greater than one have the same asymptotic class and are //! intentionally represented by the single `log(v)` factor. //! +//! # Comparison +//! +//! Two terms are compared through their log-space difference `Q`. A weight's +//! sign is decided directly when every summand has the same sign. Mixed signs +//! use exact rational powers with a certified rational interval for `e`; +//! overlapping intervals or exhausted arithmetic budgets remain unknown. `Q` is +//! bounded below when every atom without a proven nonnegative weight is +//! strictly dominated (componentwise in `(power, log power)`, strictly in one +//! variable) by an atom with a proven positive weight. +//! //! # `Pow` note //! //! `Pow(base, k)` for a nonnegative constant `k` raises **each** antichain term -//! of `base` to the power `k` (scaling its exponents). This is the tight +//! of `base` to the power `k` (scaling its weights). This is the tight //! asymptotic answer — `(n + m)^2 ≍ max(n, m)^2 = max(n^2, m^2)` by AM-GM, so no //! binomial cross term is introduced — and it is what makes the widening chain //! `sqrt((n − m)^2) ≍ n + m` hold exactly. use crate::expr::{AlgebraicAnalysis, BigInt, Expr, ExprNode, ExprNodeId}; +use crate::parameters::ParameterRelation; use num_rational::BigRational; use num_traits::{One, Signed, ToPrimitive, Zero}; use std::cmp::Ordering; @@ -62,6 +93,12 @@ use std::collections::{BTreeMap, HashMap}; /// Maximum number of incomparable terms retained in one Big-O normal form. const ANTICHAIN_CAP: usize = 32; +/// Maximum number of monomials in one expanded exponent. +const EXPONENT_TERM_CAP: usize = 64; +/// Largest integer power of a sum expanded exactly inside an exponent. +const EXPONENT_POWER_CAP: u32 = 16; +/// Bit budget for exact rational powers (weight signs and exponent coefficients). +const RATIONAL_POWER_BIT_CAP: u64 = 1 << 16; /// An exact fixed exponential base. #[derive(Clone, Debug, PartialEq, Eq, PartialOrd, Ord)] @@ -73,241 +110,212 @@ enum ExpBase { } impl ExpBase { - fn directly_comparable_value(&self) -> Option<&BigRational> { - match self { - ExpBase::Rational(base) => Some(base), - ExpBase::Natural => None, - } - } - - fn direction(&self) -> Ordering { + /// Sign of `ln(base)`. + fn ln_sign(&self) -> Ordering { match self { ExpBase::Rational(base) => base.cmp(&BigRational::one()), ExpBase::Natural => Ordering::Greater, } } - fn coefficient_cmp(&self, a: &BigRational, b: &BigRational) -> Ordering { - if self.direction() == Ordering::Greater { - a.cmp(b) - } else { - a.cmp(b).reverse() + fn rational_bounds(&self) -> (BigRational, BigRational) { + match self { + ExpBase::Rational(base) => (base.clone(), base.clone()), + ExpBase::Natural => { + // S = sum_{k=0}^{10} 1/k! < e. The remaining series starts at + // 1/11!, with each later ratio at most 1/12, so its sum is + // at most 12/(11 * 11!). Cache these exact rational bounds. + static BOUNDS: std::sync::LazyLock<(BigRational, BigRational)> = + std::sync::LazyLock::new(|| { + let mut factorial = BigInt::one(); + let mut lower = BigRational::one(); + for k in 1..=10 { + factorial *= k; + lower += BigRational::new(BigInt::one(), factorial.clone()); + } + let tail = BigRational::new(BigInt::from(12), factorial * 121); + (lower.clone(), lower + tail) + }); + BOUNDS.clone() + } } } } -/// Exponential, polynomial, and logarithmic growth associated with one size -/// variable. Missing components have exponent zero. +/// `∏_v v^power · (ln v)^log_power`; absent variables have both exponents zero. +type Monomial = BTreeMap, (BigRational, BigRational)>; + +/// One log-space atom of a growth term. #[derive(Clone, Debug, PartialEq, Eq, PartialOrd, Ord)] -struct VariableGrowth { - exp: BTreeMap, - poly: BigRational, - log: u32, +enum Atom { + Mono(Monomial), + /// `ln ln v`. + LogLog(Box), } -impl VariableGrowth { - fn empty() -> Self { - Self { - exp: BTreeMap::new(), - poly: BigRational::zero(), - log: 0, - } - } +/// `Σ_b c_b · ln b` with exact coefficients and symbolic bases. +type Weight = BTreeMap; - fn exponential(base: ExpBase, coefficient: BigRational) -> Self { - Self { - exp: BTreeMap::from([(base, coefficient)]), - poly: BigRational::zero(), - log: 0, - } - } +fn single_variable(variable: &str, power: BigRational, log_power: BigRational) -> Monomial { + BTreeMap::from([(variable.into(), (power, log_power))]) +} - fn polynomial(degree: BigRational) -> Self { - Self { - exp: BTreeMap::new(), - poly: degree, - log: 0, - } +fn ln_atom(variable: &str) -> Atom { + Atom::Mono(single_variable( + variable, + BigRational::zero(), + BigRational::one(), + )) +} + +fn natural(coefficient: BigRational) -> Weight { + BTreeMap::from([(ExpBase::Natural, coefficient)]) +} + +fn monomial_mul(a: &Monomial, b: &Monomial) -> Monomial { + let mut result = a.clone(); + for (variable, (power, log_power)) in b { + let entry = result + .entry(variable.clone()) + .or_insert_with(|| (BigRational::zero(), BigRational::zero())); + entry.0 += power; + entry.1 += log_power; } + result.retain(|_, (power, log_power)| !power.is_zero() || !log_power.is_zero()); + result +} - fn logarithmic(power: u32) -> Self { - Self { - exp: BTreeMap::new(), - poly: BigRational::zero(), - log: power, +fn monomial_pow(a: &Monomial, k: &BigRational) -> Monomial { + a.iter() + .map(|(variable, (power, log_power))| (variable.clone(), (power * k, log_power * k))) + .collect() +} + +/// `b = o(a)` as all variables grow: `a ≥ b` per variable in `(power, log)` +/// lexicographic order, strictly for at least one variable. +fn monomial_dominates(a: &Monomial, b: &Monomial) -> bool { + let zero = (BigRational::zero(), BigRational::zero()); + let mut strict = a.keys().any(|variable| !b.contains_key(variable)); + for (variable, y) in b { + let x = a.get(variable).unwrap_or(&zero); + match x.0.cmp(&y.0).then(x.1.cmp(&y.1)) { + Ordering::Less => return false, + Ordering::Greater => strict = true, + Ordering::Equal => {} } } + strict +} - fn is_empty(&self) -> bool { - self.exp.is_empty() && self.poly.is_zero() && self.log == 0 +fn atom_dominates(a: &Atom, b: &Atom) -> bool { + match (a, b) { + (Atom::Mono(a), Atom::Mono(b)) => monomial_dominates(a, b), + (Atom::Mono(a), Atom::LogLog(variable)) => a + .get(variable) + .is_some_and(|(power, log_power)| power.is_positive() || log_power.is_positive()), + (Atom::LogLog(_), _) => false, } +} - fn mul(&self, other: &Self) -> Result { - let mut result = self.clone(); - for (base, coefficient) in &other.exp { - *result - .exp - .entry(base.clone()) - .or_insert_with(BigRational::zero) += coefficient; - } - result.exp.retain(|_, coefficient| !coefficient.is_zero()); - result.poly += &other.poly; - result.log = result.log.checked_add(other.log).ok_or_else(|| { - GrowthFailure::RepresentedExponentOutOfRange(format!("{} + {}", self.log, other.log)) - })?; - Ok(result) +fn weight_add(target: &mut Weight, other: &Weight, scale: &BigRational) { + for (base, coefficient) in other { + *target.entry(base.clone()).or_insert_with(BigRational::zero) += coefficient * scale; } + target.retain(|_, coefficient| !coefficient.is_zero()); +} - fn pow(&self, power: &BigRational) -> Result { - let exp = self - .exp - .iter() - .filter_map(|(base, coefficient)| { - let coefficient = coefficient * power; - (!coefficient.is_zero()).then(|| (base.clone(), coefficient)) - }) - .collect(); - let poly = &self.poly * power; - let scaled_log = BigRational::from_integer(BigInt::from(self.log)) * power; - let rounded = - (scaled_log.numer() + scaled_log.denom() - BigInt::one()) / scaled_log.denom(); - let Some(log) = rounded.to_u32() else { - return Err(GrowthFailure::RepresentedExponentOutOfRange( - scaled_log.to_string(), - )); - }; - Ok(Self { exp, poly, log }) - } - - fn cmp_exp(&self, other: &Self) -> Option { - let mut left_count = 0; - let mut right_count = 0; - let mut left_single: Option<(&ExpBase, BigRational)> = None; - let mut right_single: Option<(&ExpBase, BigRational)> = None; - - for (base, a) in &self.exp { - if let Some(b) = other.exp.get(base) { - match base.coefficient_cmp(a, b) { - Ordering::Equal => {} - Ordering::Greater => { - left_count += 1; - left_single = Some((base, a - b)); - } - Ordering::Less => { - right_count += 1; - right_single = Some((base, b - a)); - } - } - } else { - left_count += 1; - left_single = Some((base, a.clone())); - } - } +fn weight_scale(weight: &Weight, scale: &BigRational) -> Weight { + let mut result = Weight::new(); + weight_add(&mut result, weight, scale); + result +} - for (base, coefficient) in &other.exp { - if !self.exp.contains_key(base) { - right_count += 1; - right_single = Some((base, coefficient.clone())); - } - } +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +enum Sign { + Negative, + Zero, + Positive, + Unknown, +} - match (left_count, right_count) { - (0, 0) => Some(Ordering::Equal), - (0, _) => Some(Ordering::Less), - (_, 0) => Some(Ordering::Greater), - (1, 1) => { - let (a_base, a_coefficient) = left_single?; - let (b_base, b_coefficient) = right_single?; - Self::cmp_single_factor(a_base, &a_coefficient, b_base, &b_coefficient) - } - _ => None, +/// Decide the sign of `Σ_b c_b · ln b`. +fn weight_sign(weight: &Weight) -> Sign { + let (mut positive, mut negative) = (false, false); + for (base, coefficient) in weight { + let ln_sign = base.ln_sign(); + if ln_sign == Ordering::Equal || coefficient.is_zero() { + continue; } - } - - fn cmp_growth(&self, other: &Self) -> Option { - let exponential = self.cmp_exp(other)?; - Some(if exponential == Ordering::Equal { - self.poly.cmp(&other.poly).then(self.log.cmp(&other.log)) + if (ln_sign == Ordering::Greater) == coefficient.is_positive() { + positive = true; } else { - exponential - }) - } - - fn cmp_single_factor( - a_base: &ExpBase, - a_coefficient: &BigRational, - b_base: &ExpBase, - b_coefficient: &BigRational, - ) -> Option { - if a_base == b_base { - return Some(a_base.coefficient_cmp(a_coefficient, b_coefficient)); - } - - if a_coefficient == b_coefficient { - match (a_base, b_base) { - (ExpBase::Natural, ExpBase::Rational(_)) => { - let base = b_base.directly_comparable_value()?; - if base <= &BigRational::from_integer(2.into()) { - return Some(Ordering::Greater); - } - if base >= &BigRational::from_integer(3.into()) { - return Some(Ordering::Less); - } - return None; - } - (ExpBase::Rational(_), ExpBase::Natural) => { - return Self::cmp_single_factor(b_base, b_coefficient, a_base, a_coefficient) - .map(Ordering::reverse); - } - _ => {} - } + negative = true; } + } + match (positive, negative) { + (false, false) => return Sign::Zero, + (true, false) => return Sign::Positive, + (false, true) => return Sign::Negative, + (true, true) => {} + } - let (a_base, b_base) = ( - a_base.directly_comparable_value()?, - b_base.directly_comparable_value()?, - ); - if a_coefficient == b_coefficient { - let base_order = a_base.cmp(b_base); - return if a_coefficient.is_positive() { - Some(base_order) - } else { - Some(base_order.reverse()) - }; - } + exact_weight_sign(weight).unwrap_or(Sign::Unknown) +} - if a_base > &BigRational::one() && b_base > &BigRational::one() { - match (a_base.cmp(b_base), a_coefficient.cmp(b_coefficient)) { - (Ordering::Greater | Ordering::Equal, Ordering::Greater | Ordering::Equal) => { - Some(Ordering::Greater) - } - (Ordering::Less | Ordering::Equal, Ordering::Less | Ordering::Equal) => { - Some(Ordering::Less) - } - _ => None, - } - } else if a_base < &BigRational::one() && b_base < &BigRational::one() { - match (a_base.cmp(b_base), a_coefficient.cmp(b_coefficient)) { - (Ordering::Less | Ordering::Equal, Ordering::Less | Ordering::Equal) => { - Some(Ordering::Greater) - } - (Ordering::Greater | Ordering::Equal, Ordering::Greater | Ordering::Equal) => { - Some(Ordering::Less) - } - _ => None, - } +/// Compare `∏ b^(c_b·D)` with one, using exact rational interval bounds +/// and the common denominator `D`. Overlapping intervals remain unproved. +fn exact_weight_sign(weight: &Weight) -> Option { + let mut denominator = BigInt::one(); + for coefficient in weight.values() { + let d = coefficient.denom(); + denominator = &denominator * BigRational::new(denominator.clone(), d.clone()).denom(); + } + let mut bits = 0_u64; + let mut powers = Vec::new(); + for (base, coefficient) in weight { + let (lower, upper) = base.rational_bounds(); + let power = (coefficient * BigRational::from_integer(denominator.clone())).to_integer(); + let magnitude = power.abs().to_u64()?; + let base_bits = (lower.numer().bits() + lower.denom().bits()) + .max(upper.numer().bits() + upper.denom().bits()); + bits = bits.checked_add(magnitude.checked_mul(base_bits)?)?; + powers.push((lower, upper, power)); + } + if bits > RATIONAL_POWER_BIT_CAP { + return None; + } + let (mut above_lower, mut above_upper, mut below_lower, mut below_upper) = ( + BigRational::one(), + BigRational::one(), + BigRational::one(), + BigRational::one(), + ); + for (lower, upper, power) in powers { + let exponent = power.abs().to_i32()?; + if power.is_positive() { + above_lower *= lower.pow(exponent); + above_upper *= upper.pow(exponent); } else { - None + below_lower *= lower.pow(exponent); + below_upper *= upper.pow(exponent); } } + if above_lower > below_upper { + Some(Sign::Positive) + } else if above_upper < below_lower { + Some(Sign::Negative) + } else if above_lower == above_upper && below_lower == below_upper && above_lower == below_lower + { + Some(Sign::Zero) + } else { + None + } } -/// One growth monomial, e.g. `2^(3k) · n^2 · m · log(n)`. -/// -/// Empty maps represent `O(1)`. +/// One growth term, `exp(Σ weight · atom)`. An empty map represents `O(1)`. #[derive(Clone, Debug, PartialEq, Eq, PartialOrd, Ord)] struct GrowthTerm { - variables: BTreeMap, VariableGrowth>, + atoms: BTreeMap, } /// The asymptotic growth class of an [`Expr`]. @@ -316,7 +324,11 @@ pub struct Growth(GrowthState); #[derive(Clone, Debug, PartialEq, Eq)] enum GrowthState { - Known(Vec), + /// A pruned antichain and whether it is also a lower bound. + Known { + terms: Vec, + relation: ParameterRelation, + }, /// Content outside the represented growth domain, with every reason that /// contributed to the result. Unknown(Vec), @@ -329,22 +341,20 @@ pub enum GrowthFailure { InvalidConstantDomain { expression: String }, #[error("negative exponent is unsupported: {0}")] NegativeExponent(String), - #[error("nonlinear exponent is unsupported: {0}")] - NonlinearExponent(String), - #[error("variable base and exponent are unsupported: {0}")] - VariableBaseAndExponent(String), - #[error("factorial of a nonconstant expression is unsupported: {0}")] - FactorialOfNonconstant(String), + #[error("exponent is not a polynomial in variables and their logarithms: {0}")] + UnsupportedExponent(String), + #[error("variable base is not a polynomial in variables and their logarithms: {0}")] + UnsupportedVariableBase(String), + #[error("factorial argument is not a polynomial in variables and their logarithms: {0}")] + UnsupportedFactorial(String), #[error("invalid exponential base: {0}")] InvalidExponentialBase(String), #[error("represented exponent is outside the growth domain: {0}")] RepresentedExponentOutOfRange(String), - #[error("exponential factor {base}^({coefficient} * {variable}) decreases as {variable} grows")] - DecayingExponential { - base: String, - variable: String, - coefficient: String, - }, + #[error("decaying exponential has a widened exponent: {0}")] + DecayingWidenedExponent(String), + #[error("logarithm argument lacks a proved lower bound: {0}")] + UnprovedLogarithmBound(String), #[error("missing substitution for {0}")] MissingSubstitution(String), #[error("Big-O antichain has {terms} terms, exceeding the limit of {limit}")] @@ -352,90 +362,87 @@ pub enum GrowthFailure { } impl GrowthTerm { - /// The `O(1)` term (all maps empty). + /// The `O(1)` term. fn one() -> Self { GrowthTerm { - variables: BTreeMap::new(), + atoms: BTreeMap::new(), } } - fn insert(&mut self, variable: Box, growth: VariableGrowth) { - if !growth.is_empty() { - self.variables.insert(variable, growth); + fn from_atoms(atoms: impl IntoIterator) -> Self { + let mut term = GrowthTerm::one(); + for (atom, weight) in atoms { + term.add_weight(atom, &weight, &BigRational::one()); } + term } - /// Raise this term to a nonnegative real power `k` (scale every exponent). - /// Log powers are `u32`; a fractional result is rounded **up** (a valid - /// upper bound, since `(log v)^p ≤ (log v)^⌈p⌉` for `v ≥ 2`). - fn pow(&self, k: &BigRational) -> Result { - let mut r = GrowthTerm::one(); - for (variable, growth) in &self.variables { - r.insert(variable.clone(), growth.pow(k)?); + fn add_weight(&mut self, atom: Atom, weight: &Weight, scale: &BigRational) { + if matches!(&atom, Atom::Mono(monomial) if monomial.is_empty()) { + return; // constant factor + } + let entry = self.atoms.entry(atom.clone()).or_default(); + weight_add(entry, weight, scale); + if entry.len() > 1 && weight_sign(entry) == Sign::Zero { + entry.clear(); + } + if entry.is_empty() { + self.atoms.remove(&atom); } - Ok(r) } - /// Multiply two monomials (add matching exponents). - fn mul(&self, other: &GrowthTerm) -> Result { - let mut t = self.clone(); - for (variable, growth) in &other.variables { - let combined = match t.variables.get(variable) { - Some(current) => current.mul(growth)?, - None => growth.clone(), - }; - if combined.is_empty() { - t.variables.remove(variable); - } else { - t.variables.insert(variable.clone(), combined); - } + /// Raise this term to a nonnegative constant power (scale every weight). + fn pow(&self, k: &BigRational) -> GrowthTerm { + GrowthTerm::from_atoms( + self.atoms + .iter() + .map(|(atom, weight)| (atom.clone(), weight_scale(weight, k))), + ) + } + + /// Multiply two terms (add matching weights). + fn mul(&self, other: &GrowthTerm) -> GrowthTerm { + let mut term = self.clone(); + for (atom, weight) in &other.atoms { + term.add_weight(atom.clone(), weight, &BigRational::one()); } - Ok(t) + term } - /// Partial order on terms: `Some(Greater)` iff `self` dominates `other` - /// (`≥` on every variable and `>` on at least one). Per variable, - /// exponential products are compared only when a symbolic proof succeeds; - /// polynomial degree and log power then break proven exponential ties. - /// Returns `None` for incomparable or unproved terms. + /// Partial order on terms: `Some(Greater)` iff `self` provably grows + /// unboundedly faster than `other`, `Some(Equal)` iff their ratio is + /// provably bounded both ways. Returns `None` for incomparable or unproved + /// terms. fn cmp(&self, other: &GrowthTerm) -> Option { - let mut saw_gt = false; - let mut saw_lt = false; - let empty = VariableGrowth::empty(); - let mut left = self.variables.iter().peekable(); - let mut right = other.variables.iter().peekable(); - loop { - let (a, b) = match (left.peek(), right.peek()) { - (None, None) => break, - (Some((left_variable, _)), Some((right_variable, _))) => { - match left_variable.cmp(right_variable) { - Ordering::Less => (left.next().unwrap().1, &empty), - Ordering::Greater => (&empty, right.next().unwrap().1), - Ordering::Equal => (left.next().unwrap().1, right.next().unwrap().1), - } - } - (Some(_), None) => (left.next().unwrap().1, &empty), - (None, Some(_)) => (&empty, right.next().unwrap().1), - }; - let order = a.cmp_growth(b)?; - match order { - Ordering::Greater => saw_gt = true, - Ordering::Less => saw_lt = true, - Ordering::Equal => {} - } + if self == other { + return Some(Ordering::Equal); } - match (saw_gt, saw_lt) { - (true, true) => None, - (true, false) => Some(Ordering::Greater), - (false, true) => Some(Ordering::Less), - (false, false) => Some(Ordering::Equal), + let mut difference: BTreeMap<&Atom, Weight> = self + .atoms + .iter() + .map(|(atom, weight)| (atom, weight.clone())) + .collect(); + for (atom, weight) in &other.atoms { + weight_add( + difference.entry(atom).or_default(), + weight, + &-BigRational::one(), + ); + } + let signs: Vec<(&Atom, Sign)> = difference + .iter() + .map(|(atom, weight)| (*atom, weight_sign(weight))) + .filter(|(_, sign)| *sign != Sign::Zero) + .collect(); + if signs.is_empty() { + Some(Ordering::Equal) + } else if grows_unboundedly(&signs, Sign::Positive) { + Some(Ordering::Greater) + } else if grows_unboundedly(&signs, Sign::Negative) { + Some(Ordering::Less) + } else { + None } - } - - /// `true` iff `self` dominates `other` (grows at least as fast, and strictly - /// faster on at least one variable). - fn dominates(&self, other: &GrowthTerm) -> bool { - matches!(self.cmp(other), Some(Ordering::Greater)) } /// `true` iff `self` dominates `other` or is asymptotically equal to it. @@ -447,6 +454,19 @@ impl GrowthTerm { } } +/// The log difference tends to infinity in direction `up`: one atom with that +/// sign strictly dominates every atom whose sign is not `up`, which also bounds +/// their sum. +fn grows_unboundedly(signs: &[(&Atom, Sign)], up: Sign) -> bool { + signs.iter().any(|(high, sign)| { + *sign == up + && signs + .iter() + .filter(|(_, sign)| *sign != up) + .all(|(low, _)| atom_dominates(high, low)) + }) +} + impl Growth { pub(crate) fn unknown(failure: GrowthFailure) -> Self { Self(GrowthState::Unknown(vec![failure])) @@ -454,11 +474,20 @@ impl Growth { pub fn failures(&self) -> Option<&[GrowthFailure]> { match &self.0 { - GrowthState::Known(_) => None, + GrowthState::Known { .. } => None, GrowthState::Unknown(failures) => Some(failures), } } + /// Whether a known growth is tight (`Exact`, a Θ bound) or only an + /// `UpperBound`; `None` for unknown growth. + pub fn relation(&self) -> Option { + match &self.0 { + GrowthState::Known { relation, .. } => Some(*relation), + GrowthState::Unknown(_) => None, + } + } + /// Compute the growth class of an expression in a single bottom-up pass. pub fn from_expr(expr: &Expr) -> Growth { let analysis = AlgebraicAnalysis::new(&[expr]); @@ -476,7 +505,7 @@ impl Growth { /// some term of `self` — the standard antichain (Pareto) comparison. pub fn dominates(&self, other: &Growth) -> bool { match (&self.0, &other.0) { - (GrowthState::Known(a), GrowthState::Known(b)) => { + (GrowthState::Known { terms: a, .. }, GrowthState::Known { terms: b, .. }) => { b.iter().all(|tb| a.iter().any(|ta| ta.dominates_or_eq(tb))) } _ => false, @@ -492,23 +521,19 @@ impl Growth { pub fn to_expr(&self) -> Option { match &self.0 { GrowthState::Unknown(_) => None, - GrowthState::Known(terms) => { - if terms.is_empty() { - return Some(Expr::integer(1)); - } - let mut it = terms.iter().map(term_to_expr); - let mut acc = it.next().unwrap(); - for e in it { - acc = acc + e; - } - Some(acc) - } + GrowthState::Known { terms, .. } => Some( + terms + .iter() + .map(term_to_expr) + .reduce(|acc, e| acc + e) + .unwrap_or_else(|| Expr::integer(1)), + ), } } /// Canonical Big-O string for this growth class: `O()` for a bounded - /// class, or `O(?)` for unknown growth (no honest asymptotic bound — - /// nonlinear exponent or factorial). This is the single source of truth for + /// class, or `O(?)` for unknown growth (no honest asymptotic bound). + /// This is the single source of truth for /// how a growth is displayed as Big-O; presentation layers must call it rather /// than re-deriving the mapping (and the `Unknown` spelling) themselves. pub fn to_big_o(&self) -> String { @@ -531,7 +556,10 @@ fn growth_from_analysis( let facts = analysis.facts(expression); if facts.is_constant { let growth = if facts.constant_domain == Some(true) { - constant_growth() + known( + vec![GrowthTerm::one()], + exact_if(facts.sign == Some(Ordering::Greater)), + ) } else { unknown(GrowthFailure::InvalidConstantDomain { expression: expression.to_string(), @@ -543,14 +571,13 @@ fn growth_from_analysis( let growth = match expression.node() { ExprNode::Const(_) => unreachable!("constants are handled before node projection"), - ExprNode::Var(variable) => { - let mut term = GrowthTerm::one(); - term.insert( - variable.as_str().into(), - VariableGrowth::polynomial(BigRational::one()), - ); - exact_growth(vec![term]) - } + ExprNode::Var(variable) => known( + vec![GrowthTerm::from_atoms([( + ln_atom(variable.as_str()), + natural(BigRational::one()), + )])], + ParameterRelation::Exact, + ), ExprNode::Add(values) => values .iter() .map(|value| growth_from_analysis(value, analysis, memo)) @@ -578,78 +605,455 @@ fn growth_from_analysis( } else { pow_const(growth_from_analysis(base, analysis, memo), power) } - } else if let ExprNode::Exp(argument) = base.node() { - match analysis.facts(argument).exact_rational.as_ref() { - Some(coefficient) => exponential( - ExpBase::Natural, - scale_growth_linear(exponent_facts.linear.clone(), coefficient), - exponent, - ), - None => unknown(GrowthFailure::InvalidExponentialBase(base.to_string())), - } - } else if let Some(base_value) = base_facts.exact_rational.as_ref() { - if base_value.is_positive() { - exponential( - ExpBase::Rational(base_value.clone()), - growth_linear(exponent_facts.linear.clone()), - exponent, - ) - } else { - unknown(GrowthFailure::InvalidExponentialBase(base.to_string())) - } - } else if base_facts.is_constant { - unknown(GrowthFailure::InvalidExponentialBase(base.to_string())) } else { - unknown(GrowthFailure::VariableBaseAndExponent( - expression.to_string(), - )) + power_growth(base, exponent, analysis) } } - ExprNode::Exp(value) => { - let value_growth = growth_from_analysis(value, analysis, memo); - if value_growth.failures().is_some() { - value_growth - } else { - exponential( - ExpBase::Natural, - growth_linear(analysis.facts(value).linear.clone()), - expression, - ) - } - } - ExprNode::Log(value) => log_growth(growth_from_analysis(value, analysis, memo)), - ExprNode::Factorial(_) => unknown(GrowthFailure::FactorialOfNonconstant( - expression.to_string(), - )), + ExprNode::Exp(value) => match exponent_poly(value, analysis) { + Ok(exponent) => exponential( + LogBound::constant(natural(BigRational::one())), + &exponent, + value, + ), + Err(failure) => unknown(failure), + }, + ExprNode::Log(value) => log_growth(growth_from_analysis(value, analysis, memo), value), + ExprNode::Factorial(value) => factorial_growth(expression, value, analysis, memo), }; memo.insert(expression.node_identity(), growth.clone()); growth } -fn growth_linear( - linear: Option>, -) -> Option, BigRational>> { - Some( - linear? +/// `x!` is `x^(x+1/2)·e^(−x)` up to a constant (Stirling) when `x^x` is exact, +/// which holds when `x` expands exactly to one monomial; otherwise `x! ≤ x^x`. +fn factorial_growth( + expression: &Expr, + value: &Expr, + analysis: &AlgebraicAnalysis, + memo: &mut HashMap, +) -> Growth { + let self_power = power_growth(value, value, analysis); + match self_power.relation() { + None => unknown(GrowthFailure::UnsupportedFactorial(expression.to_string())), + Some(ParameterRelation::UpperBound) => self_power, + Some(ParameterRelation::Exact) => { + let root = pow_const( + growth_from_analysis(value, analysis, memo), + &BigRational::new(BigInt::from(1), BigInt::from(2)), + ); + let decay = match exponent_poly(value, analysis) { + Ok(poly) => exponential( + LogBound::constant(natural(-BigRational::one())), + &poly, + value, + ), + Err(failure) => unknown(failure), + }; + mul(mul(self_power, root), decay) + } + } +} + +/// `base^exponent` for a nonconstant exponent. +fn power_growth(base: &Expr, exponent: &Expr, analysis: &AlgebraicAnalysis) -> Growth { + match (log_bound(base, analysis), exponent_poly(exponent, analysis)) { + (Ok(bound), Ok(poly)) => exponential(bound, &poly, exponent), + (Err(failure), _) | (_, Err(failure)) => unknown(failure), + } +} + +/// Bound `ln base` for an exponential base: `e^c`, a positive rational, or a +/// polynomial in variables and their logarithms. +fn log_bound(base: &Expr, analysis: &AlgebraicAnalysis) -> Result { + let facts = analysis.facts(base); + let invalid = || GrowthFailure::InvalidExponentialBase(base.to_string()); + if let ExprNode::Exp(argument) = base.node() { + let coefficient = analysis.facts(argument).exact_rational.clone(); + return coefficient + .map(|c| LogBound::constant(natural(c))) + .ok_or_else(invalid); + } + if facts.is_constant { + return match facts.exact_rational.as_ref() { + Some(value) if value.is_positive() => Ok(LogBound::constant(BTreeMap::from([( + ExpBase::Rational(value.clone()), + BigRational::one(), + )]))), + _ => Err(invalid()), + }; + } + exponent_poly(base, analysis) + .ok() + .and_then(|poly| poly.monomial_bound()) + .map(LogBound::from_monomial_bound) + .ok_or_else(|| GrowthFailure::UnsupportedVariableBase(base.to_string())) +} + +/// An upper bound `ln B ≤ Σ weight · monomial` for an exponential base `B`. +struct LogBound { + terms: BTreeMap, + /// The bound equals `ln B` exactly. + exact: bool, +} + +impl LogBound { + fn constant(weight: Weight) -> Self { + LogBound { + terms: BTreeMap::from([(Monomial::new(), weight)]), + exact: true, + } + } + + /// `ln B ≤ ln C + Σ_v (power_v + log_v) · ln v` from [`MonomialBound`]. + fn from_monomial_bound(bound: MonomialBound) -> Self { + let mut terms: BTreeMap = bound + .largest .into_iter() - .map(|(symbol, coefficient)| (symbol.as_str().into(), coefficient)) - .collect(), - ) + .map(|(variable, (power, log_power))| { + ( + single_variable(&variable, BigRational::zero(), BigRational::one()), + natural(power + log_power), + ) + }) + .collect(); + if !bound.scale.is_one() { + terms.insert( + Monomial::new(), + BTreeMap::from([(ExpBase::Rational(bound.scale), BigRational::one())]), + ); + } + LogBound { + terms, + exact: bound.exact, + } + } +} + +/// `Σ c_i μ_i ≤ scale · largest` over the positive part, where `largest` takes +/// each variable's largest power and log power (`ln ln v` widens to `ln v`). +struct MonomialBound { + scale: BigRational, + largest: Monomial, + /// The polynomial is exactly one monomial without log powers, so + /// `B = scale · largest`. + exact: bool, } -fn scale_growth_linear( - linear: Option>, - coefficient: &BigRational, -) -> Option, BigRational>> { - Some( - linear? +/// `exp(E · ln B)` as a single growth term. +fn exponential(bound: LogBound, exponent: &ExponentPoly, source: &Expr) -> Growth { + let decaying = bound + .terms + .values() + .any(|weight| weight_sign(weight) != Sign::Positive); + if !exponent.exact && decaying { + return unknown(GrowthFailure::DecayingWidenedExponent(source.to_string())); + } + // Multiplying two signed upper bounds, or reversing one with a negative + // exponent, is unsound when the actual base can approach zero. + if !bound.exact && (!exponent.exact || exponent.has_negative()) { + return unknown(GrowthFailure::UnsupportedExponent(source.to_string())); + } + let exponent_terms = exponent.terms.iter(); + let mut term = GrowthTerm::one(); + for (monomial, coefficient) in exponent_terms { + for (log_monomial, weight) in &bound.terms { + term.add_weight( + Atom::Mono(monomial_mul(monomial, log_monomial)), + weight, + coefficient, + ); + } + } + known(vec![term], exact_if(bound.exact && exponent.exact)) +} + +/// An exponent expanded as `Σ coefficient · monomial`. +#[derive(Clone, Debug)] +struct ExponentPoly { + terms: BTreeMap, + /// The expansion equals the exponent; otherwise it is an upper bound that + /// may only be combined monotonically. + exact: bool, +} + +impl ExponentPoly { + fn constant(value: BigRational) -> Self { + Self::from_terms([(Monomial::new(), value)], true) + } + + fn from_terms(terms: impl IntoIterator, exact: bool) -> Self { + let mut result = ExponentPoly { + terms: BTreeMap::new(), + exact, + }; + for (monomial, coefficient) in terms { + *result + .terms + .entry(monomial) + .or_insert_with(BigRational::zero) += coefficient; + } + result.terms.retain(|_, coefficient| !coefficient.is_zero()); + result + } + + fn has_negative(&self) -> bool { + self.terms.values().any(Signed::is_negative) + } + + fn checked(self, source: &Expr) -> Result { + if self.terms.len() > EXPONENT_TERM_CAP { + return Err(GrowthFailure::RepresentedExponentOutOfRange( + source.to_string(), + )); + } + Ok(self) + } + + fn add(&self, other: &Self) -> Self { + Self::from_terms( + self.terms + .iter() + .chain(&other.terms) + .map(|(m, c)| (m.clone(), c.clone())), + self.exact && other.exact, + ) + } + + fn mul(&self, other: &Self, source: &Expr) -> Result { + // An upper bound times a factor that may be negative is not an upper bound. + if (!self.exact && !other.exact) + || (!self.exact && other.has_negative()) + || (!other.exact && self.has_negative()) + { + return Err(GrowthFailure::UnsupportedExponent(source.to_string())); + } + let mut terms = Vec::with_capacity(self.terms.len() * other.terms.len()); + for (a, x) in &self.terms { + for (b, y) in &other.terms { + terms.push((monomial_mul(a, b), x * y)); + } + } + Self::from_terms(terms, self.exact && other.exact).checked(source) + } + + fn pow(&self, k: &BigRational, source: &Expr) -> Result { + if k.is_zero() { + return Ok(Self::constant(BigRational::one())); + } + if k.is_one() { + return Ok(self.clone()); + } + // A signed upper bound may hide a large negative value. Its square or + // other nonlinear power cannot be bounded by powering that upper bound. + if !self.exact { + return Err(GrowthFailure::UnsupportedExponent(source.to_string())); + } + let out_of_range = || GrowthFailure::RepresentedExponentOutOfRange(source.to_string()); + if let (1, Some((monomial, coefficient))) = (self.terms.len(), self.terms.iter().next()) { + let exact_coefficient = if k.is_integer() { + Some(bounded_pow(coefficient, &k.to_integer()).ok_or_else(out_of_range)?) + } else { + coefficient.is_one().then(BigRational::one) + }; + return match exact_coefficient { + Some(value) => Ok(Self::from_terms( + [(monomial_pow(monomial, k), value)], + self.exact, + )), + None if coefficient.is_positive() => Ok(Self::from_terms( + [( + monomial_pow(monomial, k), + ceil_power_bound(coefficient, k).ok_or_else(out_of_range)?, + )], + false, + )), + None => Err(GrowthFailure::UnsupportedExponent(source.to_string())), + }; + } + if let Some(power) = k + .is_integer() + .then(|| k.to_integer().to_u32()) + .flatten() + .filter(|power| *power <= EXPONENT_POWER_CAP) + { + let mut result = Self::constant(BigRational::one()); + for _ in 0..power { + result = result.mul(self, source)?; + } + return Ok(result); + } + // |Σ x_i|^k ≤ m^max(k-1, 0) · Σ |x_i|^k. Negative terms cannot + // be discarded before an even power, including beyond the expansion cap. + let magnitudes: Vec<_> = self + .terms + .iter() + .map(|(monomial, coefficient)| (monomial, coefficient.abs())) + .collect(); + let excess = (k - BigRational::one()).max(BigRational::zero()).ceil(); + let spread = ceil_power_bound( + &BigRational::from_integer(BigInt::from(magnitudes.len())), + &excess, + ) + .ok_or_else(out_of_range)?; + let mut terms = Vec::with_capacity(magnitudes.len()); + for (monomial, coefficient) in magnitudes { + let bound = ceil_power_bound(&coefficient, k).ok_or_else(out_of_range)?; + terms.push((monomial_pow(monomial, k), &spread * bound)); + } + Self::from_terms(terms, false).checked(source) + } + + fn monomial_bound(&self) -> Option { + let positive: Vec<(&Monomial, &BigRational)> = self + .terms + .iter() + .filter(|(_, coefficient)| coefficient.is_positive()) + .collect(); + if positive.is_empty() { + return None; + } + let scale: BigRational = positive.iter().map(|(_, coefficient)| *coefficient).sum(); + let mut largest = Monomial::new(); + for (monomial, _) in &positive { + for (variable, (power, log_power)) in *monomial { + let entry = largest + .entry(variable.clone()) + .or_insert_with(|| (BigRational::zero(), BigRational::zero())); + entry.0 = entry.0.clone().max(power.clone()); + entry.1 = entry.1.clone().max(log_power.clone()); + } + } + let exact = self.exact + && self.terms.len() == 1 + && largest.values().all(|(_, log_power)| log_power.is_zero()); + Some(MonomialBound { + scale, + largest, + exact, + }) + } + + /// `ln(Σ c_i μ_i) ≤ ln C + Σ_v (power_v + log_v) · ln v` over the positive part. + fn log(&self, source: &Expr) -> Result { + // Without a positive-coefficient exact argument, a signed logarithmic + // bound may hide arbitrary decay and cannot safely enter exponent algebra. + if !self.exact || self.has_negative() { + return Err(GrowthFailure::UnprovedLogarithmBound(source.to_string())); + } + let bound = self + .monomial_bound() + .ok_or_else(|| GrowthFailure::UnsupportedExponent(source.to_string()))?; + let logs = bound + .largest .into_iter() - .map(|(symbol, value)| (symbol.as_str().into(), coefficient * value)) - .collect(), - ) + .map(|(variable, (power, log_power))| { + ( + single_variable(&variable, BigRational::zero(), BigRational::one()), + power + log_power, + ) + }); + let constant = (Monomial::new(), ln_upper_bound(&bound.scale)); + let exact = bound.exact && bound.scale.is_one(); + Ok(Self::from_terms(logs.chain([constant]), exact)) + } +} + +/// A rational upper bound on `c^k` for `c > 0`, `k ≥ 0`, within the bit budget. +fn ceil_power_bound(c: &BigRational, k: &BigRational) -> Option { + if c <= &BigRational::one() { + return Some(BigRational::one()); + } + bounded_pow(c, &k.ceil().to_integer()) } -fn constant_growth() -> Growth { - exact_growth(vec![GrowthTerm::one()]) + +/// `c^power` for a nonnegative integer power, or `None` past the bit budget. +fn bounded_pow(c: &BigRational, power: &BigInt) -> Option { + if c.is_zero() || c.abs().is_one() { + let odd = (power % BigInt::from(2)).is_one(); + return Some(if c.is_negative() && odd { + -BigRational::one() + } else { + c.abs() + }); + } + let power = power.to_i32()?; + let bits = + (c.numer().bits() + c.denom().bits()).checked_mul(u64::from(power.unsigned_abs()))?; + (bits <= RATIONAL_POWER_BIT_CAP).then(|| c.pow(power)) +} + +/// A rational upper bound on `ln c` for `c > 0`. +fn ln_upper_bound(c: &BigRational) -> BigRational { + if c <= &BigRational::one() { + return BigRational::zero(); + } + let scale = 1_000_000_i64; + match c.to_f64().map(f64::ln).filter(|value| value.is_finite()) { + Some(value) => BigRational::new( + BigInt::from((value * scale as f64).floor() as i64 + 2), + BigInt::from(scale), + ), + // ln c ≤ bits(c) · ln 2 < bits(c) · 7/10. + None => BigRational::new(BigInt::from(c.numer().bits() * 7), BigInt::from(10)), + } +} + +/// Expand an exponent into a polynomial over variables and their logarithms. +fn exponent_poly( + expression: &Expr, + analysis: &AlgebraicAnalysis, +) -> Result { + let facts = analysis.facts(expression); + if facts.is_constant { + return facts + .exact_rational + .clone() + .map(ExponentPoly::constant) + .ok_or_else(|| GrowthFailure::UnsupportedExponent(expression.to_string())); + } + match expression.node() { + ExprNode::Var(variable) => Ok(ExponentPoly::from_terms( + [( + single_variable(variable.as_str(), BigRational::one(), BigRational::zero()), + BigRational::one(), + )], + true, + )), + ExprNode::Add(values) => { + let mut result = ExponentPoly::constant(BigRational::zero()); + for value in values.iter() { + result = result + .add(&exponent_poly(value, analysis)?) + .checked(expression)?; + } + Ok(result) + } + ExprNode::Mul(values) => { + let mut result = ExponentPoly::constant(BigRational::one()); + for value in values.iter() { + result = result.mul(&exponent_poly(value, analysis)?, expression)?; + } + Ok(result) + } + ExprNode::Pow(base, power) => match analysis.facts(power).exact_rational.as_ref() { + Some(k) if k.is_negative() => Err(GrowthFailure::NegativeExponent(power.to_string())), + Some(k) => exponent_poly(base, analysis)?.pow(k, expression), + None => Err(GrowthFailure::UnsupportedExponent(expression.to_string())), + }, + ExprNode::Log(value) => exponent_poly(value, analysis)?.log(expression), + ExprNode::Exp(_) | ExprNode::Factorial(_) => { + Err(GrowthFailure::UnsupportedExponent(expression.to_string())) + } + ExprNode::Const(_) => unreachable!("constants are handled above"), + } +} + +fn exact_if(exact: bool) -> ParameterRelation { + if exact { + ParameterRelation::Exact + } else { + ParameterRelation::UpperBound + } } fn unknown(failure: GrowthFailure) -> Growth { @@ -669,59 +1073,64 @@ fn merge_unknown(left: Growth, right: Growth) -> Growth { Growth(GrowthState::Unknown(failures)) } -/// Render one monomial as a product of its factors (or `Const(1)` when empty). +/// Render one term as a product of its factors (`Const(1)` when empty). +/// `Expr` products sort their factors canonically. fn term_to_expr(t: &GrowthTerm) -> Expr { - let mut factors: Vec = Vec::new(); - for (variable, growth) in &t.variables { - factors.extend( - growth - .exp + t.atoms + .iter() + .flat_map(|(atom, weight)| { + weight .iter() - .map(|(base, coefficient)| exp_factor(variable, base, coefficient)), - ); - if !growth.poly.is_zero() { - factors.push(poly_factor(variable, &growth.poly)); - } - if growth.log != 0 { - factors.push(log_factor(variable, growth.log)); - } - } - let mut it = factors.into_iter(); - match it.next() { - None => Expr::integer(1), - Some(first) => it.fold(first, |acc, f| acc * f), - } -} - -/// Render a stored exponential factor without changing its base or coefficient. -fn exp_factor(v: &str, base: &ExpBase, coefficient: &BigRational) -> Expr { - let exponent = if coefficient.is_one() { - Expr::variable(v) - } else { - Expr::constant(coefficient.clone()) * Expr::variable(v) - }; - match base { - ExpBase::Rational(base) => Expr::pow(Expr::constant(base.clone()), exponent), - ExpBase::Natural => Expr::exp(exponent), - } + .map(move |(base, coefficient)| atom_factor(atom, base, coefficient)) + }) + .fold(Expr::integer(1), |acc, f| acc * f) } -/// Render `v^degree` (`Display` turns degree `0.5` into `sqrt(v)`). -fn poly_factor(v: &str, degree: &BigRational) -> Expr { - if degree.is_one() { - Expr::variable(v) - } else { - Expr::pow(Expr::variable(v), Expr::constant(degree.clone())) - } +fn monomial_expr(monomial: &Monomial) -> Expr { + monomial + .iter() + .flat_map(|(variable, (power, log_power))| { + let variable = Expr::variable(variable.as_ref()); + [ + Expr::pow(variable.clone(), Expr::constant(power.clone())), + Expr::pow(Expr::log(variable), Expr::constant(log_power.clone())), + ] + }) + .fold(Expr::integer(1), |acc, f| acc * f) } -/// Render `(log v)^power`. -fn log_factor(v: &str, power: u32) -> Expr { - let log = Expr::log(Expr::variable(v)); - if power == 1 { - log - } else { - Expr::pow(log, Expr::integer(power)) +/// Render `exp(coefficient · ln(base) · atom)` with its authoritative base. +fn atom_factor(atom: &Atom, base: &ExpBase, coefficient: &BigRational) -> Expr { + let c = Expr::constant(coefficient.clone()); + match (atom, base) { + (Atom::LogLog(variable), ExpBase::Natural) => { + Expr::pow(Expr::log(Expr::variable(variable.as_ref())), c) + } + (Atom::LogLog(variable), ExpBase::Rational(value)) => Expr::pow( + Expr::constant(value.clone()), + c * Expr::log(Expr::log(Expr::variable(variable.as_ref()))), + ), + (Atom::Mono(monomial), ExpBase::Natural) => { + let log_variable = monomial + .iter() + .find(|(_, (_, log_power))| log_power >= &BigRational::one()) + .map(|(variable, _)| variable.clone()); + match log_variable { + // exp(c · R · ln w) = w^(c · R). + Some(variable) => { + let mut rest = monomial.clone(); + let entry = rest + .get_mut(&variable) + .expect("variable is in the monomial"); + entry.1 -= BigRational::one(); + Expr::pow(Expr::variable(variable.as_ref()), c * monomial_expr(&rest)) + } + None => Expr::exp(c * monomial_expr(monomial)), + } + } + (Atom::Mono(monomial), ExpBase::Rational(value)) => { + Expr::pow(Expr::constant(value.clone()), c * monomial_expr(monomial)) + } } } @@ -735,20 +1144,29 @@ fn prune(mut terms: Vec) -> Vec { terms.sort(); let mut result: Vec = Vec::new(); for t in terms { - if result.iter().any(|r| r.dominates_or_eq(&t)) { - continue; + let mut kept = true; + result.retain(|r| match r.cmp(&t) { + Some(Ordering::Greater | Ordering::Equal) => { + kept = false; + true + } + Some(Ordering::Less) => false, + None => true, + }); + if kept { + result.push(t); } - result.retain(|r| !t.dominates(r)); - result.push(t); } result } -fn exact_growth(terms: Vec) -> Growth { - finish_growth(terms) +fn both_exact(a: &Growth, b: &Growth) -> bool { + [a, b] + .iter() + .all(|growth| growth.relation() == Some(ParameterRelation::Exact)) } -fn finish_growth(terms: Vec) -> Growth { +fn known(terms: Vec, relation: ParameterRelation) -> Growth { let terms = prune(terms); if terms.len() > ANTICHAIN_CAP { return unknown(GrowthFailure::AntichainLimitExceeded { @@ -756,17 +1174,28 @@ fn finish_growth(terms: Vec) -> Growth { terms: terms.len(), }); } - Growth(GrowthState::Known(terms)) + Growth(GrowthState::Known { terms, relation }) } -/// Antichain union (asymptotic `+ ≍ max`). +/// Antichain union (asymptotic `+ ≍ max`). The sum stays exact when every +/// widened summand is strictly dominated by an exact one: growths bound +/// magnitudes, so a widened summand, even a negative one, is `o` of the exact part. fn add(a: Growth, b: Growth) -> Growth { if a.failures().is_some() || b.failures().is_some() { return merge_unknown(a, b); } - let mut terms = into_terms(a); - terms.extend(into_terms(b)); - finish_growth(terms) + let covers = |high: &Growth, low: &Growth| { + high.relation() == Some(ParameterRelation::Exact) + && terms(low).iter().all(|l| { + terms(high) + .iter() + .any(|h| h.cmp(l) == Some(Ordering::Greater)) + }) + }; + let relation = exact_if(both_exact(&a, &b) || covers(&a, &b) || covers(&b, &a)); + let mut sum = into_terms(a); + sum.extend(into_terms(b)); + known(sum, relation) } /// Pairwise product of two antichains. @@ -774,115 +1203,100 @@ fn mul(a: Growth, b: Growth) -> Growth { if a.failures().is_some() || b.failures().is_some() { return merge_unknown(a, b); } - let x = into_terms(a); - let y = into_terms(b); - let mut product = Vec::with_capacity(x.len() * y.len()); - for tx in &x { - for ty in &y { - match tx.mul(ty) { - Ok(term) => product.push(term), - Err(failure) => return unknown(failure), - } + let relation = exact_if(both_exact(&a, &b)); + let mut product = Vec::with_capacity(terms(&a).len() * terms(&b).len()); + for tx in terms(&a) { + for ty in terms(&b) { + product.push(tx.mul(ty)); } } - finish_growth(product) + known(product, relation) } /// Raise a whole antichain to a nonnegative real power `k` (raise each term). fn pow_const(g: Growth, k: &BigRational) -> Growth { match g.0 { GrowthState::Unknown(failures) => Growth(GrowthState::Unknown(failures)), - GrowthState::Known(terms) => match terms.iter().map(|term| term.pow(k)).collect() { - Ok(terms) => finish_growth(terms), - Err(failure) => unknown(failure), - }, + GrowthState::Known { terms, relation } => { + known(terms.iter().map(|term| term.pow(k)).collect(), relation) + } } } -fn into_terms(growth: Growth) -> Vec { - match growth.0 { - GrowthState::Known(terms) => terms, +fn terms(growth: &Growth) -> &[GrowthTerm] { + match &growth.0 { + GrowthState::Known { terms, .. } => terms, GrowthState::Unknown(_) => unreachable!("unknown growth is handled before term access"), } } -/// Transfer function for a symbolic fixed-base exponential. -fn exponential( - base: ExpBase, - linear: Option, BigRational>>, - exponent: &Expr, -) -> Growth { - let direction = base.direction(); - if direction == Ordering::Equal { - // 1^x = 1 for every x: bounded by O(1). - return exact_growth(vec![GrowthTerm::one()]); - } - match linear { - None => unknown(GrowthFailure::NonlinearExponent(exponent.to_string())), - Some(coeffs) => { - let mut term = GrowthTerm::one(); - for (v, coeff) in coeffs { - if (direction == Ordering::Greater && coeff.is_positive()) - || (direction == Ordering::Less && coeff.is_negative()) - { - term.insert(v, VariableGrowth::exponential(base.clone(), coeff)); - } else if !coeff.is_zero() { - return unknown(GrowthFailure::DecayingExponential { - base: match &base { - ExpBase::Rational(value) => value.to_string(), - ExpBase::Natural => "e".to_string(), - }, - variable: v.to_string(), - coefficient: coeff.to_string(), - }); - } - } - exact_growth(vec![term]) - } +fn into_terms(growth: Growth) -> Vec { + match growth.0 { + GrowthState::Known { terms, .. } => terms, + GrowthState::Unknown(_) => unreachable!("unknown growth is handled before term access"), } } /// Transfer function for `Log(a)`: `log` of an antichain is `log` of its /// dominant term(s), unioned. Uses `log(n^a · m^b) ≍ log n + log m` and -/// `log(2^(r·n)) ≍ n`. -fn log_growth(g: Growth) -> Growth { +/// `log(2^(r·n)) ≍ n`. Since `ln max = max ln`, the result is exact when the +/// argument is exact and every term's logarithm is exact. An upper bound alone +/// cannot bound the magnitude of a logarithm: its positive argument may approach +/// zero arbitrarily faster than the represented upper bound. +fn log_growth(g: Growth, argument: &Expr) -> Growth { match g.0 { GrowthState::Unknown(failures) => Growth(GrowthState::Unknown(failures)), - GrowthState::Known(terms) => { - let mut out = Vec::new(); - for t in &terms { - out.extend(log_term(t)); - } - if out.is_empty() { - out.push(GrowthTerm::one()); // log(O(1)) = O(1) - } - finish_growth(out) - } + GrowthState::Known { + relation: ParameterRelation::UpperBound, + .. + } => unknown(GrowthFailure::UnprovedLogarithmBound(argument.to_string())), + GrowthState::Known { terms, .. } => known( + terms.iter().flat_map(log_term).collect(), + exact_if(terms.iter().all(log_is_exact)), + ), } } -/// `log` of a single monomial, returned as its own (small) antichain of -/// summands. `log(∏ baseᵢ^(rᵢ·vᵢ) · ∏vⱼ^aⱼ · ∏(log vₖ)^sₖ)` distributes over the -/// product into a *sum* of the log of each factor, so every factor class of the -/// monomial contributes its own summand — none may be dropped (e.g. `log(2^n·m)` -/// is `n + log m`, not `n`). `finish_growth`/`prune` then collapse any dominated -/// summands (so `log(2^n·n^2)` reduces back to `n`). +/// `ln t = Σ w · atom` is `Θ` of its positive atoms when one of them strictly +/// dominates every atom not proved positive (so `ln t` is unbounded), and no +/// positive `ln ln v` atom was widened to `ln v`. +fn log_is_exact(t: &GrowthTerm) -> bool { + let signs: Vec<(&Atom, Sign)> = t + .atoms + .iter() + .map(|(atom, weight)| (atom, weight_sign(weight))) + .collect(); + grows_unboundedly(&signs, Sign::Positive) + && !signs + .iter() + .any(|(atom, sign)| matches!(atom, Atom::LogLog(_)) && *sign == Sign::Positive) +} + +/// `log` of a single term, returned as its own (small) antichain of summands: +/// `|ln exp(Σ w · atom)| ≤ Σ |w| · atom`, so every atom contributes the growth +/// of the atom itself, whatever its sign: a decaying term has a large negative +/// logarithm (`log(exp(-n)) = -n`), and growths bound magnitudes so that sums +/// and powers stay sound. Dominated summands are pruned. None may be dropped (e.g. +/// `log(2^n·m)` is `n + log m`, not `n`); `prune` then collapses +/// any dominated summands (so `log(2^n·n^2)` reduces back to `n`). `ln ln ln v` +/// widens to `ln ln v`. fn log_term(t: &GrowthTerm) -> Vec { + let one = BigRational::one(); let mut out = Vec::new(); - for (variable, growth) in &t.variables { - if !growth.exp.is_empty() { - let mut term = GrowthTerm::one(); - term.insert( - variable.clone(), - VariableGrowth::polynomial(BigRational::one()), - ); - out.push(term); - } - if growth.poly.is_positive() || growth.log != 0 { - let mut term = GrowthTerm::one(); - term.insert(variable.clone(), VariableGrowth::logarithmic(1)); - out.push(term); - } + for atom in t.atoms.keys() { + out.push(match atom { + Atom::Mono(monomial) => GrowthTerm::from_atoms(monomial.iter().flat_map( + |(variable, (power, log_power))| { + [ + (ln_atom(variable), natural(power.clone())), + (Atom::LogLog(variable.clone()), natural(log_power.clone())), + ] + }, + )), + Atom::LogLog(variable) => { + GrowthTerm::from_atoms([(Atom::LogLog(variable.clone()), natural(one.clone()))]) + } + }); } // Empty term: log(O(1)) = O(1). if out.is_empty() { diff --git a/src/unit_tests/growth.rs b/src/unit_tests/growth.rs index d47f33d98..d4186a86e 100644 --- a/src/unit_tests/growth.rs +++ b/src/unit_tests/growth.rs @@ -1,55 +1,59 @@ //! Unit tests for the symbolic growth domain (`src/growth.rs`). use super::{ - ANTICHAIN_CAP, ExpBase, Growth, GrowthFailure, GrowthState, GrowthTerm, VariableGrowth, add, - exact_growth, mul, + ANTICHAIN_CAP, Atom, ExpBase, Growth, GrowthFailure, GrowthState, GrowthTerm, add, known, + ln_atom, mul, natural, single_variable, }; use crate::expr::{ AlgebraicAnalysis, Expr, ExprNode, evaluate_approximate, expression_from_approximation, }; +use crate::parameters::ParameterRelation; use crate::registry::variant_entries; use num_rational::BigRational; use num_traits::{FromPrimitive, One, Signed, Zero}; use serde::Deserialize; use std::cmp::Ordering; +use std::collections::BTreeMap; fn rat(value: f64) -> BigRational { BigRational::from_f64(value).unwrap() } -/// Build a term from `(exp, poly, logs)` entry lists. +/// Build a term from `(2^rate·v, v^degree, log(v)^power)` entry lists. fn term(exp: &[(&str, f64)], poly: &[(&str, f64)], logs: &[(&str, u32)]) -> GrowthTerm { - let mut result = GrowthTerm::one(); - for (variable, rate) in exp { - result.insert( - (*variable).into(), - VariableGrowth::exponential(ExpBase::Rational(rat(2.0)), rat(*rate)), - ); - } - for (variable, degree) in poly { - let growth = result - .variables - .entry((*variable).into()) - .or_insert_with(VariableGrowth::empty); - growth.poly = rat(*degree); - } - for (variable, power) in logs { - let growth = result - .variables - .entry((*variable).into()) - .or_insert_with(VariableGrowth::empty); - growth.log = *power; - } - result + let exp = exp.iter().map(|(variable, rate)| { + ( + Atom::Mono(single_variable( + variable, + BigRational::one(), + BigRational::zero(), + )), + BTreeMap::from([(ExpBase::Rational(rat(2.0)), rat(*rate))]), + ) + }); + let poly = poly + .iter() + .map(|(variable, degree)| (ln_atom(variable), natural(rat(*degree)))); + let logs = logs.iter().map(|(variable, power)| { + ( + Atom::LogLog((*variable).into()), + natural(BigRational::from_integer((*power).into())), + ) + }); + GrowthTerm::from_atoms(exp.chain(poly).chain(logs)) } fn terms_of(g: &Growth) -> &[GrowthTerm] { match &g.0 { - GrowthState::Known(terms) => terms, + GrowthState::Known { terms, .. } => terms, GrowthState::Unknown(failures) => panic!("expected known growth, got {failures:?}"), } } +fn exact(terms: Vec) -> Growth { + known(terms, ParameterRelation::Exact) +} + fn g(s: &str) -> Growth { Growth::from_expr(&Expr::parse(s)) } @@ -82,7 +86,7 @@ fn test_growth_relations_against_sympy_limits() { "../../problemreductions-expr/tests/fixtures/sympy_oracle.json" )) .unwrap(); - assert_eq!(fixture.growth_cases.len(), 14); + assert_eq!(fixture.growth_cases.len(), 29); for case in fixture.growth_cases { let left = g(&case.left); @@ -101,17 +105,6 @@ fn test_growth_relations_against_sympy_limits() { } } -fn exponential_growth(factors: &[(f64, f64)]) -> VariableGrowth { - VariableGrowth { - exp: factors - .iter() - .map(|(base, coefficient)| (ExpBase::Rational(rat(*base)), rat(*coefficient))) - .collect(), - poly: BigRational::zero(), - log: 0, - } -} - // --- Core verification cases --- /// 1. No-expansion regression: the nested sum-of-squares shape that OOM'd in @@ -214,9 +207,20 @@ fn test_every_registered_complexity_uses_the_shared_analysis() { let expression = Expr::parse(entry.complexity); let growth = Growth::from_expr(&expression); if let Some(rendered) = growth.to_expr() { + let reparsed = Growth::from_expr(&rendered); + // A rendered upper bound may legitimately re-parse as exact, but an + // exact normal form must stay exact. + if growth.relation() == Some(ParameterRelation::Exact) { + assert_eq!( + reparsed.relation(), + growth.relation(), + "{}", + entry.complexity + ); + } assert_eq!( - Growth::from_expr(&rendered), - growth, + terms_of(&reparsed), + terms_of(&growth), "{}: {}", entry.name, entry.complexity @@ -225,93 +229,69 @@ fn test_every_registered_complexity_uses_the_shared_analysis() { } } -/// Multi-base products remain incomparable when the conservative symbolic -/// rules cannot prove an ordering, even when a stronger algebra system could. +/// Mixed-base products are compared through their log weights: +/// `2^(2n)·3^n = 12^n` dominates `2^n·4^n = 8^n`. #[test] -fn test_growth_unproved_multi_base_comparison_is_retained() { - let left = g("2^(2*n) * 3^n"); - let right = g("2^n * 4^n"); +fn test_growth_mixed_base_products_are_compared() { + assert!(g("2^(2*n) * 3^n").dominates(&g("2^n * 4^n"))); + assert!(!g("2^n * 4^n").dominates(&g("2^(2*n) * 3^n"))); + assert_eq!(terms_of(&g("2^(2*n) * 3^n + 2^n * 4^n")).len(), 1); +} + +/// A rational base inside the certified interval for `e` stays unproved, +/// while a base just outside that interval has a provable order. +#[test] +fn test_growth_unproved_near_e_comparison_is_retained() { + let left = g("exp(n)"); + let right = g("2.71828182846^n"); assert!(!left.dominates(&right)); assert!(!right.dominates(&left)); - assert_eq!(terms_of(&g("2^(2*n) * 3^n + 2^n * 4^n")).len(), 2); + assert_eq!(terms_of(&g("exp(n) + 2.71828182846^n")).len(), 2); + + let above_interval = g("2.7182818285^n"); + assert!(above_interval.dominates(&left)); + assert!(!left.dominates(&above_interval)); + assert_eq!(terms_of(&g("exp(n) + 2.7182818285^n")).len(), 1); } #[test] fn test_exponential_product_proof_rules() { - let empty = VariableGrowth::empty(); - let two = exponential_growth(&[(2.0, 1.0)]); - let two_squared = exponential_growth(&[(2.0, 2.0)]); - let three = exponential_growth(&[(3.0, 1.0)]); - - assert_eq!(empty.cmp_exp(&empty), Some(Ordering::Equal)); - assert_eq!(empty.cmp_exp(&two), Some(Ordering::Less)); - assert_eq!(two.cmp_exp(&empty), Some(Ordering::Greater)); - assert_eq!(two_squared.cmp_exp(&two), Some(Ordering::Greater)); - assert_eq!(two.cmp_exp(&two_squared), Some(Ordering::Less)); - assert_eq!(three.cmp_exp(&two), Some(Ordering::Greater)); - assert_eq!(two.cmp_exp(&three), Some(Ordering::Less)); - - assert_eq!( - exponential_growth(&[(3.0, 2.0)]).cmp_exp(&exponential_growth(&[(2.0, 1.0)])), - Some(Ordering::Greater) - ); - assert_eq!( - exponential_growth(&[(2.0, 1.0)]).cmp_exp(&exponential_growth(&[(3.0, 2.0)])), - Some(Ordering::Less) - ); - assert_eq!( - exponential_growth(&[(2.0, 3.0)]).cmp_exp(&exponential_growth(&[(3.0, 1.0)])), - None - ); - - assert_eq!( - exponential_growth(&[(0.25, -1.0)]).cmp_exp(&exponential_growth(&[(0.5, -1.0)])), - Some(Ordering::Greater) - ); - assert_eq!( - exponential_growth(&[(0.5, -1.0)]).cmp_exp(&exponential_growth(&[(0.25, -1.0)])), - Some(Ordering::Less) - ); - assert_eq!( - exponential_growth(&[(0.25, -2.0)]).cmp_exp(&exponential_growth(&[(0.5, -1.0)])), - Some(Ordering::Greater) - ); - assert_eq!( - exponential_growth(&[(0.5, -1.0)]).cmp_exp(&exponential_growth(&[(0.25, -2.0)])), - Some(Ordering::Less) - ); - assert_eq!( - exponential_growth(&[(0.25, -1.0)]).cmp_exp(&exponential_growth(&[(0.5, -2.0)])), - None - ); - assert_eq!(two.cmp_exp(&exponential_growth(&[(0.5, -1.0)])), None); - - let natural = VariableGrowth::exponential(ExpBase::Natural, BigRational::one()); - assert_eq!(natural.cmp_exp(&two), Some(Ordering::Greater)); - assert_eq!(two.cmp_exp(&natural), Some(Ordering::Less)); - + let cases = [ + ("2^(2*n)", "2^n", Some(Ordering::Greater)), + ("3^n", "2^n", Some(Ordering::Greater)), + ("3^(2*n)", "2^n", Some(Ordering::Greater)), + ("2^(3*n)", "3^n", Some(Ordering::Greater)), + ("0.25^(-n)", "0.5^(-n)", Some(Ordering::Greater)), + ("0.25^(-2*n)", "0.5^(-n)", Some(Ordering::Greater)), + ("0.25^(-n)", "0.5^(-2*n)", Some(Ordering::Equal)), + ("2^n", "0.5^(-n)", Some(Ordering::Equal)), + ("4^n", "2^(2*n)", Some(Ordering::Equal)), + ("exp(n)", "2^n", Some(Ordering::Greater)), + ("2^(2*n) * 3^n", "2^n * 4^n", Some(Ordering::Greater)), + ("2^n * m", "2^n * n", None), + ("2^(n*m)", "2^(n^2)", None), + ]; + for (left, right, expected) in cases { + let (a, b) = ( + terms_of(&g(left))[0].clone(), + terms_of(&g(right))[0].clone(), + ); + assert_eq!(a.cmp(&b), expected, "{left} vs {right}"); + assert_eq!( + b.cmp(&a), + expected.map(Ordering::reverse), + "{right} vs {left}" + ); + } // Constant subtrees normalize before growth comparison. assert_eq!(g("(1 + 2)^n"), g("3^n")); - - // Two residual products with no factorwise proof remain incomparable. - assert_eq!( - exponential_growth(&[(2.0, 2.0), (3.0, 1.0)]) - .cmp_exp(&exponential_growth(&[(2.0, 1.0), (4.0, 1.0)])), - None - ); } #[test] fn test_exponential_product_canonicalization() { - let combined = exponential_growth(&[(2.0, 1.0)]) - .mul(&exponential_growth(&[(2.0, 2.0)])) - .unwrap(); - assert_eq!(combined, exponential_growth(&[(2.0, 3.0)])); - - let cancelled = exponential_growth(&[(2.0, 1.0)]) - .mul(&exponential_growth(&[(2.0, -1.0)])) - .unwrap(); - assert!(cancelled.is_empty()); + assert_eq!(g("2^n * 2^(2*n)"), g("2^(3*n)")); + assert_eq!(g("2^n * 2^(-n)"), g("1")); + assert_eq!(g("4^n * 2^(-2*n)"), g("1")); } #[test] @@ -333,8 +313,204 @@ fn test_proven_equal_exponential_spelling_is_deterministic() { /// absolute-value idiom. #[test] fn test_growth_widening() { - assert_eq!(g("n - m"), g("n + m")); - assert_eq!(g("sqrt((n - m)^2)"), g("n + m")); + assert_eq!(terms_of(&g("n - m")), terms_of(&g("n + m"))); + assert_eq!(terms_of(&g("sqrt((n - m)^2)")), terms_of(&g("n + m"))); + assert_eq!(g("n - m").relation(), Some(ParameterRelation::UpperBound)); +} + +/// Exactness: a known growth is tight unless a step widened it. +#[test] +fn test_growth_relation() { + for source in [ + "n", + "3", + "n^2 + m", + "n - 1", + "n^2 - n*log(n)", + "2^n * n^3 - n^5", + "sqrt(n*m)", + "log(n)", + "log(n + 1)", + "log(2^n * m)", + "2^(n - log(n))", + "(2*n)^n", + "n^(n*m)", + "factorial(n)", + "factorial(n^2/2)", + ] { + assert_eq!( + g(source).relation(), + Some(ParameterRelation::Exact), + "{source}" + ); + } + for source in [ + "n - m", + "-n + 2", + "n*m - n^2", + "(n + 1)^n", + "2^log(n + m)", + "2^((n + m)^(1/2))", + "log(2^n * log(m))", + "log(1)", + "log(1 + 2^(-n))", + "factorial(n + m)", + "factorial(n * log(n))", + ] { + assert_eq!( + g(source).relation(), + Some(ParameterRelation::UpperBound), + "{source}" + ); + } + assert_eq!(g("2^(2^n)").relation(), None); + // The inner iterated logarithm is widened; its upper bound alone cannot + // bound the magnitude after applying another logarithm. + let nested_log = g("log(log(log(n)))"); + assert_eq!(nested_log.relation(), None); + assert!( + nested_log + .failures() + .expect("unsupported nested logarithm") + .iter() + .any(|reason| matches!(reason, GrowthFailure::UnprovedLogarithmBound(_))) + ); +} + +/// `log` of a decaying term is large and negative; its growth must bound the +/// magnitude, so it can neither be dropped from a sum nor vanish under a power. +#[test] +fn test_growth_log_of_decaying_term_bounds_magnitude() { + assert_eq!(terms_of(&g("log(exp(-n))")), terms_of(&g("n"))); + assert_eq!(terms_of(&g("log(exp(-n))^2")), terms_of(&g("n^2"))); + assert_eq!(terms_of(&g("log(2^n * 3^(-m))")), terms_of(&g("n + m"))); + for source in [ + "n + log(exp(-n)) + 1", + "n^2 + n*log(exp(-n))", + "log(n) + log(exp(-log(n)))", + "n^(1/2) + log(exp(-(n^(1/2))))", + "n + log(exp(-n))^2", + ] { + assert_eq!( + g(source).relation(), + Some(ParameterRelation::UpperBound), + "{source}" + ); + } + // A dominated decaying logarithm still leaves the sum exact. + assert_eq!( + g("n^2 + log(exp(-n))").relation(), + Some(ParameterRelation::Exact) + ); +} + +/// Algebraic cancellation can leave an exponentially small positive argument. +/// An upper bound on that argument is not a magnitude bound for its logarithm. +#[test] +fn test_growth_log_of_widened_argument_is_unknown() { + let argument = "n*(n+1) - n^2 - n + exp(-n)"; + for source in [ + format!("log({argument})"), + format!("log({argument})^2"), // Exactly n^2. + format!("n + log({argument}) + 1"), // Exactly 1, never Theta(n). + "log(n - m)^2".to_string(), + "log(sqrt((n - m)^2))".to_string(), + ] { + let result = g(&source); + assert!( + matches!( + result.failures(), + Some([GrowthFailure::UnprovedLogarithmBound(_)]) + ), + "{source}: {result:?}" + ); + assert_eq!(result.relation(), None, "{source}"); + assert_eq!(result.to_expr(), None, "{source}"); + assert!(!result.dominates(&g("n")), "{source}"); + assert!(!g("n").dominates(&result), "{source}"); + } +} + +/// A two-sided argument bound still permits growing, bounded, and decaying logs. +#[test] +fn test_growth_log_of_exact_argument_remains_supported() { + for (source, bound) in [ + ("log(n + m)", "log(n) + log(m)"), + ("log(n^2 - n)", "log(n)"), + ("log(1 + exp(-n))", "1"), + ("log(exp(-n))^2", "n^2"), + ] { + assert_eq!(terms_of(&g(source)), terms_of(&g(bound)), "{source}"); + } +} + +/// Signed upper bounds cannot be squared or multiplied as magnitude bounds. +#[test] +fn test_growth_rejects_unproved_signed_exponent_operations() { + for source in [ + "exp(log(sqrt(n^2+1)-n)^2)", // exp(asinh(n)^2), increasing without bound. + "exp((sqrt(n^2+n)-n-sqrt(n))^2)", // Exponent is Theta(n), not zero. + "exp((sqrt(n^2+n)-n-sqrt(n))*(sqrt(n^2+4*n)-n-4*sqrt(n)))", + "(sqrt(n^2+1)-n)^(-n)", // exp(n*asinh(n)), not constant. + "(sqrt(n^2+1)-n)^(sqrt(n^2+1)-n-1)", + ] { + let result = g(source); + assert!(result.failures().is_some(), "{source}: {result:?}"); + assert_eq!(result.relation(), None, "{source}"); + assert_eq!(result.to_expr(), None, "{source}"); + assert!(!g("1").dominates(&result), "{source}"); + } +} + +/// An even power is invariant under negation, also past the exact expansion cap. +#[test] +fn test_growth_large_even_exponent_power_keeps_negative_terms() { + let negative = g("exp((-n-1)^18)"); + assert_eq!(negative, g("exp((n+1)^18)")); + assert!(negative.dominates(&g("exp(n^18)"))); + assert!(!g("1").dominates(&negative)); +} + +/// Binary64 rounds the first base to one; it must not reverse or prune this order. +#[test] +fn test_growth_near_one_bases_use_certified_weight_signs() { + let larger = "(9007199254740993/9007199254740992)^(3*n)"; + let smaller = "(9007199254740994/9007199254740992)^n"; + assert!(g(larger).dominates(&g(smaller))); + assert!(!g(smaller).dominates(&g(larger))); + assert_eq!(g(&format!("{larger}+{smaller}")), g(larger)); + + // This rational base grows faster than exp(n/2^54). If the exact-product + // budget is exceeded, neither an incorrect order nor pruning is allowed. + let tiny_exp = g("exp(n/18014398509481984)"); + let rational = g("(9007199254740993/9007199254740992)^n"); + assert!(!tiny_exp.dominates(&rational)); + assert_eq!(terms_of(&add(tiny_exp, rational.clone())).len(), 2); + + assert!(g("exp(n)").dominates(&g("2^n"))); + assert!(g("3^n").dominates(&g("exp(n)"))); + assert!(!g("exp(n)").dominates(&g("3^n"))); +} + +/// Stirling: `n! ≍ n^(n+1/2)·e^(−n)`, strictly between `(n/e)^n` and `n^n`. +#[test] +fn test_growth_factorial_stirling() { + assert_eq!(g("factorial(n)"), g("n^n * sqrt(n) * exp(-n)")); + assert_eq!(g("factorial(2*n)"), g("(2*n)^(2*n) * sqrt(n) * exp(-2*n)")); + for (lower, upper) in [ + ("(n/3)^n", "factorial(n)"), + ("factorial(n)", "n^n"), + ("factorial(n)", "factorial(n + 1)"), + ("2^(n*log(n))", "factorial(2*n)"), + ] { + assert!(g(upper).dominates(&g(lower)), "{upper} ≻ {lower}"); + assert!(!g(lower).dominates(&g(upper)), "{lower} ⊁ {upper}"); + } + // Without an exact single-monomial argument, x! ≤ x^x is kept. + assert_eq!( + terms_of(&g("factorial(n + m)")), + terms_of(&g("(n + m)^(n + m)")) + ); } /// 6. Determinism: the antichain is canonically sorted, so structurally @@ -351,12 +527,12 @@ fn test_growth_determinism() { #[test] fn test_growth_unknown_negative_control() { assert_eq!( - g("2^(n*k)").failures(), - Some([GrowthFailure::NonlinearExponent("k * n".to_string())].as_slice()) + g("2^(2^n)").failures(), + Some([GrowthFailure::UnsupportedExponent("2^n".to_string())].as_slice()) ); assert!(matches!( - g("factorial(n)").failures(), - Some([GrowthFailure::FactorialOfNonconstant(_)]) + g("factorial(2^n)").failures(), + Some([GrowthFailure::UnsupportedFactorial(_)]) )); assert!(matches!( Growth::from_expr(&Expr::factorial(Expr::rational(7, 2))).failures(), @@ -367,22 +543,22 @@ fn test_growth_unknown_negative_control() { Some([GrowthFailure::InvalidConstantDomain { .. }]) )); assert_eq!( - g("factorial(n) + 2^(n*k)").failures(), + g("factorial(2^n) + 2^(2^n)").failures(), Some( [ - GrowthFailure::NonlinearExponent("k * n".to_string()), - GrowthFailure::FactorialOfNonconstant("factorial(n)".to_string()), + GrowthFailure::UnsupportedExponent("2^n".to_string()), + GrowthFailure::UnsupportedFactorial("factorial(2^n)".to_string()), ] .as_slice() ) ); // Absorption through the real `from_expr` add/mul paths. - let factorial_failure = g("factorial(n)"); - assert_eq!(g("factorial(n) + n^2"), factorial_failure); - assert_eq!(g("n^2 + factorial(n)"), factorial_failure); - assert_eq!(g("factorial(n) * n^2"), factorial_failure); - assert_eq!(g("n^2 * factorial(n)"), factorial_failure); + let factorial_failure = g("factorial(2^n)"); + assert_eq!(g("factorial(2^n) + n^2"), factorial_failure); + assert_eq!(g("n^2 + factorial(2^n)"), factorial_failure); + assert_eq!(g("factorial(2^n) * n^2"), factorial_failure); + assert_eq!(g("n^2 * factorial(2^n)"), factorial_failure); // Absorption at the operation level too. let n2 = g("n^2"); @@ -406,18 +582,18 @@ fn test_growth_reports_nested_and_numeric_failures() { let huge_constant = Expr::parse(&format!("1{}", "0".repeat(400))); assert_eq!(Growth::from_expr(&huge_constant), g("1")); - let unsupported = Expr::factorial(Expr::variable("n")); + let factorial = Expr::factorial(Expr::variable("n")); assert!(matches!( - Growth::from_expr(&Expr::exp(unsupported.clone())).failures(), - Some([GrowthFailure::FactorialOfNonconstant(_)]) + Growth::from_expr(&Expr::exp(factorial.clone())).failures(), + Some([GrowthFailure::UnsupportedExponent(_)]) )); assert!(matches!( - Growth::from_expr(&Expr::factorial(unsupported)).failures(), - Some([GrowthFailure::FactorialOfNonconstant(_)]) + Growth::from_expr(&Expr::factorial(factorial)).failures(), + Some([GrowthFailure::UnsupportedFactorial(_)]) )); assert_eq!(Growth::from_expr(&Expr::variable("n")).failures(), None); - assert_eq!(exact_growth(Vec::new()).to_expr(), Some(Expr::integer(1))); + assert_eq!(exact(Vec::new()).to_expr(), Some(Expr::integer(1))); } /// Pure constants, constant factors, and constant division are all O(1) / dropped. @@ -443,10 +619,10 @@ fn test_growth_pow_special_cases() { g("n^(-1)").failures(), Some([GrowthFailure::NegativeExponent(_)]) )); - // Variable base with variable exponent is not representable. + // A variable base must expand to a polynomial in variables and logs. assert!(matches!( - g("n^m").failures(), - Some([GrowthFailure::VariableBaseAndExponent(_)]) + g("(2^n)^m").failures(), + Some([GrowthFailure::UnsupportedVariableBase(_)]) )); } @@ -457,7 +633,7 @@ fn test_growth_to_big_o() { assert_eq!(g("n^2 + n").to_big_o(), "O(n^2)"); assert_eq!(g("2^n").to_big_o(), "O(2^n)"); assert_eq!(g("5").to_big_o(), "O(1)"); - assert_eq!(g("factorial(n)").to_big_o(), "O(?)"); + assert_eq!(g("2^(2^n)").to_big_o(), "O(?)"); // Renders exactly `O()` for bounded classes. let bounded = g("n * m"); assert_eq!( @@ -495,31 +671,83 @@ fn test_growth_exponential_roundtrip_is_exact() { } } -/// `exp(n)` uses base e; unit bases are constant, while decaying directions -/// remain explicit analysis failures rather than silently widening to O(1). +/// `exp(n)` uses base e; unit bases are constant, and decaying factors are +/// represented exactly so they cancel against growing ones. #[test] fn test_growth_exponential_variants() { // exp(n) is represented directly as e^n: exponential, dominates any polynomial. let en = g("exp(n)"); assert!(en.dominates(&g("n^5"))); - assert!(matches!( - g("2^(n - m)").failures(), - Some([GrowthFailure::DecayingExponential { variable, .. }]) if variable == "m" - )); // Unit base is exactly O(1). assert_eq!(g("1^n"), g("7")); - assert!(matches!( - g("0.5^n").failures(), - Some([GrowthFailure::DecayingExponential { - variable, - coefficient, - .. - }]) if variable == "n" && coefficient == "1" - )); + // A decaying factor is bounded by O(1) and does not survive a sum with it. + assert!(g("1").dominates(&g("0.5^n"))); + assert_eq!(g("0.5^n + 1"), g("1")); + assert!(g("2^n").dominates(&g("2^(n - m)"))); // A fractional base with a negative exponent grows and retains that exact // symbolic base instead of being translated through a common logarithm. assert_eq!(g("0.5^(-n)").to_big_o(), "O(0.5^(-1 * n))"); assert!(g("0.5^(-n)").dominates(&g("n^100"))); + // A decaying base cannot take a widened exponent. + assert!(matches!( + g("0.5^sqrt(n + m)").failures(), + Some([GrowthFailure::DecayingWidenedExponent(_)]) + )); +} + +/// Huge constant powers inside an exponent fail fast instead of computing +/// enormous rationals. +#[test] +fn test_growth_huge_exponent_powers_fail_fast() { + let start = std::time::Instant::now(); + for source in [ + "2^((n + m)^3000000000.5)", + "2^((3*n)^2000000000)", + "2^((n + m)^99999999999)", + ] { + assert!( + matches!( + g(source).failures(), + Some([GrowthFailure::RepresentedExponentOutOfRange(_)]) + ), + "{source}: {:?}", + g(source) + ); + } + assert!( + start.elapsed().as_millis() < 100, + "took {:?}", + start.elapsed() + ); + // Unit coefficients need no budget. + assert_eq!(g("2^(n^2000000000)").failures(), None); +} + +/// Nonlinear exponents, variable bases, and factorials are represented. +#[test] +fn test_growth_superexponential_shapes() { + assert_eq!(g("2^(n*m)").to_big_o(), "O(2^(m * n))"); + assert_eq!(g("k^n").to_big_o(), "O(k^n)"); + assert_eq!( + g("num_groups^(2*num_groups)").to_big_o(), + "O(num_groups^(2 * num_groups))" + ); + // (k + 1)^n ≤ (2k)^n keeps the base-2 factor, which is not O(k^n). + assert_eq!(g("(k + 1)^n").to_big_o(), "O(2^n * k^n)"); + assert!(g("(k + 1)^n").dominates(&g("k^n"))); + assert!(!g("k^n").dominates(&g("(k + 1)^n"))); + + assert!(g("2^(n^2)").dominates(&g("3^n"))); + assert!(g("2^sqrt(n)").dominates(&g("n^100"))); + assert!(g("1.0001^n").dominates(&g("2^sqrt(n)"))); + assert!(g("n^n").dominates(&g("2^(n*log(n))"))); + assert!(g("2^(0.773*n)").dominates(&g("1.1996^n"))); + assert!(!g("2^(n*m)").dominates(&g("2^(n^2)"))); + assert!(!g("2^(n^2)").dominates(&g("2^(n*m)"))); + // Exponent polynomials expand exactly, including cancellation. + assert_eq!(g("2^((n + m)^2 - 2*n*m)"), g("2^(n^2 + m^2)")); + // A logarithm of a sum widens to an upper bound with an explicit constant. + assert!(g("2^(n*log(n + 1))").dominates(&g("2^(n*log(n))"))); } /// `log` lowers each level: log of an exponential is linear, log of a @@ -531,7 +759,7 @@ fn test_growth_log_levels() { assert_eq!(g("log(3^n)"), g("n")); assert_eq!(g("log(exp(n))"), g("n")); // log(n) is a single log term. - assert_eq!(g("log(n)"), exact_growth(vec![term(&[], &[], &[("n", 1)])])); + assert_eq!(g("log(n)"), exact(vec![term(&[], &[], &[("n", 1)])])); // log(n*m) ≍ log n + log m (two summands, not a product). assert_eq!(terms_of(&g("log(n*m)")).len(), 2); // log of a constant is O(1). @@ -540,7 +768,7 @@ fn test_growth_log_levels() { // A mixed monomial's log keeps *every* factor class: log(2^n * m) ≍ n + log m. // The exponential factor must not swallow the polynomial one. let mixed = g("log(2^n * m)"); - let expected = exact_growth(vec![ + let expected = exact(vec![ term(&[], &[("n", 1.0)], &[]), term(&[], &[], &[("m", 1)]), ]); @@ -558,7 +786,7 @@ fn test_growth_log_levels() { #[test] fn test_growth_unknown_dominance() { let n2 = g("n^2"); - let unknown = g("factorial(n)"); + let unknown = g("2^(2^n)"); assert!(!unknown.dominates(&n2)); assert!(!n2.dominates(&unknown)); assert!(!unknown.dominates(&unknown)); @@ -576,17 +804,17 @@ fn test_growth_antichain_cap_reports_overflow() { .map(|variable| term(&[], &[(variable, 1.0)], &[])) .collect(); - let at_cap = exact_growth(terms[..ANTICHAIN_CAP].to_vec()); + let at_cap = exact(terms[..ANTICHAIN_CAP].to_vec()); assert_eq!(terms_of(&at_cap).len(), ANTICHAIN_CAP); - let overflow = exact_growth(terms.clone()); + let overflow = exact(terms.clone()); assert!(matches!( overflow.failures(), Some([GrowthFailure::AntichainLimitExceeded { limit, terms }]) if *limit == ANTICHAIN_CAP && *terms == ANTICHAIN_CAP + 1 )); - let reversed = exact_growth( + let reversed = exact( vars.iter() .rev() .map(|variable| term(&[], &[(variable, 1.0)], &[])) @@ -595,21 +823,20 @@ fn test_growth_antichain_cap_reports_overflow() { assert_eq!(overflow, reversed); } -/// Unproved exponential comparisons obey the same explicit resource limit. +/// Incomparable exponential terms obey the same explicit resource limit. #[test] -fn test_growth_rejects_large_unproved_exponential_antichain() { +fn test_growth_rejects_large_exponential_antichain() { let terms = (1..=ANTICHAIN_CAP + 1) .map(|i| { - let mut term = GrowthTerm::one(); - term.insert( - "n".into(), - exponential_growth(&[(2.0, i as f64), (3.0, 1.0 / i as f64)]), - ); - term + term( + &[("n", i as f64), ("m", (ANTICHAIN_CAP + 2 - i) as f64)], + &[], + &[], + ) }) .collect::>(); - let growth = exact_growth(terms.clone()); + let growth = exact(terms); assert!(matches!( growth.failures(), Some([GrowthFailure::AntichainLimitExceeded { limit, terms }]) @@ -743,35 +970,72 @@ fn gen_exponential_base(rng: &mut SplitMix64, bases: &[f64]) -> Expr { expression_from_approximation(bases[rng.below(bases.len() as u64) as usize]) } -fn gen_expr(rng: &mut SplitMix64, depth: u32, exponential_bases: &[f64]) -> Expr { +/// Generator knobs: the exponential bases and the leaf generator. +#[derive(Clone, Copy)] +struct GenConfig { + exponential_bases: &'static [f64], + leaf: fn(&mut SplitMix64) -> Expr, +} + +const STABLE: GenConfig = GenConfig { + exponential_bases: STABLE_EXPONENTIAL_BASES, + leaf: gen_leaf, +}; + +/// Constant leaves can hold off an exact term past the finite observation +/// window (`(192 + sqrt(m))^3` looks like `m^0` until `m ≈ 2^15`), which +/// breaks the numeric lower-bound check but not the upper-bound one. +const STABLE_VARIABLE_LEAVES: GenConfig = GenConfig { + exponential_bases: STABLE_EXPONENTIAL_BASES, + leaf: gen_var, +}; + +const ADVERSARIAL: GenConfig = GenConfig { + exponential_bases: ADVERSARIAL_EXPONENTIAL_BASES, + leaf: gen_leaf, +}; + +fn gen_expr(rng: &mut SplitMix64, depth: u32, config: GenConfig) -> Expr { + let exponential_bases = config.exponential_bases; if depth == 0 { - return gen_leaf(rng); + return (config.leaf)(rng); } match rng.below(100) { - 0..=19 => gen_leaf(rng), - 20..=39 => { - gen_expr(rng, depth - 1, exponential_bases) - + gen_expr(rng, depth - 1, exponential_bases) - } - 40..=54 => { - gen_expr(rng, depth - 1, exponential_bases) - * gen_expr(rng, depth - 1, exponential_bases) - } + 0..=19 => (config.leaf)(rng), + 20..=39 => gen_expr(rng, depth - 1, config) + gen_expr(rng, depth - 1, config), + 40..=54 => gen_expr(rng, depth - 1, config) * gen_expr(rng, depth - 1, config), 55..=69 => Expr::pow( - gen_expr(rng, depth - 1, exponential_bases), + gen_expr(rng, depth - 1, config), Expr::integer(1 + rng.below(3)), ), - 70..=79 => Expr::sqrt(gen_expr(rng, depth - 1, exponential_bases)), - 80..=89 => Expr::log(gen_expr(rng, depth - 1, exponential_bases)), - 90..=96 => Expr::pow( + 70..=79 => Expr::sqrt(gen_expr(rng, depth - 1, config)), + 80..=89 => Expr::log(gen_expr(rng, depth - 1, config)), + 90..=95 => Expr::pow( gen_exponential_base(rng, exponential_bases), gen_linear(rng), ), - 97..=98 => Expr::exp(gen_var(rng)), - // ~1% per node: a nonlinear exponent → Unknown (a minority of trees). + 96 => Expr::exp(gen_var(rng)), + // Signed shapes: a subtracted lower-order term (exact through the + // dominance rule) and the negative logarithm of a decaying term. The + // numeric harness cannot witness an undominated decaying logarithm — + // `exp(-x)` underflows f64 past `x ≈ 745` — so that case is pinned by + // `test_growth_log_of_decaying_term_bounds_magnitude`. + 97 => match rng.below(2) { + 0 => gen_var(rng) * gen_var(rng) - gen_var(rng), + _ => Expr::log(Expr::exp(Expr::integer(0) - gen_linear(rng))), + }, + 98 => match rng.below(3) { + 0 => Expr::pow( + gen_exponential_base(rng, exponential_bases), + gen_nonlinear(rng), + ), + 1 => Expr::pow(gen_var(rng), gen_var(rng)), + _ => Expr::factorial(gen_var(rng)), + }, + // ~1% per node: an iterated exponential → Unknown (a minority of trees). _ => Expr::pow( gen_exponential_base(rng, exponential_bases), - gen_nonlinear(rng), + Expr::pow(Expr::integer(2), gen_var(rng)), ), } } @@ -828,6 +1092,13 @@ struct UbResult { /// `broken_from_expr`. Parameterizing here is what gives the harness teeth: the /// exact same code must accept the sound transfer and reject the broken one. fn run_upper_bound(transfer: fn(&Expr) -> Growth, seed: u64, iters: usize) -> UbResult { + run_bound(transfer, seed, iters, false) +} + +/// With `lower`, check only growths reported [`ParameterRelation::Exact`], and +/// check them in the other direction: `growth ≤ C · e` (non-exact results are +/// counted as `unknown`). +fn run_bound(transfer: fn(&Expr) -> Growth, seed: u64, iters: usize, lower: bool) -> UbResult { // Anchor 2^6; check at 2^8, 2^10, 2^12 — all *larger* than the anchor. let anchor = 64_u64; let large = [256_u64, 1024, 4096]; @@ -837,11 +1108,16 @@ fn run_upper_bound(transfer: fn(&Expr) -> Growth, seed: u64, iters: usize) -> Ub let mut r = UbResult::default(); for _ in 0..iters { - let e = gen_expr(&mut rng, MAX_DEPTH, STABLE_EXPONENTIAL_BASES); + let config = if lower { + STABLE_VARIABLE_LEAVES + } else { + STABLE + }; + let e = gen_expr(&mut rng, MAX_DEPTH, config); let g = transfer(&e); let gexpr = match g.to_expr() { - Some(x) => x, - None => { + Some(x) if !lower || g.relation() == Some(ParameterRelation::Exact) => x, + _ => { r.unknown += 1; continue; } @@ -864,7 +1140,8 @@ fn run_upper_bound(transfer: fn(&Expr) -> Growth, seed: u64, iters: usize) -> Ub r.skipped += 1; continue; } - let c = (ve0 / vg0) * slack; + let (bounded0, bound0) = if lower { (vg0, ve0) } else { (ve0, vg0) }; + let c = (bounded0 / bound0) * slack; let mut conclusive = false; for &s in &large { @@ -881,12 +1158,12 @@ fn run_upper_bound(transfer: fn(&Expr) -> Growth, seed: u64, iters: usize) -> Ub } // Both finite and positive: a real, decidable comparison. conclusive = true; - let bound = c * vg; - if ve > bound { + let (bounded, bound) = if lower { (vg, c * ve) } else { (ve, c * vg) }; + if bounded > bound { r.violations += 1; if r.first_violation.is_none() { r.first_violation = Some(format!( - "e = {e} | g = {gexpr} | s = {s}: eval(e) = {ve} > {c} * {vg} = {bound}" + "e = {e} | g = {gexpr} | s = {s}: {bounded} > {bound} (eval(e) = {ve}, eval(g) = {vg}, C = {c})" )); } } @@ -925,7 +1202,7 @@ fn broken_from_expr(expression: &Expr) -> Growth { None => Growth::from_expr(expression), } } - ExprNode::Log(value) => log_growth(broken_from_expr(value)), + ExprNode::Log(value) => log_growth(broken_from_expr(value), value), _ => Growth::from_expr(expression), } } @@ -971,6 +1248,56 @@ fn test_growth_property_upper_bound_negative_control() { ); } +/// Every growth reported exact is also a lower bound: `growth ≤ C · e`. +#[test] +fn test_growth_property_exact_lower_bound_sound() { + let r = run_bound(Growth::from_expr, MASTER_SEED ^ 0x04, UB_ITERS, true); + assert_eq!( + r.violations, + 0, + "lower-bound violation ({} total); first: {}", + r.violations, + r.first_violation.as_deref().unwrap_or("") + ); + assert!( + r.meaningful >= 3000, + "need >= 3000 meaningful exact checks, got {} (non-exact {}, skipped {})", + r.meaningful, + r.unknown, + r.skipped + ); +} + +/// A deliberately broken transfer function: a top-level `Add` returns the +/// *product* of its summands' growths and claims it exact. This is still an +/// upper bound, but not a lower one, so the lower-bound harness must fail. +fn overgrown_from_expr(expression: &Expr) -> Growth { + match expression.node() { + ExprNode::Add(values) => { + let product = values + .iter() + .map(Growth::from_expr) + .reduce(mul) + .expect("normalized sum has at least two terms"); + match product.failures() { + Some(_) => product, + None => exact(terms_of(&product).to_vec()), + } + } + _ => Growth::from_expr(expression), + } +} + +#[test] +fn test_growth_property_exact_lower_bound_negative_control() { + let r = run_bound(overgrown_from_expr, MASTER_SEED ^ 0x04, UB_ITERS, true); + assert!( + r.violations > 0, + "harness failed to catch the seeded Add bug (meaningful {})", + r.meaningful + ); +} + // --- Contract 2: idempotence --- fn term_approx_eq(x: &GrowthTerm, y: &GrowthTerm) -> bool { @@ -980,7 +1307,7 @@ fn term_approx_eq(x: &GrowthTerm, y: &GrowthTerm) -> bool { fn growth_approx_eq(a: &Growth, b: &Growth) -> bool { match (&a.0, &b.0) { (GrowthState::Unknown(_), GrowthState::Unknown(_)) => true, - (GrowthState::Known(ta), GrowthState::Known(tb)) => { + (GrowthState::Known { terms: ta, .. }, GrowthState::Known { terms: tb, .. }) => { ta.len() == tb.len() && ta.iter().all(|t| tb.iter().any(|u| term_approx_eq(t, u))) && tb.iter().all(|u| ta.iter().any(|t| term_approx_eq(t, u))) @@ -998,7 +1325,7 @@ fn test_growth_property_idempotence() { for _ in 0..UB_ITERS { // Idempotence is purely symbolic, so it can safely exercise bases near // one whose numeric crossover lies far beyond the f64 test window. - let e = gen_expr(&mut rng, MAX_DEPTH, ADVERSARIAL_EXPONENTIAL_BASES); + let e = gen_expr(&mut rng, MAX_DEPTH, ADVERSARIAL); let g = Growth::from_expr(&e); let rendered = match g.to_expr() { Some(x) => x,