qualia_core_db/solvers/
polynomial.rs

1//! Polynomial & complex algebra — the engine's home for complex arithmetic, real
2//! quadratic solving, and dependency-free polynomial root finding.
3//!
4//! This is *not* linear algebra (it was previously co-located with it in a specialized
5//! lib); it is the computer-algebra primitive that matrix-spectral routines
6//! (`solvers::linear_algebra::spectral`) build on. Allocating where outputs are
7//! inherently dynamic (root vectors); all scratch is local.
8
9use crate::solvers::SolversError;
10
11/// A complex number `re + im·i`. Minimal arithmetic for polynomial root finding.
12#[derive(Debug, Clone, Copy, PartialEq)]
13pub struct Complex {
14    pub re: f64,
15    pub im: f64,
16}
17
18impl Complex {
19    pub const fn new(re: f64, im: f64) -> Self {
20        Self { re, im }
21    }
22    pub const fn real(re: f64) -> Self {
23        Self { re, im: 0.0 }
24    }
25
26    #[inline]
27    pub fn add(self, o: Complex) -> Complex {
28        Complex::new(self.re + o.re, self.im + o.im)
29    }
30    #[inline]
31    pub fn sub(self, o: Complex) -> Complex {
32        Complex::new(self.re - o.re, self.im - o.im)
33    }
34    #[inline]
35    pub fn mul(self, o: Complex) -> Complex {
36        Complex::new(
37            self.re * o.re - self.im * o.im,
38            self.re * o.im + self.im * o.re,
39        )
40    }
41    #[inline]
42    pub fn div(self, o: Complex) -> Complex {
43        let d = o.re * o.re + o.im * o.im;
44        Complex::new(
45            (self.re * o.re + self.im * o.im) / d,
46            (self.im * o.re - self.re * o.im) / d,
47        )
48    }
49    /// Modulus |z|.
50    #[inline]
51    pub fn abs(self) -> f64 {
52        self.re.hypot(self.im)
53    }
54    /// True if within `tol` of the real axis.
55    #[inline]
56    pub fn is_real(self, tol: f64) -> bool {
57        self.im.abs() <= tol
58    }
59}
60
61/// The roots of a real quadratic `a·x² + b·x + c = 0`.
62#[derive(Debug, Clone, Copy, PartialEq)]
63pub enum QuadraticRoots {
64    /// Two distinct real roots, ascending.
65    TwoReal(f64, f64),
66    /// One repeated real root (discriminant ≈ 0).
67    DoubleReal(f64),
68    /// A complex conjugate pair `re ± im·i` (im > 0).
69    ComplexPair { re: f64, im: f64 },
70    /// Degenerate leading coefficient (a ≈ 0): the single linear root of `b·x + c = 0`.
71    Linear(f64),
72}
73
74/// Solve `a·x² + b·x + c = 0` over the reals, numerically stably.
75///
76/// Uses the cancellation-avoiding form `q = -(b + sign(b)·√Δ)/2`, roots `q/a` and `c/q`,
77/// for `Δ > 0`; classifies `Δ ≈ 0` as a double root and `Δ < 0` as a complex pair. Falls
78/// back to the linear root when `a ≈ 0`. Returns [`SolversError::ComputationError`] for a
79/// non-finite or fully-degenerate input.
80pub fn solve_quadratic(a: f64, b: f64, c: f64) -> Result<QuadraticRoots, SolversError> {
81    if !(a.is_finite() && b.is_finite() && c.is_finite()) {
82        return Err(SolversError::ComputationError);
83    }
84    let scale = a.abs().max(b.abs()).max(c.abs()).max(1.0);
85
86    // Degenerate leading coefficient → linear (or no/everywhere solution).
87    if a.abs() <= f64::EPSILON * scale {
88        if b.abs() <= f64::EPSILON * scale {
89            return Err(SolversError::ComputationError);
90        }
91        return Ok(QuadraticRoots::Linear(-c / b));
92    }
93
94    let disc = b * b - 4.0 * a * c;
95    let disc_scale = (b * b).max((4.0 * a * c).abs()).max(1.0);
96    if disc.abs() <= 1e-12 * disc_scale {
97        return Ok(QuadraticRoots::DoubleReal(-b / (2.0 * a)));
98    }
99
100    if disc > 0.0 {
101        let sqrt_d = disc.sqrt();
102        let sign_b = if b >= 0.0 { 1.0 } else { -1.0 };
103        let q = -0.5 * (b + sign_b * sqrt_d);
104        let r1 = q / a;
105        let r2 = c / q;
106        let (lo, hi) = if r1 <= r2 { (r1, r2) } else { (r2, r1) };
107        Ok(QuadraticRoots::TwoReal(lo, hi))
108    } else {
109        let re = -b / (2.0 * a);
110        let im = (-disc).sqrt() / (2.0 * a.abs());
111        Ok(QuadraticRoots::ComplexPair { re, im })
112    }
113}
114
115/// Evaluate a polynomial at a complex point via Horner's method.
116/// `coeffs` are in DESCENDING order: `coeffs[0]·x^n + … + coeffs[n]`.
117fn poly_eval_complex(coeffs: &[f64], x: Complex) -> Complex {
118    let mut acc = Complex::real(0.0);
119    for &c in coeffs {
120        acc = acc.mul(x).add(Complex::real(c));
121    }
122    acc
123}
124
125/// Find all complex roots of a real polynomial (DESCENDING coefficients,
126/// `coeffs[0]·x^n + … + coeffs[n]`) via the Durand–Kerner iteration.
127///
128/// Dependency-free and finds all `n` roots simultaneously; suitable for moderate degree.
129/// Leading/trailing zeros are trimmed. Returns `n` roots (real roots have `im ≈ 0`).
130/// Returns [`SolversError::ComputationError`] for a zero or non-finite polynomial.
131pub fn polynomial_roots(coeffs: &[f64]) -> Result<Vec<Complex>, SolversError> {
132    // Trim leading zeros (they do not change the polynomial's degree meaningfully).
133    let start = coeffs
134        .iter()
135        .position(|c| c.abs() > 0.0)
136        .ok_or(SolversError::ComputationError)?;
137    let coeffs = &coeffs[start..];
138    if coeffs.len() == 1 {
139        return Ok(Vec::new()); // a nonzero constant: no roots
140    }
141    if coeffs.iter().any(|c| !c.is_finite()) {
142        return Err(SolversError::ComputationError);
143    }
144
145    // Normalise to monic.
146    let lead = coeffs[0];
147    let monic: Vec<f64> = coeffs.iter().map(|c| c / lead).collect();
148    let degree = monic.len() - 1;
149
150    // Distinct complex initial guesses on a spiral (the classic 0.4 + 0.9i seed).
151    let seed = Complex::new(0.4, 0.9);
152    let mut roots: Vec<Complex> = (0..degree)
153        .map(|k| {
154            let mut z = Complex::real(1.0);
155            for _ in 0..k {
156                z = z.mul(seed);
157            }
158            z
159        })
160        .collect();
161
162    const MAX_ITERS: usize = 500;
163    const TOL: f64 = 1e-14;
164    for _ in 0..MAX_ITERS {
165        let mut max_delta = 0.0_f64;
166        for i in 0..degree {
167            let zi = roots[i];
168            // denominator = Π_{j≠i} (zi - zj)
169            let mut denom = Complex::real(1.0);
170            for j in 0..degree {
171                if j != i {
172                    denom = denom.mul(zi.sub(roots[j]));
173                }
174            }
175            if denom.abs() == 0.0 {
176                continue; // coincident guesses; perturb on the next sweep
177            }
178            let delta = poly_eval_complex(&monic, zi).div(denom);
179            roots[i] = zi.sub(delta);
180            max_delta = max_delta.max(delta.abs());
181        }
182        if max_delta < TOL {
183            break;
184        }
185    }
186
187    // Snap near-real roots to the real axis for clean output.
188    for r in roots.iter_mut() {
189        if r.im.abs() < 1e-9 * (1.0 + r.re.abs()) {
190            r.im = 0.0;
191        }
192    }
193    Ok(roots)
194}
195
196#[cfg(test)]
197mod tests {
198    use super::*;
199
200    #[test]
201    fn quadratic_two_real() {
202        // x² − 3x + 2 = 0 → 1, 2
203        assert_eq!(
204            solve_quadratic(1.0, -3.0, 2.0).unwrap(),
205            QuadraticRoots::TwoReal(1.0, 2.0)
206        );
207    }
208
209    #[test]
210    fn quadratic_double_and_complex_and_linear() {
211        // x² − 2x + 1 → double 1
212        assert_eq!(
213            solve_quadratic(1.0, -2.0, 1.0).unwrap(),
214            QuadraticRoots::DoubleReal(1.0)
215        );
216        // x² + 1 → ±i
217        match solve_quadratic(1.0, 0.0, 1.0).unwrap() {
218            QuadraticRoots::ComplexPair { re, im } => {
219                assert!(re.abs() < 1e-12 && (im - 1.0).abs() < 1e-12);
220            }
221            other => panic!("expected complex pair, got {other:?}"),
222        }
223        // 0·x² + 2x + 4 → linear root −2
224        assert_eq!(
225            solve_quadratic(0.0, 2.0, 4.0).unwrap(),
226            QuadraticRoots::Linear(-2.0)
227        );
228    }
229
230    #[test]
231    fn quadratic_rejects_degenerate() {
232        assert!(matches!(
233            solve_quadratic(0.0, 0.0, 1.0),
234            Err(SolversError::ComputationError)
235        ));
236        assert!(matches!(
237            solve_quadratic(f64::NAN, 1.0, 1.0),
238            Err(SolversError::ComputationError)
239        ));
240    }
241
242    #[test]
243    fn roots_of_known_polynomial() {
244        // (x−1)(x−2)(x−3) = x³ − 6x² + 11x − 6
245        let roots = polynomial_roots(&[1.0, -6.0, 11.0, -6.0]).unwrap();
246        assert_eq!(roots.len(), 3);
247        let mut reals: Vec<f64> = roots.iter().map(|r| r.re).collect();
248        reals.sort_by(|a, b| a.partial_cmp(b).unwrap());
249        for (got, want) in reals.iter().zip([1.0, 2.0, 3.0]) {
250            assert!((got - want).abs() < 1e-6, "{got} != {want}");
251        }
252        assert!(roots.iter().all(|r| r.is_real(1e-6)));
253    }
254
255    #[test]
256    fn roots_complex_pair() {
257        // x² + 1 → ±i
258        let roots = polynomial_roots(&[1.0, 0.0, 1.0]).unwrap();
259        assert_eq!(roots.len(), 2);
260        assert!(roots.iter().any(|r| (r.im.abs() - 1.0).abs() < 1e-6));
261    }
262
263    #[test]
264    fn roots_reject_zero_polynomial() {
265        assert!(matches!(
266            polynomial_roots(&[0.0, 0.0]),
267            Err(SolversError::ComputationError)
268        ));
269    }
270}
qualia_core_db/solvers/polynomial.rs

qualia_core_db/solvers/
polynomial.rs
