zenith_float_num/
orthopoly.rs

1//! Classical orthogonal polynomials on [`ExactNum`] via three-term recurrences.
2
3use crate::defs::WORD_BIT_SIZE;
4use crate::Error;
5use crate::ExactNum;
6use crate::RoundingMode;
7
8/// Maximum degree for every family in this module. Larger `n` is `NaN`.
9pub const ORTHOPOLY_N_MAX: usize = 256;
10
11fn work_p(p: usize) -> usize {
12    p.saturating_add(WORD_BIT_SIZE)
13}
14
15fn op_nan() -> ExactNum {
16    ExactNum::nan(Some(Error::InvalidArgument))
17}
18
19fn finite(x: &ExactNum) -> bool {
20    !x.is_nan() && !x.is_inf()
21}
22
23fn n_ok(n: usize) -> bool {
24    n <= ORTHOPOLY_N_MAX
25}
26
27fn finish(mut y: ExactNum, p: usize, rm: RoundingMode) -> ExactNum {
28    let _ = y.set_precision(p, rm);
29    y
30}
31
32impl ExactNum {
33    /// Probabilist's Hermite polynomial `He_n(self)`.
34    ///
35    /// `He_0 = 1`, `He_1 = x`, `He_{n+1} = x He_n − n He_{n−1}`.
36    /// `n > ORTHOPOLY_N_MAX` or non-finite `self` is `NaN`.
37    pub fn hermite_he(&self, n: usize, p: usize, rm: RoundingMode) -> Self {
38        if !n_ok(n) || !finite(self) {
39            return op_nan();
40        }
41        let wrk = work_p(p);
42        if n == 0 {
43            return finish(ExactNum::from_u8(1, wrk), p, rm);
44        }
45        if n == 1 {
46            return finish(self.clone(), p, rm);
47        }
48        let mut prev = ExactNum::from_u8(1, wrk);
49        let mut cur = self.clone();
50        let _ = cur.set_precision(wrk, RoundingMode::None);
51        for k in 1..n {
52            let kf = ExactNum::from_u32(k as u32, wrk);
53            let next = self.mul(&cur, wrk, RoundingMode::None).sub(
54                &kf.mul(&prev, wrk, RoundingMode::None),
55                wrk,
56                RoundingMode::None,
57            );
58            prev = cur;
59            cur = next;
60        }
61        finish(cur, p, rm)
62    }
63
64    /// Physicist's Hermite polynomial `H_n(self)`.
65    ///
66    /// `H_0 = 1`, `H_1 = 2x`, `H_{n+1} = 2x H_n − 2n H_{n−1}`.
67    /// `n > ORTHOPOLY_N_MAX` or non-finite `self` is `NaN`.
68    pub fn hermite_h(&self, n: usize, p: usize, rm: RoundingMode) -> Self {
69        if !n_ok(n) || !finite(self) {
70            return op_nan();
71        }
72        let wrk = work_p(p);
73        let two = ExactNum::from_u8(2, wrk);
74        if n == 0 {
75            return finish(ExactNum::from_u8(1, wrk), p, rm);
76        }
77        if n == 1 {
78            return finish(two.mul(self, wrk, RoundingMode::None), p, rm);
79        }
80        let mut prev = ExactNum::from_u8(1, wrk);
81        let mut cur = two.mul(self, wrk, RoundingMode::None);
82        for k in 1..n {
83            let two_k = ExactNum::from_u32((2 * k) as u32, wrk);
84            let next = two
85                .mul(self, wrk, RoundingMode::None)
86                .mul(&cur, wrk, RoundingMode::None)
87                .sub(
88                    &two_k.mul(&prev, wrk, RoundingMode::None),
89                    wrk,
90                    RoundingMode::None,
91                );
92            prev = cur;
93            cur = next;
94        }
95        finish(cur, p, rm)
96    }
97
98    /// Laguerre polynomial `L_n(self)`.
99    ///
100    /// `L_0 = 1`, `L_1 = 1 − x`,
101    /// `L_{n+1} = ((2n+1−x) L_n − n L_{n−1}) / (n+1)`.
102    /// `n > ORTHOPOLY_N_MAX` or non-finite `self` is `NaN`.
103    pub fn laguerre(&self, n: usize, p: usize, rm: RoundingMode) -> Self {
104        self.gen_laguerre(n, &ExactNum::new(p), p, rm)
105    }
106
107    /// Generalized Laguerre `L_n^{(α)}(self)`.
108    ///
109    /// `L_0^{(α)} = 1`, `L_1^{(α)} = 1+α−x`,
110    /// `L_{n+1}^{(α)} = (((2n+1+α−x) L_n − (n+α) L_{n−1}) / (n+1)`.
111    /// `n > ORTHOPOLY_N_MAX` or a non-finite argument is `NaN`.
112    pub fn gen_laguerre(&self, n: usize, alpha: &Self, p: usize, rm: RoundingMode) -> Self {
113        if !n_ok(n) || !finite(self) || !finite(alpha) {
114            return op_nan();
115        }
116        let wrk = work_p(p);
117        let one = ExactNum::from_u8(1, wrk);
118        if n == 0 {
119            return finish(one, p, rm);
120        }
121        if n == 1 {
122            return finish(
123                one.add(alpha, wrk, RoundingMode::None)
124                    .sub(self, wrk, RoundingMode::None),
125                p,
126                rm,
127            );
128        }
129        let mut prev = one.clone();
130        let mut cur = one
131            .add(alpha, wrk, RoundingMode::None)
132            .sub(self, wrk, RoundingMode::None);
133        for k in 1..n {
134            let kf = ExactNum::from_u32(k as u32, wrk);
135            let two_k_1 = ExactNum::from_u32((2 * k + 1) as u32, wrk);
136            let coeff =
137                two_k_1
138                    .add(alpha, wrk, RoundingMode::None)
139                    .sub(self, wrk, RoundingMode::None);
140            let k_a = kf.add(alpha, wrk, RoundingMode::None);
141            let num = coeff.mul(&cur, wrk, RoundingMode::None).sub(
142                &k_a.mul(&prev, wrk, RoundingMode::None),
143                wrk,
144                RoundingMode::None,
145            );
146            let den = ExactNum::from_u32((k + 1) as u32, wrk);
147            let next = num.div(&den, wrk, RoundingMode::None);
148            prev = cur;
149            cur = next;
150        }
151        finish(cur, p, rm)
152    }
153
154    /// Chebyshev polynomial of the first kind `T_n(self)`.
155    ///
156    /// `T_0 = 1`, `T_1 = x`, `T_{n+1} = 2x T_n − T_{n−1}`.
157    /// `n > ORTHOPOLY_N_MAX` or non-finite `self` is `NaN`.
158    pub fn chebyshev_t(&self, n: usize, p: usize, rm: RoundingMode) -> Self {
159        if !n_ok(n) || !finite(self) {
160            return op_nan();
161        }
162        let wrk = work_p(p);
163        if n == 0 {
164            return finish(ExactNum::from_u8(1, wrk), p, rm);
165        }
166        if n == 1 {
167            return finish(self.clone(), p, rm);
168        }
169        let two = ExactNum::from_u8(2, wrk);
170        let mut prev = ExactNum::from_u8(1, wrk);
171        let mut cur = self.clone();
172        let _ = cur.set_precision(wrk, RoundingMode::None);
173        for _ in 1..n {
174            let next = two
175                .mul(self, wrk, RoundingMode::None)
176                .mul(&cur, wrk, RoundingMode::None)
177                .sub(&prev, wrk, RoundingMode::None);
178            prev = cur;
179            cur = next;
180        }
181        finish(cur, p, rm)
182    }
183
184    /// Chebyshev polynomial of the second kind `U_n(self)`.
185    ///
186    /// `U_0 = 1`, `U_1 = 2x`, `U_{n+1} = 2x U_n − U_{n−1}`.
187    /// `n > ORTHOPOLY_N_MAX` or non-finite `self` is `NaN`.
188    pub fn chebyshev_u(&self, n: usize, p: usize, rm: RoundingMode) -> Self {
189        if !n_ok(n) || !finite(self) {
190            return op_nan();
191        }
192        let wrk = work_p(p);
193        let two = ExactNum::from_u8(2, wrk);
194        if n == 0 {
195            return finish(ExactNum::from_u8(1, wrk), p, rm);
196        }
197        if n == 1 {
198            return finish(two.mul(self, wrk, RoundingMode::None), p, rm);
199        }
200        let mut prev = ExactNum::from_u8(1, wrk);
201        let mut cur = two.mul(self, wrk, RoundingMode::None);
202        for _ in 1..n {
203            let next = two
204                .mul(self, wrk, RoundingMode::None)
205                .mul(&cur, wrk, RoundingMode::None)
206                .sub(&prev, wrk, RoundingMode::None);
207            prev = cur;
208            cur = next;
209        }
210        finish(cur, p, rm)
211    }
212
213    /// Gegenbauer (ultraspherical) polynomial `C_n^{(λ)}(self)`.
214    ///
215    /// `C_0 = 1`, `C_1 = 2λ x`,
216    /// `C_{n+1} = (2(n+λ) x C_n − (n+2λ−1) C_{n−1}) / (n+1)`.
217    /// Standard `C_n^{(λ)}`: `C_2^{(1)} = 4x² − 1 = U_2`, and
218    /// `C_n^{(1/2)} = P_n` (Legendre). `n > ORTHOPOLY_N_MAX` or a
219    /// non-finite argument is `NaN`.
220    pub fn gegenbauer(&self, n: usize, lambda: &Self, p: usize, rm: RoundingMode) -> Self {
221        if !n_ok(n) || !finite(self) || !finite(lambda) {
222            return op_nan();
223        }
224        let wrk = work_p(p);
225        let one = ExactNum::from_u8(1, wrk);
226        if n == 0 {
227            return finish(one, p, rm);
228        }
229        let two = ExactNum::from_u8(2, wrk);
230        if n == 1 {
231            return finish(
232                two.mul(lambda, wrk, RoundingMode::None)
233                    .mul(self, wrk, RoundingMode::None),
234                p,
235                rm,
236            );
237        }
238        let mut prev = one.clone();
239        let mut cur = two
240            .mul(lambda, wrk, RoundingMode::None)
241            .mul(self, wrk, RoundingMode::None);
242        for k in 1..n {
243            let kf = ExactNum::from_u32(k as u32, wrk);
244            let two_n_l = two.mul(
245                &kf.add(lambda, wrk, RoundingMode::None),
246                wrk,
247                RoundingMode::None,
248            );
249            let n_2l_1 = kf
250                .add(
251                    &two.mul(lambda, wrk, RoundingMode::None),
252                    wrk,
253                    RoundingMode::None,
254                )
255                .sub(&one, wrk, RoundingMode::None);
256            let num = two_n_l
257                .mul(self, wrk, RoundingMode::None)
258                .mul(&cur, wrk, RoundingMode::None)
259                .sub(
260                    &n_2l_1.mul(&prev, wrk, RoundingMode::None),
261                    wrk,
262                    RoundingMode::None,
263                );
264            let den = ExactNum::from_u32((k + 1) as u32, wrk);
265            let next = num.div(&den, wrk, RoundingMode::None);
266            prev = cur;
267            cur = next;
268        }
269        finish(cur, p, rm)
270    }
271}
272
273#[cfg(test)]
274mod tests {
275    use super::*;
276    use crate::Consts;
277
278    fn gold_p() -> (usize, RoundingMode) {
279        (256, RoundingMode::ToEven)
280    }
281
282    #[test]
283    fn orthopoly_he4_laguerre_t5_gegenbauer_recurrence() {
284        let (p, rm) = gold_p();
285        let mut cc = Consts::new().expect("consts");
286        let zero = ExactNum::new(p);
287        let one = ExactNum::from_u8(1, p);
288        let two = ExactNum::from_u8(2, p);
289        let three = ExactNum::from_u8(3, p);
290        let half = one.div(&two, p, rm);
291
292        assert_eq!(zero.hermite_he(4, p, rm).cmp(&three), Some(0));
293        assert_eq!(zero.laguerre(3, p, rm).cmp(&one), Some(0));
294        assert_eq!(zero.gen_laguerre(3, &zero, p, rm).cmp(&one), Some(0));
295
296        let pi = cc.pi(p, rm);
297        let five = ExactNum::from_u8(5, p);
298        let c = pi.div(&five, p, rm).cos(p, rm, &mut cc);
299        let t5 = c.chebyshev_t(5, p, rm);
300        let neg_one = one.neg();
301        assert_eq!(t5.cmp(&neg_one), Some(0));
302
303        // Standard C_2^{(1)} = 4x² − 1 = U_2. The plan's "3x²−1 at λ=1"
304        // is twice Legendre P_2 = 2 C_2^{(1/2)}.
305        let g1 = half.gegenbauer(2, &one, p, rm);
306        let u2 = half.chebyshev_u(2, p, rm);
307        let four_x2_m1 = ExactNum::from_u8(4, p)
308            .mul(&half, p, rm)
309            .mul(&half, p, rm)
310            .sub(&one, p, rm);
311        assert_eq!(g1.cmp(&four_x2_m1), Some(0));
312        assert_eq!(g1.cmp(&u2), Some(0));
313        let lam_half = half.clone();
314        let g_leg = half.gegenbauer(2, &lam_half, p, rm);
315        let three_x2_m1 = three.mul(&half, p, rm).mul(&half, p, rm).sub(&one, p, rm);
316        let two_p2 = g_leg.mul(&two, p, rm);
317        assert_eq!(two_p2.cmp(&three_x2_m1), Some(0));
318
319        let t6 = half.chebyshev_t(6, p, rm);
320        let t5h = half.chebyshev_t(5, p, rm);
321        let t4h = half.chebyshev_t(4, p, rm);
322        let rec = two.mul(&half, p, rm).mul(&t5h, p, rm).sub(&t4h, p, rm);
323        assert_eq!(t6.cmp(&rec), Some(0));
324
325        assert!(zero.hermite_he(ORTHOPOLY_N_MAX + 1, p, rm).is_nan());
326        assert_eq!(
327            zero.hermite_h(4, p, rm).cmp(&ExactNum::from_u8(12, p)),
328            Some(0)
329        );
330    }
331}
zenith_float_num/orthopoly.rs

zenith_float_num/
orthopoly.rs
