qualia_core_db/solvers/special_functions/
orthogonal.rs

1//! Classical orthogonal polynomials by their three-term recurrences. Each evaluates
2//! `P_n(x)` in `O(n)` with no allocation.
3
4/// Legendre polynomial `P_n(x)`. `(n+1)P_{n+1} = (2n+1)x P_n − n P_{n-1}`.
5pub fn legendre(n: u32, x: f64) -> f64 {
6    if n == 0 {
7        return 1.0;
8    }
9    let (mut p0, mut p1) = (1.0, x);
10    for k in 1..n {
11        let kf = k as f64;
12        let p2 = ((2.0 * kf + 1.0) * x * p1 - kf * p0) / (kf + 1.0);
13        p0 = p1;
14        p1 = p2;
15    }
16    p1
17}
18
19/// Chebyshev polynomial of the first kind `T_n(x)`. `T_{n+1} = 2x T_n − T_{n-1}`.
20pub fn chebyshev_t(n: u32, x: f64) -> f64 {
21    if n == 0 {
22        return 1.0;
23    }
24    let (mut t0, mut t1) = (1.0, x);
25    for _ in 1..n {
26        let t2 = 2.0 * x * t1 - t0;
27        t0 = t1;
28        t1 = t2;
29    }
30    t1
31}
32
33/// Chebyshev polynomial of the second kind `U_n(x)`. `U_0 = 1`, `U_1 = 2x`,
34/// `U_{n+1} = 2x U_n − U_{n-1}`.
35pub fn chebyshev_u(n: u32, x: f64) -> f64 {
36    if n == 0 {
37        return 1.0;
38    }
39    let (mut u0, mut u1) = (1.0, 2.0 * x);
40    for _ in 1..n {
41        let u2 = 2.0 * x * u1 - u0;
42        u0 = u1;
43        u1 = u2;
44    }
45    u1
46}
47
48/// Physicists' Hermite polynomial `H_n(x)`. `H_0 = 1`, `H_1 = 2x`,
49/// `H_{n+1} = 2x H_n − 2n H_{n-1}`.
50pub fn hermite(n: u32, x: f64) -> f64 {
51    if n == 0 {
52        return 1.0;
53    }
54    let (mut h0, mut h1) = (1.0, 2.0 * x);
55    for k in 1..n {
56        let h2 = 2.0 * x * h1 - 2.0 * k as f64 * h0;
57        h0 = h1;
58        h1 = h2;
59    }
60    h1
61}
62
63/// Laguerre polynomial `L_n(x)`. `L_0 = 1`, `L_1 = 1 − x`,
64/// `(n+1)L_{n+1} = (2n+1−x)L_n − n L_{n-1}`.
65pub fn laguerre(n: u32, x: f64) -> f64 {
66    if n == 0 {
67        return 1.0;
68    }
69    let (mut l0, mut l1) = (1.0, 1.0 - x);
70    for k in 1..n {
71        let kf = k as f64;
72        let l2 = ((2.0 * kf + 1.0 - x) * l1 - kf * l0) / (kf + 1.0);
73        l0 = l1;
74        l1 = l2;
75    }
76    l1
77}
78
79#[cfg(test)]
80mod tests {
81    use super::*;
82    const EPS: f64 = 1e-12;
83
84    #[test]
85    fn legendre_known() {
86        // P_2(x) = (3x²−1)/2, P_3(x) = (5x³−3x)/2
87        assert!((legendre(2, 0.5) - (3.0 * 0.25 - 1.0) / 2.0).abs() < EPS);
88        assert!((legendre(3, 0.5) - (5.0 * 0.125 - 1.5) / 2.0).abs() < EPS);
89        assert!((legendre(5, 1.0) - 1.0).abs() < EPS); // P_n(1) = 1
90        assert!((legendre(0, 7.0) - 1.0).abs() < EPS);
91    }
92
93    #[test]
94    fn chebyshev_known() {
95        // T_n(cos θ) = cos(n θ)
96        let theta = 0.7_f64;
97        for n in 0..6u32 {
98            assert!((chebyshev_t(n, theta.cos()) - (n as f64 * theta).cos()).abs() < 1e-10);
99        }
100        // T_2 = 2x²−1, U_2 = 4x²−1
101        assert!((chebyshev_t(2, 0.3) - (2.0 * 0.09 - 1.0)).abs() < EPS);
102        assert!((chebyshev_u(2, 0.3) - (4.0 * 0.09 - 1.0)).abs() < EPS);
103    }
104
105    #[test]
106    fn hermite_known() {
107        // H_2 = 4x²−2, H_3 = 8x³−12x
108        assert!((hermite(2, 1.5) - (4.0 * 2.25 - 2.0)).abs() < EPS);
109        assert!((hermite(3, 1.5) - (8.0 * 3.375 - 18.0)).abs() < 1e-10);
110    }
111
112    #[test]
113    fn laguerre_known() {
114        // L_2 = (x²−4x+2)/2, L_n(0)=1
115        assert!((laguerre(2, 1.0) - (1.0 - 4.0 + 2.0) / 2.0).abs() < EPS);
116        assert!((laguerre(4, 0.0) - 1.0).abs() < EPS);
117    }
118}
qualia_core_db/solvers/special_functions/orthogonal.rs

qualia_core_db/solvers/special_functions/
orthogonal.rs
