peroxide/special/
lanczos.rs

1//! Lanczos approximation Coefficient generator
2
3use crate::statistics::ops::{double_factorial, factorial, C};
4use crate::structure::matrix::Matrix;
5use crate::traits::matrix::MatrixTrait;
6use crate::traits::pointer::{Oxide, RedoxCommon};
7use crate::util::non_macro::zeros;
8use crate::util::useful::sgn;
9use std::f64::consts::PI;
10
11const G: f64 = 5f64;
12const N: usize = 7;
13
14// Lanczos g=5, n=7
15const LG5N7: [f64; 7] = [
16    1.000000000189712,
17    76.18009172948503,
18    -86.50532032927205,
19    24.01409824118972,
20    -1.2317395783752254,
21    0.0012086577526594748,
22    -0.00000539702438713199,
23];
24
25pub fn ln_gamma_approx(z: f64) -> f64 {
26    let z = z - 1f64;
27    let base = z + G + 0.5;
28    let mut s = 0f64;
29    for i in 1..N {
30        s += LG5N7[i] / (z + i as f64);
31    }
32    s += LG5N7[0];
33    (2f64 * PI).sqrt().ln() + s.ln() - base + base.ln() * (z + 0.5)
34}
35
36pub fn gamma_approx(z: f64) -> f64 {
37    if z > 1f64 {
38        let z_int = z as usize;
39        if z - (z_int as f64) == 0f64 {
40            return factorial(z_int - 1) as f64;
41        }
42    }
43
44    if z < 0.5 {
45        PI / ((PI * z).sin() * gamma_approx(1f64 - z))
46    } else {
47        ln_gamma_approx(z).exp()
48    }
49}
50
51/// Lanczos Approximation Coefficient
52pub fn tlg1(g: f64, n: usize) -> Vec<f64> {
53    (lanczos_coeff(g, n - 1).ox() * g.exp() / (2f64 * PI).sqrt()).red()
54}
55
56fn lanczos_coeff(g: f64, n: usize) -> Vec<f64> {
57    let m = dr_gen(n) * b_gen(n) * (c_gen(n) * dc_gen(n));
58    let f = f_gen(g, n);
59    m * f
60}
61
62fn b_gen(n: usize) -> Matrix {
63    Matrix::from_index(
64        |i, j| {
65            if i == 0 {
66                1f64
67            } else if j >= i {
68                sgn(j - i) * C(i + j - 1, j - i) as f64
69            } else {
70                0f64
71            }
72        },
73        (n + 1, n + 1),
74    )
75}
76
77fn c_gen(n: usize) -> Matrix {
78    Matrix::from_index(
79        |i, j| {
80            if i == 0 && j == 0 {
81                0.5
82            } else if j > i {
83                0f64
84            } else {
85                sgn(i - j) * 4f64.powi(j as i32) * (i as f64) * (C(i + j, 2 * j) as f64)
86                    / (i + j) as f64
87            }
88        },
89        (n + 1, n + 1),
90    )
91}
92
93fn dc_gen(n: usize) -> Matrix {
94    let mut m = zeros(n + 1, n + 1);
95    m[(0, 0)] = 2f64;
96    for i in 1..n + 1 {
97        m[(i, i)] = 2f64 * double_factorial(2 * i - 1) as f64;
98    }
99    m
100}
101
102fn dr_gen(n: usize) -> Matrix {
103    let mut m = zeros(n + 1, n + 1);
104    m[(0, 0)] = 1f64;
105    for i in 1..n + 1 {
106        m[(i, i)] = -((i * C(2 * i - 1, i)) as f64);
107    }
108    m
109}
110
111fn f(g: f64, n: usize) -> f64 {
112    2f64.sqrt() * (n as f64 + 0.5).exp() / (2f64 * (n as f64 + g) + 1f64).powf(n as f64 + 0.5)
113}
114
115fn f_gen(g: f64, n: usize) -> Vec<f64> {
116    let mut v = vec![0f64; n + 1];
117    for i in 0..n + 1 {
118        v[i] = f(g, i);
119    }
120    v
121}
peroxide/special/lanczos.rs

peroxide/special/
lanczos.rs
