1use tiny_solver::na;
2
3pub const fn binomial_coefficient(n: usize, k: usize) -> usize {
4 if k > n {
5 0
6 } else {
7 let mut r = 1;
8 let mut d = 1;
9 let mut n = n;
10 while d < k + 1 {
11 r *= n;
12 n -= 1;
13 r /= d;
14 d += 1;
15 }
16 r
17 }
18}
19
20const fn power(num: f64, i: usize) -> f64 {
21 if i == 0 {
22 return 1.0;
23 }
24 let mut n = num;
25 let mut c = 1;
26 while c < i {
27 n *= num;
28 c += 1;
29 }
30 n
31}
32
33pub const fn compute_blending_matrix_c<const N: usize>(cumulative: bool) -> [[f64; N]; N] {
34 let mut m = [[0.0; N]; N];
35 let mut n = 2;
36 let mut factorial = 1.0;
37 while n < N {
38 factorial *= n as f64;
39 n += 1;
40 }
41 let mut i = 0;
42 while i < N {
43 let mut j = 0;
44 while j < N {
45 let mut s = j;
46 let mut sum = 0.0;
47 while s < N {
48 let sign = power(-1.0, s - j);
49 sum += sign
50 * binomial_coefficient(N, s - j) as f64
51 * power((N - s) as f64 - 1.0, N - 1 - i);
52 s += 1;
53 }
54 m[j][i] = binomial_coefficient(N - 1, N - 1 - i) as f64 * sum / factorial;
55 j += 1;
56 }
57 i += 1;
58 }
59 let mut i = 0;
60 if cumulative {
61 while i < N {
62 let mut j = i + 1;
63 while j < N {
64 let mut col = 0;
65 while col < N {
66 m[i][col] += m[j][col];
67 col += 1;
68 }
69 j += 1;
70 }
71 i += 1;
72 }
73 }
74 m
75}
76
77pub fn compute_blending_matrix<T: na::RealField>(n: usize, cumulative: bool) -> na::DMatrix<T> {
78 let mut m = na::DMatrix::zeros(n, n);
79 let factorial = (1..n).reduce(|acc, e| acc * e).unwrap() as f64;
80 for i in 0..n {
81 for j in 0..n {
82 let sum: f64 = (j..n)
83 .map(|s| {
84 let sign = (-1.0_f64).powi((s - j) as i32);
85 sign * binomial_coefficient(n, s - j) as f64
86 * ((n - s) as f64 - 1.0).powi((n - 1 - i) as i32)
87 })
88 .sum();
89 m[(j, i)] = binomial_coefficient(n - 1, n - 1 - i) as f64 * sum / factorial;
90 }
91 }
92 if cumulative {
93 for i in 0..n {
94 for j in i + 1..n {
95 m.set_row(i, &(m.row(i) + m.row(j)));
96 }
97 }
98 }
99 m.cast()
100}
101
102pub fn compute_base_coefficients(n: usize) -> na::DMatrix<f64> {
103 let mut base_coefficients = na::DMatrix::zeros(n, n);
104 base_coefficients.row_mut(0).add_scalar_mut(1.0);
105 let deg = n - 1;
106 let mut order = deg;
107 for row in 1..n {
108 for i in (deg - order)..n {
109 base_coefficients[(row, i)] =
110 (order + i - deg) as f64 * base_coefficients[(row - 1, i)];
111 }
112 order -= 1;
113 }
114 base_coefficients
115}
116
117pub const fn compute_base_coefficients_c<const N: usize>() -> [[f64; N]; N] {
118 let mut base_coefficients = [[0.0; N]; N];
119 let mut col = 0;
120 while col < N {
121 base_coefficients[0][col] = 1.0;
122 col += 1;
123 }
124 let deg = N - 1;
125 let mut order = deg;
126
127 let mut row = 1;
128 while row < N {
129 let mut i = deg - order;
130 while i < N {
131 base_coefficients[row][i] = (order + i - deg) as f64 * base_coefficients[row - 1][i];
132 i += 1;
133 }
134 order -= 1;
135 row += 1;
136 }
137 base_coefficients
138}
139
140#[cfg(test)]
141mod tests {
142 use crate::common::{
143 binomial_coefficient, compute_base_coefficients, compute_base_coefficients_c,
144 compute_blending_matrix, compute_blending_matrix_c,
145 };
146
147 #[test]
148 fn test_binomial_coefficient() {
149 assert!(binomial_coefficient(1, 1) == 1);
150 assert!(binomial_coefficient(3, 2) == 3);
151 assert!(binomial_coefficient(4, 2) == 6);
152 assert!(binomial_coefficient(5, 2) == 10);
153 }
154 #[test]
155 fn test_compute_blending_matrix() {
156 const N: usize = 5;
157 let m0 = compute_blending_matrix::<f64>(N, false);
158 let m1 = compute_blending_matrix_c::<N>(false);
159 for i in 0..N {
160 for j in 0..N {
161 assert!((m0[(i, j)] - m1[i][j]).abs() < 1e-10);
162 }
163 }
164 }
165
166 #[test]
167 fn test_compute_base_coefficients() {
168 const N: usize = 5;
169 let m0 = compute_base_coefficients(N);
170 let m1 = compute_base_coefficients_c::<N>();
171 for i in 0..N {
172 for j in 0..N {
173 assert!((m0[(i, j)] - m1[i][j]).abs() < 1e-10);
174 }
175 }
176 }
177}