use_polynomial/
lib.rs

1#![forbid(unsafe_code)]
2#![doc = include_str!("../README.md")]
3
4//! Small polynomial primitives and operations for `RustUse`.
5
6#[doc(inline)]
7pub use polynomial::Polynomial;
8#[doc(inline)]
9pub use roots::{linear_root, quadratic_roots};
10
11/// Polynomial primitives and direct operations.
12pub mod polynomial {
13    use core::ops::{Add, Div, Mul, Neg, Sub};
14
15    use crate::roots::{linear_root, quadratic_roots};
16
17    /// A polynomial whose coefficients are stored in ascending power order.
18    ///
19    /// `coefficients[i]` is the coefficient for `x^i`. For example,
20    /// `3 + 2x + x^2` is represented as `vec![3.0, 2.0, 1.0]`.
21    ///
22    /// The zero polynomial is represented as an empty coefficient vector.
23    #[derive(Debug, Clone, PartialEq)]
24    pub struct Polynomial {
25        coefficients: Vec<f64>,
26    }
27
28    impl Polynomial {
29        /// Creates a polynomial from coefficients in ascending power order.
30        ///
31        /// Trailing zero coefficients are trimmed. The zero polynomial is
32        /// stored as an empty coefficient vector.
33        #[must_use]
34        pub fn new(mut coefficients: Vec<f64>) -> Self {
35            while coefficients
36                .last()
37                .is_some_and(|coefficient| *coefficient == 0.0)
38            {
39                coefficients.pop();
40            }
41
42            Self { coefficients }
43        }
44
45        /// Creates a constant polynomial.
46        #[must_use]
47        pub fn constant(value: f64) -> Self {
48            Self::new(vec![value])
49        }
50
51        /// Creates the zero polynomial.
52        #[must_use]
53        pub fn zero() -> Self {
54            Self::new(Vec::new())
55        }
56
57        /// Creates the constant polynomial `1`.
58        #[must_use]
59        pub fn one() -> Self {
60            Self::constant(1.0)
61        }
62
63        /// Returns the coefficients in ascending power order.
64        #[must_use]
65        pub fn coefficients(&self) -> &[f64] {
66            &self.coefficients
67        }
68
69        /// Returns the polynomial degree, or `None` for the zero polynomial.
70        #[must_use]
71        pub const fn degree(&self) -> Option<usize> {
72            if self.coefficients.is_empty() {
73                None
74            } else {
75                Some(self.coefficients.len() - 1)
76            }
77        }
78
79        /// Returns the leading coefficient, or `None` for the zero polynomial.
80        #[must_use]
81        pub fn leading_coefficient(&self) -> Option<f64> {
82            self.coefficients.last().copied()
83        }
84
85        /// Returns `true` when the polynomial is the zero polynomial.
86        #[must_use]
87        pub const fn is_zero(&self) -> bool {
88            self.coefficients.is_empty()
89        }
90
91        /// Evaluates the polynomial at `x` using Horner's method.
92        #[must_use]
93        pub fn evaluate(&self, x: f64) -> f64 {
94            self.coefficients
95                .iter()
96                .rev()
97                .fold(0.0, |accumulator, coefficient| {
98                    accumulator * x + coefficient
99                })
100        }
101
102        /// Returns the derivative of the polynomial.
103        #[must_use]
104        #[allow(clippy::cast_precision_loss)]
105        pub fn derivative(&self) -> Self {
106            if self.coefficients.len() < 2 {
107                return Self::zero();
108            }
109
110            let coefficients = self
111                .coefficients
112                .iter()
113                .enumerate()
114                .skip(1)
115                .map(|(index, coefficient)| *coefficient * index as f64)
116                .collect();
117
118            Self::new(coefficients)
119        }
120
121        /// Returns the `n`th derivative of the polynomial.
122        #[must_use]
123        pub fn nth_derivative(&self, n: usize) -> Self {
124            if n == 0 {
125                return self.clone();
126            }
127
128            let mut derivative = self.clone();
129            for _ in 0..n {
130                derivative = derivative.derivative();
131                if derivative.is_zero() {
132                    break;
133                }
134            }
135
136            derivative
137        }
138
139        /// Returns the indefinite integral with the provided constant term.
140        #[must_use]
141        #[allow(clippy::cast_precision_loss)]
142        pub fn integral(&self, constant: f64) -> Self {
143            let mut coefficients = Vec::with_capacity(self.coefficients.len() + 1);
144            coefficients.push(constant);
145            coefficients.extend(
146                self.coefficients
147                    .iter()
148                    .enumerate()
149                    .map(|(index, coefficient)| *coefficient / (index + 1) as f64),
150            );
151
152            Self::new(coefficients)
153        }
154
155        /// Returns real roots for degree 0, 1, or 2 polynomials.
156        ///
157        /// Higher-degree polynomials return `None`. Constant and zero
158        /// polynomials return `Some(vec![])`.
159        #[must_use]
160        pub fn real_roots_degree_1_or_2(&self) -> Option<Vec<f64>> {
161            match self.degree() {
162                None | Some(0) => Some(Vec::new()),
163                Some(1) => Some(
164                    linear_root(self.coefficients[1], self.coefficients[0])
165                        .into_iter()
166                        .collect(),
167                ),
168                Some(2) => Some(quadratic_roots(
169                    self.coefficients[2],
170                    self.coefficients[1],
171                    self.coefficients[0],
172                )),
173                Some(_) => None,
174            }
175        }
176    }
177
178    impl Add for Polynomial {
179        type Output = Self;
180
181        fn add(self, rhs: Self) -> Self::Output {
182            let max_len = self.coefficients.len().max(rhs.coefficients.len());
183            let coefficients = (0..max_len)
184                .map(|index| {
185                    self.coefficients.get(index).copied().unwrap_or(0.0)
186                        + rhs.coefficients.get(index).copied().unwrap_or(0.0)
187                })
188                .collect();
189
190            Self::new(coefficients)
191        }
192    }
193
194    impl Sub for Polynomial {
195        type Output = Self;
196
197        fn sub(self, rhs: Self) -> Self::Output {
198            let max_len = self.coefficients.len().max(rhs.coefficients.len());
199            let coefficients = (0..max_len)
200                .map(|index| {
201                    self.coefficients.get(index).copied().unwrap_or(0.0)
202                        - rhs.coefficients.get(index).copied().unwrap_or(0.0)
203                })
204                .collect();
205
206            Self::new(coefficients)
207        }
208    }
209
210    impl Mul for Polynomial {
211        type Output = Self;
212
213        fn mul(self, rhs: Self) -> Self::Output {
214            if self.is_zero() || rhs.is_zero() {
215                return Self::zero();
216            }
217
218            let mut coefficients = vec![0.0; self.coefficients.len() + rhs.coefficients.len() - 1];
219
220            for (left_index, left_coefficient) in self.coefficients.iter().enumerate() {
221                for (right_index, right_coefficient) in rhs.coefficients.iter().enumerate() {
222                    coefficients[left_index + right_index] += left_coefficient * right_coefficient;
223                }
224            }
225
226            Self::new(coefficients)
227        }
228    }
229
230    impl Neg for Polynomial {
231        type Output = Self;
232
233        fn neg(self) -> Self::Output {
234            Self::new(
235                self.coefficients
236                    .into_iter()
237                    .map(|coefficient| -coefficient)
238                    .collect(),
239            )
240        }
241    }
242
243    impl Mul<f64> for Polynomial {
244        type Output = Self;
245
246        fn mul(self, rhs: f64) -> Self::Output {
247            Self::new(
248                self.coefficients
249                    .into_iter()
250                    .map(|coefficient| coefficient * rhs)
251                    .collect(),
252            )
253        }
254    }
255
256    impl Div<f64> for Polynomial {
257        type Output = Self;
258
259        fn div(self, rhs: f64) -> Self::Output {
260            Self::new(
261                self.coefficients
262                    .into_iter()
263                    .map(|coefficient| coefficient / rhs)
264                    .collect(),
265            )
266        }
267    }
268}
269
270/// Small real-root helpers for low-degree polynomials.
271pub mod roots {
272    /// Returns the real root of `ax + b = 0`.
273    ///
274    /// Returns `None` when `a == 0.0`.
275    #[must_use]
276    pub fn linear_root(a: f64, b: f64) -> Option<f64> {
277        (a != 0.0).then_some(-b / a)
278    }
279
280    /// Returns the real roots of `ax^2 + bx + c = 0`.
281    ///
282    /// When `a == 0.0`, this falls back to the linear case. Only real roots are
283    /// returned. Distinct real roots are returned in ascending order.
284    #[must_use]
285    pub fn quadratic_roots(a: f64, b: f64, c: f64) -> Vec<f64> {
286        if a == 0.0 {
287            return linear_root(b, c).into_iter().collect();
288        }
289
290        let discriminant = b.mul_add(b, -(4.0 * a * c));
291
292        if discriminant.is_nan() || discriminant < 0.0 {
293            return Vec::new();
294        }
295
296        if discriminant == 0.0 {
297            return vec![-b / (2.0 * a)];
298        }
299
300        let square_root = discriminant.sqrt();
301        let first = (-b - square_root) / (2.0 * a);
302        let second = (-b + square_root) / (2.0 * a);
303
304        if first <= second {
305            vec![first, second]
306        } else {
307            vec![second, first]
308        }
309    }
310}
311
312#[cfg(test)]
313mod tests {
314    use super::{Polynomial, linear_root, quadratic_roots};
315
316    #[test]
317    fn constructor_trims_trailing_zeros() {
318        let polynomial = Polynomial::new(vec![1.0, 2.0, 0.0, 0.0]);
319
320        assert_eq!(polynomial.coefficients(), &[1.0, 2.0]);
321    }
322
323    #[test]
324    fn zero_polynomial_uses_empty_representation() {
325        let zero = Polynomial::new(vec![0.0, 0.0]);
326
327        assert!(zero.coefficients().is_empty());
328        assert_eq!(zero.degree(), None);
329        assert_eq!(zero.leading_coefficient(), None);
330        assert!(zero.is_zero());
331    }
332
333    #[test]
334    fn constant_constructors_work() {
335        assert_eq!(Polynomial::constant(4.0).coefficients(), &[4.0]);
336        assert!(Polynomial::zero().is_zero());
337        assert_eq!(Polynomial::one().coefficients(), &[1.0]);
338    }
339
340    #[test]
341    fn reports_degree_and_leading_coefficient() {
342        let polynomial = Polynomial::new(vec![3.0, 2.0, 1.0]);
343
344        assert_eq!(polynomial.degree(), Some(2));
345        assert_eq!(polynomial.leading_coefficient(), Some(1.0));
346    }
347
348    #[test]
349    fn evaluates_with_horner_method() {
350        let polynomial = Polynomial::new(vec![3.0, 2.0, 1.0]);
351
352        assert!((polynomial.evaluate(2.0) - 11.0).abs() <= f64::EPSILON);
353        assert!(Polynomial::zero().evaluate(3.0).abs() <= f64::EPSILON);
354    }
355
356    #[test]
357    fn derivative_handles_zero_constant_linear_and_quadratic() {
358        assert_eq!(Polynomial::zero().derivative(), Polynomial::zero());
359        assert_eq!(Polynomial::constant(5.0).derivative(), Polynomial::zero());
360        assert_eq!(
361            Polynomial::new(vec![3.0, 2.0]).derivative(),
362            Polynomial::constant(2.0)
363        );
364        assert_eq!(
365            Polynomial::new(vec![3.0, 2.0, 1.0]).derivative(),
366            Polynomial::new(vec![2.0, 2.0])
367        );
368    }
369
370    #[test]
371    fn nth_derivative_behaves_as_expected() {
372        let polynomial = Polynomial::new(vec![1.0, 2.0, 3.0, 4.0]);
373
374        assert_eq!(
375            polynomial.nth_derivative(0),
376            Polynomial::new(vec![1.0, 2.0, 3.0, 4.0])
377        );
378        assert_eq!(
379            polynomial.nth_derivative(1),
380            Polynomial::new(vec![2.0, 6.0, 12.0])
381        );
382        assert_eq!(
383            polynomial.nth_derivative(2),
384            Polynomial::new(vec![6.0, 24.0])
385        );
386        assert_eq!(polynomial.nth_derivative(4), Polynomial::zero());
387    }
388
389    #[test]
390    fn computes_integral_with_constant() {
391        let polynomial = Polynomial::new(vec![2.0, 6.0]);
392
393        assert_eq!(
394            polynomial.integral(5.0),
395            Polynomial::new(vec![5.0, 2.0, 3.0])
396        );
397    }
398
399    #[test]
400    fn adds_polynomials() {
401        let left = Polynomial::new(vec![1.0, 2.0, 3.0]);
402        let right = Polynomial::new(vec![4.0, -2.0]);
403
404        assert_eq!(left + right, Polynomial::new(vec![5.0, 0.0, 3.0]));
405    }
406
407    #[test]
408    fn subtracts_polynomials() {
409        let left = Polynomial::new(vec![5.0, 1.0]);
410        let right = Polynomial::new(vec![2.0, 3.0, 4.0]);
411
412        assert_eq!(left - right, Polynomial::new(vec![3.0, -2.0, -4.0]));
413    }
414
415    #[test]
416    fn multiplies_polynomials() {
417        let left = Polynomial::new(vec![1.0, 1.0]);
418        let right = Polynomial::new(vec![1.0, -1.0]);
419
420        assert_eq!(left * right, Polynomial::new(vec![1.0, 0.0, -1.0]));
421    }
422
423    #[test]
424    fn negates_polynomials() {
425        let polynomial = Polynomial::new(vec![1.0, -2.0, 3.0]);
426
427        assert_eq!(-polynomial, Polynomial::new(vec![-1.0, 2.0, -3.0]));
428    }
429
430    #[test]
431    fn multiplies_by_scalar() {
432        let polynomial = Polynomial::new(vec![1.0, -2.0, 3.0]);
433
434        assert_eq!(polynomial * 2.0, Polynomial::new(vec![2.0, -4.0, 6.0]));
435        assert_eq!(Polynomial::new(vec![1.0, 2.0]) * 0.0, Polynomial::zero());
436    }
437
438    #[test]
439    fn divides_by_scalar() {
440        let polynomial = Polynomial::new(vec![2.0, -4.0, 6.0]);
441
442        assert_eq!(polynomial / 2.0, Polynomial::new(vec![1.0, -2.0, 3.0]));
443    }
444
445    #[test]
446    fn solves_linear_roots() {
447        assert_eq!(linear_root(2.0, -4.0), Some(2.0));
448        assert_eq!(linear_root(0.0, 3.0), None);
449    }
450
451    #[test]
452    fn solves_quadratic_roots_with_two_real_roots() {
453        assert_eq!(quadratic_roots(1.0, -3.0, 2.0), vec![1.0, 2.0]);
454    }
455
456    #[test]
457    fn solves_quadratic_roots_with_one_repeated_root() {
458        assert_eq!(quadratic_roots(1.0, -2.0, 1.0), vec![1.0]);
459    }
460
461    #[test]
462    fn solves_quadratic_roots_with_no_real_roots() {
463        assert!(quadratic_roots(1.0, 0.0, 1.0).is_empty());
464    }
465
466    #[test]
467    fn quadratic_roots_fall_back_to_linear_case() {
468        assert_eq!(quadratic_roots(0.0, 2.0, -4.0), vec![2.0]);
469        assert!(quadratic_roots(0.0, 0.0, 1.0).is_empty());
470    }
471
472    #[test]
473    fn low_degree_real_root_helper_handles_supported_cases() {
474        assert_eq!(Polynomial::zero().real_roots_degree_1_or_2(), Some(vec![]));
475        assert_eq!(
476            Polynomial::constant(5.0).real_roots_degree_1_or_2(),
477            Some(vec![])
478        );
479        assert_eq!(
480            Polynomial::new(vec![-4.0, 2.0]).real_roots_degree_1_or_2(),
481            Some(vec![2.0])
482        );
483        assert_eq!(
484            Polynomial::new(vec![2.0, -3.0, 1.0]).real_roots_degree_1_or_2(),
485            Some(vec![1.0, 2.0])
486        );
487        assert_eq!(
488            Polynomial::new(vec![1.0, 0.0, 0.0, 1.0]).real_roots_degree_1_or_2(),
489            None
490        );
491    }
492}
use_polynomial/lib.rs

use_polynomial/
lib.rs
