hmath/
poly.rs

1use crate::Ratio;
2
3mod from_points;
4
5pub use from_points::{from_points, from_points_generic, cubic_2_points, quadratic_3_points, linear_2_points};
6
7/// [3, 4, 5] -> 3x^2 + 4x + 5
8#[derive(Clone, Debug, PartialEq)]
9pub struct Polynomial {
10    coeffs: Vec<Ratio>,
11}
12
13impl Polynomial {
14    pub fn from_vec(mut coeffs: Vec<Ratio>) -> Self {
15        let mut i = 0;
16
17        while let Some(n) = coeffs.get(i) {
18            if n.is_zero() {
19                i += 1;
20            } else {
21                break;
22            }
23        }
24
25        // [0, 0, 3, 4] = 0x^3 + 0x^2 + 3x + 4 = 3x + 4 = [3, 4]
26        if i != 0 {
27            coeffs = coeffs[i..].to_vec();
28        }
29
30        if coeffs.is_empty() { coeffs = vec![Ratio::zero()]; }
31
32        Polynomial { coeffs }
33    }
34
35    pub fn from_vec_generic<T: Into<Ratio>>(coeffs: Vec<T>) -> Self {
36        if coeffs.is_empty() { return Polynomial::from_vec(vec![]); }
37
38        Polynomial { coeffs: coeffs.into_iter().map(|n| n.into()).collect() }
39    }
40
41    pub fn to_vec(&self) -> &Vec<Ratio> {
42        &self.coeffs
43    }
44
45    /// It panics if a coefficient cannot be converted to an f32. (Eg. not in range)
46    pub fn to_vec_f32(&self) -> Vec<f32> {
47        self.coeffs.iter().map(|n| n.to_ieee754_f32().unwrap()).collect()
48    }
49
50    /// It panics if a coefficient cannot be converted to an f64. (Eg. not in range)
51    pub fn to_vec_f64(&self) -> Vec<f64> {
52        self.coeffs.iter().map(|n| n.to_ieee754_f64().unwrap()).collect()
53    }
54
55    #[must_use = "method returns a new number and does not mutate the original value"]
56    pub fn add(&self, other: &Polynomial) -> Self {
57        todo!()
58    }
59
60    pub fn add_mut(&mut self, other: &Polynomial) {
61        todo!()
62    }
63
64    #[must_use = "method returns a new number and does not mutate the original value"]
65    pub fn mul_k<T: Into<Ratio>>(&self, k: T) -> Self {
66        let k = k.into();
67
68        let result = Polynomial::from_vec(self.coeffs.iter().map(
69            |n| n.mul_rat(&k)
70        ).collect());
71
72        #[cfg(test)] {
73            let mut s = self.clone();
74            s.mul_k_mut(k);
75
76            assert_eq!(s, result);
77        }
78
79        result
80    }
81
82    pub fn mul_k_mut<T: Into<Ratio>>(&mut self, k: T) {
83        let k = k.into();
84
85        if k.is_zero() {
86            *self = Polynomial::from_vec(vec![]);
87            return;
88        }
89
90        for n in self.coeffs.iter_mut() {
91            n.mul_rat_mut(&k);
92        }
93    }
94
95    #[must_use = "method returns a new number and does not mutate the original value"]
96    pub fn mul(&self, other: &Polynomial) -> Self {
97        todo!()
98    }
99
100    pub fn mul_mut(&mut self, other: &Polynomial) {
101        todo!()
102    }
103
104    #[must_use = "method returns a new number and does not mutate the original value"]
105    pub fn differentiate(&self) -> Self {
106        let mut result = Vec::with_capacity(self.coeffs.len());
107
108        for (ind, value) in self.coeffs.iter().rev().enumerate().rev() {
109            result.push(value.mul_i32(ind as i32));
110        }
111
112        result.pop().unwrap();
113
114        if result.is_empty() {
115            result.push(Ratio::zero());
116        }
117
118        let result = Polynomial::from_vec(result);
119
120        #[cfg(test)] {
121            let mut p = self.clone();
122            p.differentiate_mut();
123
124            assert_eq!(result, p);
125        }
126
127        result
128    }
129
130    pub fn differentiate_mut(&mut self) {
131        for (ind, value) in self.coeffs.iter_mut().rev().enumerate() {
132            value.mul_i32_mut(ind as i32);
133        }
134
135        self.coeffs.pop().unwrap();
136
137        if self.coeffs.is_empty() {
138            self.coeffs.push(Ratio::zero());
139        }
140    }
141
142    /// f(x)
143    pub fn calc(&self, x: &Ratio) -> Ratio {
144        let mut result = Ratio::zero();
145
146        for coeff in self.coeffs.iter() {
147            result.mul_rat_mut(x);
148            result.add_rat_mut(coeff);
149        }
150
151        result
152    }
153
154    /// It returns `x - (f(x)/f'(x))`. If you pre-calculated `f'`, pass it to `prime`. Set `prime` to `None` otherwise.
155    pub fn newton_method(&self, x: &Ratio, prime: &Option<Polynomial>) -> Ratio {
156        let fpx = match prime {
157            Some(p) => p.calc(x),
158            None => {
159                let fp = self.differentiate();
160                fp.calc(x)
161            }
162        };
163
164        x.sub_rat(&self.calc(x).div_rat(&fpx))
165    }
166
167    pub fn to_approx_string(&self, max_len: usize) -> String {
168        self.to_string(Some(max_len))
169    }
170
171    pub fn to_ratio_string(&self) -> String {
172        self.to_string(None)
173    }
174
175    fn to_string(&self, approx: Option<usize>) -> String {
176        let mut result = Vec::with_capacity(self.coeffs.len());
177
178        for (ind, value) in self.coeffs.iter().rev().enumerate().rev() {
179            if value.is_zero() {
180                continue;
181            }
182
183            let abs_val = value.abs();
184
185            if !result.is_empty() {
186                if value.is_neg() {
187                    result.push(" - ".to_string());
188                } else {
189                    result.push(" + ".to_string());
190                }
191            } else if value.is_neg() {
192                result.push("-".to_string());
193            }
194
195            if !abs_val.is_one() {
196                if let Some(n) = approx {
197                    result.push(abs_val.to_approx_string(n));
198                } else {
199                    result.push(abs_val.to_ratio_string());
200                }
201
202                if ind > 0 {
203                    result.push(" * ".to_string());
204                }
205            } else if ind == 0 {
206                if let Some(n) = approx {
207                    result.push(abs_val.to_approx_string(n));
208                } else {
209                    result.push(abs_val.to_ratio_string());
210                }
211            }
212
213            if ind != 0 {
214                result.push("x".to_string());
215            }
216
217            if ind > 1 {
218                result.push(format!("^{ind}"));
219            }
220        }
221
222        if result.is_empty() {
223            result.push("0".to_string());
224        }
225
226        result.concat()
227    }
228}
229
230#[cfg(test)]
231mod tests {
232    use crate::{Polynomial, Ratio};
233
234    #[test]
235    fn newtons_method() {
236        let f = Polynomial::from_vec_generic(vec![1, 0, -10]);
237        let fp = Some(f.differentiate());
238        let mut n = Ratio::from_i32(3);
239
240        for _ in 0..6 {
241            n = f.newton_method(&n, &fp);
242        }
243
244        assert_eq!("3.162277660168379331998893544432", n.to_approx_string(32));
245    }
246
247    #[test]
248    fn diff_test() {
249        // 3x^3 + 4x^2 + 5x + 6
250        let p1 = Polynomial::from_vec_generic(vec![3, 4, 5, 6]);
251
252        // 9x^2 + 8x + 5
253        let p2 = Polynomial::from_vec_generic(vec![9, 8, 5]);
254
255        // 18x + 8
256        let p3 = Polynomial::from_vec_generic(vec![18, 8]);
257
258        // 18
259        let p4 = Polynomial::from_vec_generic(vec![18]);
260
261        // 0
262        let p5 = Polynomial::from_vec_generic(vec![0]);
263
264        assert_eq!(p1.differentiate(), p2);
265        assert_eq!(p2.differentiate(), p3);
266        assert_eq!(p3.differentiate(), p4);
267        assert_eq!(p4.differentiate(), p5);
268        assert_eq!(p5.differentiate(), p5);
269    }
270
271    #[test]
272    fn to_string_test() {
273        let v = vec![3.5, 4.25, 5.0, 6.5, 0.0, 1.0, 2.0, -3.0, -4.0];
274        let v = v.into_iter().map(|n| Ratio::from_ieee754_f32(n).unwrap()).collect();
275    
276        assert_eq!("7/2 * x^8 + 17/4 * x^7 + 5 * x^6 + 13/2 * x^5 + x^3 + 2 * x^2 - 3 * x - 4", Polynomial::from_vec(v).to_ratio_string());
277        assert_eq!("0", Polynomial::from_vec_generic(vec![0]).to_ratio_string());
278        assert_eq!("1", Polynomial::from_vec_generic(vec![1]).to_ratio_string());
279        assert_eq!("-1", Polynomial::from_vec_generic(vec![-1]).to_ratio_string());
280        assert_eq!("0", Polynomial::from_vec_generic(vec![0; 3]).to_ratio_string());
281        assert_eq!("x^2 + x + 1", Polynomial::from_vec_generic(vec![1; 3]).to_ratio_string());
282        assert_eq!("-x^2 - x - 1", Polynomial::from_vec_generic(vec![-1; 3]).to_ratio_string());
283        assert_eq!("x^3 + x", Polynomial::from_vec_generic(vec![1, 0, 1, 0]).to_ratio_string());
284        assert_eq!("x^2 + 1", Polynomial::from_vec_generic(vec![0, 1, 0, 1]).to_ratio_string());
285        assert_eq!("x^3 + x", Polynomial::from_vec_generic(vec![0, 1, 0, 1, 0]).to_ratio_string());
286        assert_eq!("x^4 + x^2 + 1", Polynomial::from_vec_generic(vec![1, 0, 1, 0, 1]).to_ratio_string());
287        assert_eq!("-x^3 - x", Polynomial::from_vec_generic(vec![-1, 0, -1, 0]).to_ratio_string());
288        assert_eq!("-x^2 - 1", Polynomial::from_vec_generic(vec![0, -1, 0, -1]).to_ratio_string());
289        assert_eq!("-x^3 - x", Polynomial::from_vec_generic(vec![0, -1, 0, -1, 0]).to_ratio_string());
290        assert_eq!("-x^4 - x^2 - 1", Polynomial::from_vec_generic(vec![-1, 0, -1, 0, -1]).to_ratio_string());
291    }
292}
hmath/poly.rs

hmath/
poly.rs
