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 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 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        pub fn derivative(&self) -> Self {
105            if self.coefficients.len() < 2 {
106                return Self::zero();
107            }
108
109            let coefficients = self
110                .coefficients
111                .iter()
112                .enumerate()
113                .skip(1)
114                .map(|(index, coefficient)| *coefficient * index as f64)
115                .collect();
116
117            Self::new(coefficients)
118        }
119
120        /// Returns the `n`th derivative of the polynomial.
121        #[must_use]
122        pub fn nth_derivative(&self, n: usize) -> Self {
123            if n == 0 {
124                return self.clone();
125            }
126
127            let mut derivative = self.clone();
128            for _ in 0..n {
129                derivative = derivative.derivative();
130                if derivative.is_zero() {
131                    break;
132                }
133            }
134
135            derivative
136        }
137
138        /// Returns the indefinite integral with the provided constant term.
139        #[must_use]
140        pub fn integral(&self, constant: f64) -> Self {
141            let mut coefficients = Vec::with_capacity(self.coefficients.len() + 1);
142            coefficients.push(constant);
143            coefficients.extend(
144                self.coefficients
145                    .iter()
146                    .enumerate()
147                    .map(|(index, coefficient)| *coefficient / (index + 1) as f64),
148            );
149
150            Self::new(coefficients)
151        }
152
153        /// Returns real roots for degree 0, 1, or 2 polynomials.
154        ///
155        /// Higher-degree polynomials return `None`. Constant and zero
156        /// polynomials return `Some(vec![])`.
157        #[must_use]
158        pub fn real_roots_degree_1_or_2(&self) -> Option<Vec<f64>> {
159            match self.degree() {
160                None | Some(0) => Some(Vec::new()),
161                Some(1) => Some(
162                    linear_root(self.coefficients[1], self.coefficients[0])
163                        .into_iter()
164                        .collect(),
165                ),
166                Some(2) => Some(quadratic_roots(
167                    self.coefficients[2],
168                    self.coefficients[1],
169                    self.coefficients[0],
170                )),
171                Some(_) => None,
172            }
173        }
174    }
175
176    impl Add for Polynomial {
177        type Output = Self;
178
179        fn add(self, rhs: Self) -> Self::Output {
180            let max_len = self.coefficients.len().max(rhs.coefficients.len());
181            let coefficients = (0..max_len)
182                .map(|index| {
183                    self.coefficients.get(index).copied().unwrap_or(0.0)
184                        + rhs.coefficients.get(index).copied().unwrap_or(0.0)
185                })
186                .collect();
187
188            Self::new(coefficients)
189        }
190    }
191
192    impl Sub for Polynomial {
193        type Output = Self;
194
195        fn sub(self, rhs: Self) -> Self::Output {
196            let max_len = self.coefficients.len().max(rhs.coefficients.len());
197            let coefficients = (0..max_len)
198                .map(|index| {
199                    self.coefficients.get(index).copied().unwrap_or(0.0)
200                        - rhs.coefficients.get(index).copied().unwrap_or(0.0)
201                })
202                .collect();
203
204            Self::new(coefficients)
205        }
206    }
207
208    impl Mul for Polynomial {
209        type Output = Self;
210
211        fn mul(self, rhs: Self) -> Self::Output {
212            if self.is_zero() || rhs.is_zero() {
213                return Self::zero();
214            }
215
216            let mut coefficients = vec![0.0; self.coefficients.len() + rhs.coefficients.len() - 1];
217
218            for (left_index, left_coefficient) in self.coefficients.iter().enumerate() {
219                for (right_index, right_coefficient) in rhs.coefficients.iter().enumerate() {
220                    coefficients[left_index + right_index] += left_coefficient * right_coefficient;
221                }
222            }
223
224            Self::new(coefficients)
225        }
226    }
227
228    impl Neg for Polynomial {
229        type Output = Self;
230
231        fn neg(self) -> Self::Output {
232            Self::new(
233                self.coefficients
234                    .into_iter()
235                    .map(|coefficient| -coefficient)
236                    .collect(),
237            )
238        }
239    }
240
241    impl Mul<f64> for Polynomial {
242        type Output = Self;
243
244        fn mul(self, rhs: f64) -> Self::Output {
245            Self::new(
246                self.coefficients
247                    .into_iter()
248                    .map(|coefficient| coefficient * rhs)
249                    .collect(),
250            )
251        }
252    }
253
254    impl Div<f64> for Polynomial {
255        type Output = Self;
256
257        fn div(self, rhs: f64) -> Self::Output {
258            Self::new(
259                self.coefficients
260                    .into_iter()
261                    .map(|coefficient| coefficient / rhs)
262                    .collect(),
263            )
264        }
265    }
266}
267
268/// Small real-root helpers for low-degree polynomials.
269pub mod roots {
270    /// Returns the real root of `ax + b = 0`.
271    ///
272    /// Returns `None` when `a == 0.0`.
273    #[must_use]
274    pub fn linear_root(a: f64, b: f64) -> Option<f64> {
275        (a != 0.0).then_some(-b / a)
276    }
277
278    /// Returns the real roots of `ax^2 + bx + c = 0`.
279    ///
280    /// When `a == 0.0`, this falls back to the linear case. Only real roots are
281    /// returned. Distinct real roots are returned in ascending order.
282    #[must_use]
283    pub fn quadratic_roots(a: f64, b: f64, c: f64) -> Vec<f64> {
284        if a == 0.0 {
285            return linear_root(b, c).into_iter().collect();
286        }
287
288        let discriminant = b * b - 4.0 * a * c;
289
290        if discriminant.is_nan() || discriminant < 0.0 {
291            return Vec::new();
292        }
293
294        if discriminant == 0.0 {
295            return vec![-b / (2.0 * a)];
296        }
297
298        let square_root = discriminant.sqrt();
299        let first = (-b - square_root) / (2.0 * a);
300        let second = (-b + square_root) / (2.0 * a);
301
302        if first <= second {
303            vec![first, second]
304        } else {
305            vec![second, first]
306        }
307    }
308}
309
310#[cfg(test)]
311mod tests {
312    use super::{Polynomial, linear_root, quadratic_roots};
313
314    #[test]
315    fn constructor_trims_trailing_zeros() {
316        let polynomial = Polynomial::new(vec![1.0, 2.0, 0.0, 0.0]);
317
318        assert_eq!(polynomial.coefficients(), &[1.0, 2.0]);
319    }
320
321    #[test]
322    fn zero_polynomial_uses_empty_representation() {
323        let zero = Polynomial::new(vec![0.0, 0.0]);
324
325        assert!(zero.coefficients().is_empty());
326        assert_eq!(zero.degree(), None);
327        assert_eq!(zero.leading_coefficient(), None);
328        assert!(zero.is_zero());
329    }
330
331    #[test]
332    fn constant_constructors_work() {
333        assert_eq!(Polynomial::constant(4.0).coefficients(), &[4.0]);
334        assert!(Polynomial::zero().is_zero());
335        assert_eq!(Polynomial::one().coefficients(), &[1.0]);
336    }
337
338    #[test]
339    fn reports_degree_and_leading_coefficient() {
340        let polynomial = Polynomial::new(vec![3.0, 2.0, 1.0]);
341
342        assert_eq!(polynomial.degree(), Some(2));
343        assert_eq!(polynomial.leading_coefficient(), Some(1.0));
344    }
345
346    #[test]
347    fn evaluates_with_horner_method() {
348        let polynomial = Polynomial::new(vec![3.0, 2.0, 1.0]);
349
350        assert_eq!(polynomial.evaluate(2.0), 11.0);
351        assert_eq!(Polynomial::zero().evaluate(3.0), 0.0);
352    }
353
354    #[test]
355    fn derivative_handles_zero_constant_linear_and_quadratic() {
356        assert_eq!(Polynomial::zero().derivative(), Polynomial::zero());
357        assert_eq!(Polynomial::constant(5.0).derivative(), Polynomial::zero());
358        assert_eq!(
359            Polynomial::new(vec![3.0, 2.0]).derivative(),
360            Polynomial::constant(2.0)
361        );
362        assert_eq!(
363            Polynomial::new(vec![3.0, 2.0, 1.0]).derivative(),
364            Polynomial::new(vec![2.0, 2.0])
365        );
366    }
367
368    #[test]
369    fn nth_derivative_behaves_as_expected() {
370        let polynomial = Polynomial::new(vec![1.0, 2.0, 3.0, 4.0]);
371
372        assert_eq!(polynomial.nth_derivative(0), polynomial.clone());
373        assert_eq!(
374            polynomial.nth_derivative(1),
375            Polynomial::new(vec![2.0, 6.0, 12.0])
376        );
377        assert_eq!(
378            polynomial.nth_derivative(2),
379            Polynomial::new(vec![6.0, 24.0])
380        );
381        assert_eq!(polynomial.nth_derivative(4), Polynomial::zero());
382    }
383
384    #[test]
385    fn computes_integral_with_constant() {
386        let polynomial = Polynomial::new(vec![2.0, 6.0]);
387
388        assert_eq!(
389            polynomial.integral(5.0),
390            Polynomial::new(vec![5.0, 2.0, 3.0])
391        );
392    }
393
394    #[test]
395    fn adds_polynomials() {
396        let left = Polynomial::new(vec![1.0, 2.0, 3.0]);
397        let right = Polynomial::new(vec![4.0, -2.0]);
398
399        assert_eq!(left + right, Polynomial::new(vec![5.0, 0.0, 3.0]));
400    }
401
402    #[test]
403    fn subtracts_polynomials() {
404        let left = Polynomial::new(vec![5.0, 1.0]);
405        let right = Polynomial::new(vec![2.0, 3.0, 4.0]);
406
407        assert_eq!(left - right, Polynomial::new(vec![3.0, -2.0, -4.0]));
408    }
409
410    #[test]
411    fn multiplies_polynomials() {
412        let left = Polynomial::new(vec![1.0, 1.0]);
413        let right = Polynomial::new(vec![1.0, -1.0]);
414
415        assert_eq!(left * right, Polynomial::new(vec![1.0, 0.0, -1.0]));
416    }
417
418    #[test]
419    fn negates_polynomials() {
420        let polynomial = Polynomial::new(vec![1.0, -2.0, 3.0]);
421
422        assert_eq!(-polynomial, Polynomial::new(vec![-1.0, 2.0, -3.0]));
423    }
424
425    #[test]
426    fn multiplies_by_scalar() {
427        let polynomial = Polynomial::new(vec![1.0, -2.0, 3.0]);
428
429        assert_eq!(polynomial * 2.0, Polynomial::new(vec![2.0, -4.0, 6.0]));
430        assert_eq!(Polynomial::new(vec![1.0, 2.0]) * 0.0, Polynomial::zero());
431    }
432
433    #[test]
434    fn divides_by_scalar() {
435        let polynomial = Polynomial::new(vec![2.0, -4.0, 6.0]);
436
437        assert_eq!(polynomial / 2.0, Polynomial::new(vec![1.0, -2.0, 3.0]));
438    }
439
440    #[test]
441    fn solves_linear_roots() {
442        assert_eq!(linear_root(2.0, -4.0), Some(2.0));
443        assert_eq!(linear_root(0.0, 3.0), None);
444    }
445
446    #[test]
447    fn solves_quadratic_roots_with_two_real_roots() {
448        assert_eq!(quadratic_roots(1.0, -3.0, 2.0), vec![1.0, 2.0]);
449    }
450
451    #[test]
452    fn solves_quadratic_roots_with_one_repeated_root() {
453        assert_eq!(quadratic_roots(1.0, -2.0, 1.0), vec![1.0]);
454    }
455
456    #[test]
457    fn solves_quadratic_roots_with_no_real_roots() {
458        assert!(quadratic_roots(1.0, 0.0, 1.0).is_empty());
459    }
460
461    #[test]
462    fn quadratic_roots_fall_back_to_linear_case() {
463        assert_eq!(quadratic_roots(0.0, 2.0, -4.0), vec![2.0]);
464        assert!(quadratic_roots(0.0, 0.0, 1.0).is_empty());
465    }
466
467    #[test]
468    fn low_degree_real_root_helper_handles_supported_cases() {
469        assert_eq!(Polynomial::zero().real_roots_degree_1_or_2(), Some(vec![]));
470        assert_eq!(
471            Polynomial::constant(5.0).real_roots_degree_1_or_2(),
472            Some(vec![])
473        );
474        assert_eq!(
475            Polynomial::new(vec![-4.0, 2.0]).real_roots_degree_1_or_2(),
476            Some(vec![2.0])
477        );
478        assert_eq!(
479            Polynomial::new(vec![2.0, -3.0, 1.0]).real_roots_degree_1_or_2(),
480            Some(vec![1.0, 2.0])
481        );
482        assert_eq!(
483            Polynomial::new(vec![1.0, 0.0, 0.0, 1.0]).real_roots_degree_1_or_2(),
484            None
485        );
486    }
487}
use_polynomial/lib.rs

use_polynomial/
lib.rs
