qualia_core_db/specialized_libs/
symbolic_solve.rs

1//! **Equation solving** (Gap analysis §3.2) — polynomial roots (any degree), real-root
2//! extraction, linear systems, and roots of a CAS polynomial expression.
3//!
4//! Reuses the engine's real root finder ([`crate::solvers::polynomial::polynomial_roots`],
5//! Durand–Kerner) and the polynomial least-squares fit
6//! ([`crate::solvers::interpolation::poly_fit`]) — no re-implementation. (The earlier
7//! sub-agent could not reach these because it was built on the wrong branch; here they
8//! exist.)
9
10use crate::solvers::interpolation::poly_fit;
11use crate::solvers::polynomial::{polynomial_roots, Complex};
12use crate::specialized_libs::symbolic_algebra::Expr;
13
14#[derive(Debug, Clone, PartialEq, Eq)]
15pub enum SolveError {
16    /// Degenerate / non-finite polynomial, or a singular linear system.
17    NoSolution,
18}
19
20/// All complex roots of a real polynomial given in **descending** coefficients
21/// (`coeffs[0]·xⁿ + … + coeffs[n]`).
22pub fn roots(coeffs: &[f64]) -> Result<Vec<Complex>, SolveError> {
23    polynomial_roots(coeffs).map_err(|_| SolveError::NoSolution)
24}
25
26/// The **real** roots (those with `|im| < tol`), ascending.
27pub fn real_roots(coeffs: &[f64], tol: f64) -> Result<Vec<f64>, SolveError> {
28    let mut rs: Vec<f64> = roots(coeffs)?
29        .into_iter()
30        .filter(|z| z.im.abs() < tol)
31        .map(|z| z.re)
32        .collect();
33    rs.sort_by(|a, b| a.partial_cmp(b).unwrap_or(core::cmp::Ordering::Equal));
34    Ok(rs)
35}
36
37/// Solve `A x = b` (row-major `n×n`) by Gaussian elimination with partial pivoting.
38/// `None` if singular or shapes are inconsistent.
39pub fn solve_linear_system(a: &[f64], b: &[f64], n: usize) -> Option<Vec<f64>> {
40    if a.len() != n * n || b.len() != n || n == 0 {
41        return None;
42    }
43    let mut m = a.to_vec();
44    let mut rhs = b.to_vec();
45    for col in 0..n {
46        let mut piv = col;
47        let mut best = m[col * n + col].abs();
48        for r in (col + 1)..n {
49            let v = m[r * n + col].abs();
50            if v > best {
51                best = v;
52                piv = r;
53            }
54        }
55        if best < 1e-14 {
56            return None;
57        }
58        if piv != col {
59            for c in 0..n {
60                m.swap(col * n + c, piv * n + c);
61            }
62            rhs.swap(col, piv);
63        }
64        for r in (col + 1)..n {
65            let f = m[r * n + col] / m[col * n + col];
66            for c in col..n {
67                m[r * n + c] -= f * m[col * n + c];
68            }
69            rhs[r] -= f * rhs[col];
70        }
71    }
72    let mut x = vec![0.0; n];
73    for i in (0..n).rev() {
74        let mut s = rhs[i];
75        for j in (i + 1)..n {
76            s -= m[i * n + j] * x[j];
77        }
78        x[i] = s / m[i * n + i];
79    }
80    Some(x)
81}
82
83/// Real roots of a polynomial **expression** of one variable, by sampling it at
84/// `degree+1` points, recovering its power-basis coefficients via [`poly_fit`], and
85/// root-finding. (Exact for a genuine polynomial of the given degree.)
86pub fn solve_polynomial_expr(
87    expr: &Expr,
88    var: &str,
89    degree: usize,
90    tol: f64,
91) -> Result<Vec<f64>, SolveError> {
92    let xs: Vec<f64> = (0..=degree + 1)
93        .map(|i| i as f64 - (degree as f64) / 2.0)
94        .collect();
95    let ys: Vec<f64> = xs
96        .iter()
97        .map(|&x| {
98            let mut env = std::collections::HashMap::new();
99            env.insert(var.to_string(), x);
100            expr.eval(&env).ok_or(SolveError::NoSolution)
101        })
102        .collect::<Result<_, _>>()?;
103    // poly_fit → ascending coeffs; polynomial_roots wants descending.
104    let mut asc = poly_fit(&xs, &ys, degree).map_err(|_| SolveError::NoSolution)?;
105    asc.reverse();
106    real_roots(&asc, tol)
107}
108
109#[cfg(test)]
110mod tests {
111    use super::*;
112    use crate::specialized_libs::symbolic_algebra::{add, c, mul, pow, sub, var};
113
114    fn close(a: &[f64], b: &[f64]) -> bool {
115        a.len() == b.len() && a.iter().zip(b).all(|(x, y)| (x - y).abs() < 1e-6)
116    }
117
118    #[test]
119    fn quadratic_and_cubic_real_roots() {
120        // x² − 5x + 6 → {2, 3}
121        assert!(close(
122            &real_roots(&[1.0, -5.0, 6.0], 1e-6).unwrap(),
123            &[2.0, 3.0]
124        ));
125        // x³ − 6x² + 11x − 6 → {1, 2, 3}
126        assert!(close(
127            &real_roots(&[1.0, -6.0, 11.0, -6.0], 1e-6).unwrap(),
128            &[1.0, 2.0, 3.0]
129        ));
130    }
131
132    #[test]
133    fn complex_roots_filtered_out() {
134        // x² + 1 → no real roots.
135        assert!(real_roots(&[1.0, 0.0, 1.0], 1e-6).unwrap().is_empty());
136    }
137
138    #[test]
139    fn linear_system_solves() {
140        // [[2,1],[1,3]] x = [3,5] → x = [4/5, 7/5]
141        let x = solve_linear_system(&[2.0, 1.0, 1.0, 3.0], &[3.0, 5.0], 2).unwrap();
142        assert!((x[0] - 0.8).abs() < 1e-9 && (x[1] - 1.4).abs() < 1e-9);
143        // Singular → None.
144        assert!(solve_linear_system(&[1.0, 2.0, 2.0, 4.0], &[1.0, 2.0], 2).is_none());
145    }
146
147    #[test]
148    fn roots_from_a_cas_expression() {
149        // f = x² − 5x + 6 as an Expr → {2, 3}
150        let f = add(sub(pow(var("x"), 2), mul(c(5.0), var("x"))), c(6.0));
151        let r = solve_polynomial_expr(&f, "x", 2, 1e-6).unwrap();
152        assert!(close(&r, &[2.0, 3.0]));
153    }
154}
qualia_core_db/specialized_libs/symbolic_solve.rs

qualia_core_db/specialized_libs/
symbolic_solve.rs
