peroxide/structure/
polynomial.rs

1#[allow(unused_imports)]
2use crate::structure::matrix::*;
3#[allow(unused_imports)]
4use crate::structure::vector::*;
5use crate::util::useful::*;
6#[cfg(feature = "serde")]
7use serde::{Deserialize, Serialize};
8
9use crate::traits::fp::FPVector;
10use peroxide_num::PowOps;
11use std::cmp::{max, min};
12use std::fmt;
13use std::ops::{Add, Div, Mul, Neg, Sub};
14
15// =============================================================================
16// Polynomial Structure
17// =============================================================================
18/// Polynomial Structure
19#[derive(Debug, Clone, Default)]
20#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
21#[cfg_attr(
22    feature = "rkyv",
23    derive(rkyv::Archive, rkyv::Deserialize, rkyv::Serialize)
24)]
25pub struct Polynomial {
26    pub coef: Vec<f64>,
27}
28
29/// Polynomial Print
30///
31/// # Examples
32/// ```
33/// #[macro_use]
34/// extern crate peroxide;
35/// use peroxide::fuga::*;
36///
37/// fn main() {
38///     let a = poly(c!(1,3,2));
39///     a.print(); //x^2 + 3x + 2
40/// }
41/// ```
42impl fmt::Display for Polynomial {
43    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
44        let mut result = String::new();
45        let l = self.coef.len() - 1;
46
47        if l == 0 {
48            let value = self.coef[0];
49            let temp = choose_shorter_string(format!("{}", value), format!("{:.4}", value));
50            return write!(f, "{}", temp);
51        }
52
53        if l == 1 {
54            let coef_1 = self.coef[0];
55            let coef_0 = self.coef[1];
56
57            if coef_1 == 1. {
58                result.push('x');
59            } else if coef_1 == -1. {
60                result.push_str("-x");
61            } else {
62                let temp = choose_shorter_string(format!("{}x", coef_1), format!("{:.4}x", coef_1));
63                result.push_str(&temp);
64            }
65
66            if coef_0 > 0. {
67                let temp =
68                    choose_shorter_string(format!(" + {}", coef_0), format!(" + {:.4}", coef_0));
69                result.push_str(&temp);
70            } else if coef_0 < 0. {
71                let temp = choose_shorter_string(
72                    format!(" - {}", coef_0.abs()),
73                    format!(" - {:.4}", coef_0.abs()),
74                );
75                result.push_str(&temp);
76            }
77            return write!(f, "{}", result);
78        }
79
80        for i in 0..l + 1 {
81            match i {
82                0 => {
83                    let value = self.coef[i];
84                    if value == 1. {
85                        result.push_str(&format!("x^{}", l));
86                    } else if value == -1. {
87                        result.push_str(&format!("-x^{}", l));
88                    } else {
89                        let temp = choose_shorter_string(
90                            format!("{}x^{}", value, l),
91                            format!("{:.4}x^{}", value, l),
92                        );
93                        result.push_str(&temp);
94                    }
95                }
96                i if i == l => {
97                    let value = self.coef[i];
98                    if value > 0. {
99                        let temp = choose_shorter_string(
100                            format!(" + {}", value),
101                            format!(" + {:.4}", value),
102                        );
103                        result.push_str(&temp);
104                    } else if value < 0. {
105                        let temp = choose_shorter_string(
106                            format!(" - {}", value.abs()),
107                            format!(" - {:.4}", value.abs()),
108                        );
109                        result.push_str(&temp);
110                    }
111                }
112                i if i == l - 1 => {
113                    let value = self.coef[i];
114                    if value == 1. {
115                        result.push_str(" + x");
116                    } else if value > 0. {
117                        let temp = choose_shorter_string(
118                            format!(" + {}x", value),
119                            format!(" + {:.4}x", value),
120                        );
121                        result.push_str(&temp);
122                    } else if value == -1. {
123                        result.push_str(" - x");
124                    } else if value < 0. {
125                        let temp = choose_shorter_string(
126                            format!(" - {}x", value.abs()),
127                            format!(" - {:.4}x", value.abs()),
128                        );
129                        result.push_str(&temp);
130                    }
131                }
132                _ => {
133                    let value = self.coef[i];
134                    if value == 1. {
135                        result.push_str(&format!(" + x^{}", l - i));
136                    } else if value > 0. {
137                        let temp = choose_shorter_string(
138                            format!(" + {}x^{}", value, l - i),
139                            format!(" + {:.4}x^{}", value, l - i),
140                        );
141                        result.push_str(&temp);
142                    } else if value == -1. {
143                        result.push_str(&format!(" - x^{}", l - i));
144                    } else if value < 0. {
145                        let temp = choose_shorter_string(
146                            format!(" - {}x^{}", value.abs(), l - i),
147                            format!(" - {:.4}x^{}", value.abs(), l - i),
148                        );
149                        result.push_str(&temp);
150                    }
151                }
152            }
153        }
154        write!(f, "{}", result)
155    }
156}
157
158impl Polynomial {
159    /// Create Polynomial
160    pub fn new(coef: Vec<f64>) -> Self {
161        Self { coef }
162    }
163
164    /// Evaluate polynomial with value according to Horner's method
165    ///
166    /// # Examples
167    /// ```
168    /// #[macro_use]
169    /// extern crate peroxide;
170    /// use peroxide::fuga::*;
171    ///
172    /// fn main() {
173    ///     let a = poly(c!(1,3,2));
174    ///     assert_eq!(a.eval(1), 6_f64);
175    ///
176    ///     let b = poly(c!(1, 1, -2, -2));
177    ///     let x = 2_f64.sqrt();
178    ///     let horner_evaluation = b.eval(x);
179    ///     let naive_evaluation = x.powf(3.0) + x.powf(2.0) - 2.0*x - 2.0;
180    ///     assert_eq!(horner_evaluation, 0_f64);
181    ///     assert_ne!(naive_evaluation, horner_evaluation);
182    /// }
183    /// ```
184    pub fn eval<T>(&self, x: T) -> f64
185    where
186        T: Into<f64> + Copy,
187    {
188        let x = x.into();
189        let l = self.coef.len() - 1;
190        let mut s = self.coef[0];
191        for i in 0..l {
192            s = self.coef[i + 1] + x * s;
193        }
194        s
195    }
196
197    pub fn eval_vec(&self, v: Vec<f64>) -> Vec<f64> {
198        v.fmap(|t| self.eval(t))
199    }
200
201    /// Linear transformation of a polynomial by a given x according to Horner's method
202    ///
203    /// # Examples
204    /// ```
205    /// #[macro_use]
206    /// extern crate peroxide;
207    /// use peroxide::fuga::*;
208    ///
209    /// fn main() {
210    ///     let a = poly(c!(1,3,2));
211    ///     let translated = a.translate_x(2);
212    ///
213    ///     assert_eq!(translated.eval(3), 6_f64);
214    /// }
215    /// ```
216    pub fn translate_x<X>(&self, x: X) -> Self
217    where
218        X: Into<f64> + Copy,
219    {
220        let d = Self::new(vec![1f64, x.into()]);
221
222        let mut coef = vec![0f64; self.coef.len()];
223
224        let (mut p, ri) = self.horner_division(&d);
225        coef[self.coef.len() - 1] = ri;
226
227        for i in (0..(self.coef.len() - 1)).rev() {
228            if p.coef.len() == 1 {
229                coef[i] = p.coef[0];
230                break;
231            }
232
233            let t = p.horner_division(&d);
234            coef[i] = t.1;
235            p = t.0;
236        }
237
238        Self::new(coef)
239    }
240
241    pub fn horner_division(&self, other: &Self) -> (Self, f64) {
242        assert_eq!(other.coef.len(), 2usize);
243        assert_eq!(other.coef[0], 1.0f64);
244
245        let mut coef = vec![0f64; self.coef.len() - 1];
246        coef[0] = self.coef[0];
247
248        let d = other.coef[1];
249        for i in 1..coef.len() {
250            coef[i] = self.coef[i] - d * coef[i - 1];
251        }
252
253        let remainder = self.coef[self.coef.len() - 1] - d * coef[coef.len() - 1];
254        (Self::new(coef), remainder)
255    }
256}
257
258/// Convenient to declare polynomial
259pub fn poly(coef: Vec<f64>) -> Polynomial {
260    Polynomial::new(coef)
261}
262
263// =============================================================================
264// std::ops for Polynomial
265// =============================================================================
266
267impl Neg for Polynomial {
268    type Output = Self;
269
270    fn neg(self) -> Self::Output {
271        Self::new(
272            self.coef
273                .clone()
274                .into_iter()
275                .map(|x| -x)
276                .collect::<Vec<f64>>(),
277        )
278    }
279}
280
281impl Add<Polynomial> for Polynomial {
282    type Output = Self;
283    fn add(self, other: Self) -> Self {
284        let (l1, l2) = (self.coef.len(), other.coef.len());
285        let l_max = max(l1, l2);
286        let l_min = min(l1, l2);
287        let v_max = choose_longer_vec(&self.coef, &other.coef);
288        let v_min = choose_shorter_vec(&self.coef, &other.coef);
289        let mut coef = vec![0f64; l_max];
290
291        for i in 0..l_max {
292            if i < l_max - l_min {
293                coef[i] = v_max[i];
294            } else {
295                let j = i - (l_max - l_min);
296                coef[i] = v_max[i] + v_min[j];
297            }
298        }
299        Self::new(coef)
300    }
301}
302
303impl<T> Add<T> for Polynomial
304where
305    T: Into<f64> + Copy,
306{
307    type Output = Self;
308    fn add(self, other: T) -> Self {
309        let mut new_coef = self.coef.clone();
310        new_coef[self.coef.len() - 1] += other.into();
311        Self::new(new_coef)
312    }
313}
314
315impl Sub<Polynomial> for Polynomial {
316    type Output = Self;
317    fn sub(self, other: Self) -> Self {
318        self.add(other.neg())
319    }
320}
321
322impl<T> Sub<T> for Polynomial
323where
324    T: Into<f64> + Copy,
325{
326    type Output = Self;
327    fn sub(self, other: T) -> Self {
328        let mut new_coef = self.coef.clone();
329        new_coef[self.coef.len() - 1] -= other.into();
330        Self::new(new_coef)
331    }
332}
333
334impl<T> Mul<T> for Polynomial
335where
336    T: Into<f64> + Copy,
337{
338    type Output = Self;
339    fn mul(self, other: T) -> Self {
340        Self::new(
341            self.coef
342                .into_iter()
343                .map(|x| x * other.into())
344                .collect::<Vec<f64>>(),
345        )
346    }
347}
348
349impl Mul<Polynomial> for Polynomial {
350    type Output = Self;
351    fn mul(self, other: Self) -> Self {
352        let (l1, l2) = (self.coef.len(), other.coef.len());
353        let (n1, n2) = (l1 - 1, l2 - 1);
354        let n = n1 + n2;
355        let mut result = vec![0f64; n + 1];
356
357        for i in 0..l1 {
358            let fixed_val = self.coef[i];
359            let fixed_exp = n1 - i;
360
361            for j in 0..l2 {
362                let target_val = other.coef[j];
363                let target_exp = n2 - j;
364
365                let result_val = fixed_val * target_val;
366                let result_exp = fixed_exp + target_exp;
367
368                result[n - result_exp] += result_val;
369            }
370        }
371
372        Self::new(result)
373    }
374}
375
376impl<T> Div<T> for Polynomial
377where
378    T: Into<f64> + Copy,
379{
380    type Output = Self;
381    fn div(self, other: T) -> Self {
382        let val = other.into();
383        assert_ne!(val, 0f64);
384
385        Self::new(self.coef.fmap(|x| x / val))
386    }
387}
388
389impl Div<Polynomial> for Polynomial {
390    type Output = (Self, Self);
391    fn div(self, other: Self) -> Self::Output {
392        let l1 = self.coef.len();
393        let l2 = other.coef.len();
394        assert!(l1 >= l2);
395
396        let mut temp = self.clone();
397        let mut quot_vec: Vec<f64> = Vec::new();
398        let denom = other.coef[0];
399
400        while temp.coef.len() >= l2 {
401            let l = temp.coef.len();
402            let target = temp.coef[0];
403            let nom = target / denom;
404            quot_vec.push(nom);
405            let mut temp_vec = vec![0f64; l - 1];
406
407            for i in 1..l {
408                if i < l2 {
409                    temp_vec[i - 1] = temp.coef[i] - nom * other.coef[i];
410                } else {
411                    temp_vec[i - 1] = temp.coef[i];
412                }
413            }
414
415            temp = poly(temp_vec);
416        }
417
418        let rem = temp;
419
420        (poly(quot_vec), rem)
421    }
422}
423
424impl Mul<Polynomial> for usize {
425    type Output = Polynomial;
426
427    fn mul(self, rhs: Polynomial) -> Self::Output {
428        rhs.mul(self as f64)
429    }
430}
431
432impl Mul<Polynomial> for i32 {
433    type Output = Polynomial;
434
435    fn mul(self, rhs: Polynomial) -> Self::Output {
436        rhs.mul(self as f64)
437    }
438}
439
440impl Mul<Polynomial> for i64 {
441    type Output = Polynomial;
442
443    fn mul(self, rhs: Polynomial) -> Self::Output {
444        rhs.mul(self as f64)
445    }
446}
447
448impl Mul<Polynomial> for f32 {
449    type Output = Polynomial;
450
451    fn mul(self, rhs: Polynomial) -> Self::Output {
452        rhs.mul(self as f64)
453    }
454}
455
456impl Mul<Polynomial> for f64 {
457    type Output = Polynomial;
458
459    fn mul(self, rhs: Polynomial) -> Self::Output {
460        rhs.mul(self)
461    }
462}
463
464// =============================================================================
465// Extra operations for Polynomial
466// =============================================================================
467impl PowOps for Polynomial {
468    type Float = f64;
469    fn powi(&self, n: i32) -> Self {
470        let mut result = self.clone();
471        for _i in 0..n - 1 {
472            result = result * self.clone();
473        }
474        result
475    }
476
477    fn powf(&self, _f: f64) -> Self {
478        unimplemented!()
479    }
480
481    fn pow(&self, _f: Self) -> Self {
482        unimplemented!()
483    }
484
485    fn sqrt(&self) -> Self {
486        unimplemented!()
487    }
488}
489
490// =============================================================================
491// Calculus for Polynomial
492// =============================================================================
493pub trait Calculus {
494    fn derivative(&self) -> Self;
495    fn integral(&self) -> Self;
496    fn integrate<T: Into<f64> + Copy>(&self, interval: (T, T)) -> f64;
497}
498
499impl Calculus for Polynomial {
500    fn derivative(&self) -> Self {
501        let l = self.coef.len() - 1;
502        let mut result = vec![0f64; l];
503
504        for i in 0..l {
505            result[i] = self.coef[i] * (l - i) as f64;
506        }
507        Self::new(result)
508    }
509
510    fn integral(&self) -> Self {
511        let l = self.coef.len();
512        let mut result = vec![0f64; l + 1];
513
514        for i in 0..l {
515            result[i] = self.coef[i] / (l - i) as f64;
516        }
517        Self::new(result)
518    }
519
520    fn integrate<T: Into<f64> + Copy>(&self, interval: (T, T)) -> f64 {
521        let (a, b) = (interval.0.into(), interval.1.into());
522        let integral = self.integral();
523        integral.eval(b) - integral.eval(a)
524    }
525}
526
527// =============================================================================
528// Useful Polynomial
529// =============================================================================
530/// Lagrange Polynomial
531pub fn lagrange_polynomial(node_x: Vec<f64>, node_y: Vec<f64>) -> Polynomial {
532    assert_eq!(node_x.len(), node_y.len());
533    let l = node_x.len();
534
535    if l <= 1 {
536        panic!("Lagrange Polynomial needs at least 2 nodes");
537    } else if l == 2 {
538        let p0 = poly(vec![1f64, -node_x[0]]);
539        let p1 = poly(vec![1f64, -node_x[1]]);
540
541        let a = node_y[1] / (node_x[1] - node_x[0]);
542        let b = -node_y[0] / (node_x[1] - node_x[0]);
543
544        p0 * a + p1 * b
545    } else if l == 3 {
546        let p0 = poly(vec![1f64, -(node_x[0] + node_x[1]), node_x[0] * node_x[1]]);
547        let p1 = poly(vec![1f64, -(node_x[0] + node_x[2]), node_x[0] * node_x[2]]);
548        let p2 = poly(vec![1f64, -(node_x[1] + node_x[2]), node_x[1] * node_x[2]]);
549
550        let a = node_y[2] / ((node_x[2] - node_x[0]) * (node_x[2] - node_x[1]));
551        let b = node_y[1] / ((node_x[1] - node_x[0]) * (node_x[1] - node_x[2]));
552        let c = node_y[0] / ((node_x[0] - node_x[1]) * (node_x[0] - node_x[2]));
553
554        p0 * a + p1 * b + p2 * c
555    } else {
556        let mut result = Polynomial::new(vec![0f64; l]);
557
558        for i in 0..l {
559            let fixed_val = node_x[i];
560            let prod = node_y[i];
561            let mut id = poly(vec![1f64]);
562
563            for j in 0..l {
564                if j == i {
565                    continue;
566                } else {
567                    let target_val = node_x[j];
568                    let denom = fixed_val - target_val;
569                    id = id * (poly(vec![1f64, -target_val]) / denom);
570                }
571            }
572            result = result + (id * prod);
573        }
574        result
575    }
576}
577
578/// Legendre Polynomial
579///
580/// # Description
581/// : Generate `n`-th order of Legendre polynomial
582pub fn legendre_polynomial(n: usize) -> Polynomial {
583    match n {
584        0 => poly(vec![1f64]),       // 1
585        1 => poly(vec![1f64, 0f64]), // x
586        2 => poly(vec![1.5, 0f64, -0.5]),
587        3 => poly(vec![2.5, 0f64, -1.5, 0f64]),
588        _ => {
589            let k = n - 1;
590            let k_f64 = k as f64;
591            ((2f64 * k_f64 + 1f64) * poly(vec![1f64, 0f64]) * legendre_polynomial(k)
592                - k_f64 * legendre_polynomial(k - 1))
593                / (k_f64 + 1f64)
594        }
595    }
596}
597
598#[derive(Debug, Clone, Copy, PartialEq, Eq)]
599pub enum SpecialKind {
600    First,
601    Second,
602}
603
604/// Chebyshev Polynomial
605pub fn chebyshev_polynomial(n: usize, kind: SpecialKind) -> Polynomial {
606    let mut prev = Polynomial::new(vec![1f64]);
607    let mut curr = match kind {
608        SpecialKind::First => Polynomial::new(vec![1f64, 0f64]),
609        SpecialKind::Second => Polynomial::new(vec![2f64, 0f64]),
610    };
611
612    match n {
613        0 => prev,
614        1 => curr,
615        _ => {
616            for _i in 1..n {
617                std::mem::swap(&mut prev, &mut curr);
618                curr = poly(vec![2f64, 0f64]) * prev.clone() - curr;
619            }
620            curr
621        }
622    }
623}
624
625/// Hermite Polynomial
626///
627/// # Description
628/// Generate `n`-th order Hermite polynomial. The physics convention
629/// H_n(x) = (-1)^n e^{x^2}\frac{d^n}{dx^n}e^{-x^2} is used.
630pub fn hermite_polynomial(n: usize) -> Polynomial {
631    let mut prev = Polynomial::new(vec![1f64]); // 1
632    let mut curr = Polynomial::new(vec![2f64, 0f64]); // 2x
633
634    match n {
635        0 => prev,
636        1 => curr,
637        _ => {
638            for idx in 1..n {
639                let k = idx as f64;
640                std::mem::swap(&mut prev, &mut curr);
641                curr = poly(vec![2f64, 0f64]) * prev.clone() - 2.0 * k * curr;
642            }
643            curr
644        }
645    }
646}
647
648/// Bessel Polynomial
649///
650/// # Description
651/// Generate `n`-th order Bessel polynomial. Definition according to
652/// Krall and Fink (1949).
653pub fn bessel_polynomial(n: usize) -> Polynomial {
654    let mut prev = Polynomial::new(vec![1f64]); // 1
655    let mut curr = Polynomial::new(vec![1f64, 1f64]); // x + 1
656
657    match n {
658        0 => prev,
659        1 => curr,
660        _ => {
661            for idx in 1..n {
662                let k = idx as f64;
663                std::mem::swap(&mut prev, &mut curr);
664                curr = (2.0 * k + 1.0) * poly(vec![1f64, 0f64]) * prev.clone() + curr;
665            }
666            curr
667        }
668    }
669}
peroxide/structure/polynomial.rs

peroxide/structure/
polynomial.rs
