1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
use super::polynomial::Polynomial;
use alga::general::ComplexField;
mod spline;
pub use spline::*;
pub fn lagrange<N: ComplexField>(
xs: &[N],
ys: &[N],
tol: N::RealField,
) -> Result<Polynomial<N>, String> {
if xs.len() != ys.len() {
return Err("lagrange: slices have mismatched dimension".to_owned());
}
let mut qs = vec![Polynomial::with_tolerance(tol)?; xs.len() * xs.len()];
for (ind, y) in ys.iter().enumerate() {
qs[ind] = polynomial![*y];
}
for i in 1..xs.len() {
let mut poly_2 = polynomial![N::one(), -xs[i]];
poly_2.set_tolerance(tol)?;
for j in 1..=i {
let mut poly_1 = polynomial![N::one(), -xs[i - j]];
poly_1.set_tolerance(tol)?;
let idenom = N::one() / (xs[i] - xs[i - j]);
let numer =
&poly_1 * &qs[i + xs.len() * (j - 1)] - &poly_2 * &qs[(i - 1) + xs.len() * (j - 1)];
qs[i + xs.len() * j] = numer * idenom;
}
}
for i in 0..=qs[xs.len() * xs.len() - 1].order() {
if qs[xs.len() * xs.len() - 1].get_coefficient(i).abs() < tol {
qs[xs.len() * xs.len() - 1].purge_coefficient(i);
}
}
qs[xs.len() * xs.len() - 1].purge_leading();
Ok(qs[xs.len() * xs.len() - 1].clone())
}
pub fn hermite<N: ComplexField>(
xs: &[N],
ys: &[N],
derivs: &[N],
tol: N::RealField,
) -> Result<Polynomial<N>, String> {
if xs.len() != ys.len() {
return Err("hermite: slices have mismatched dimension".to_owned());
}
if xs.len() != derivs.len() {
return Err("hermite: derivatives have mismatched dimension".to_owned());
}
let mut zs = vec![N::zero(); 2 * xs.len()];
let mut qs = vec![N::zero(); 4 * xs.len() * xs.len()];
for i in 0..xs.len() {
zs[2 * i] = xs[i];
zs[2 * i + 1] = xs[i];
qs[2 * i] = ys[i];
qs[2 * i + 1] = ys[i];
qs[2 * i + 1 + (2 * xs.len())] = derivs[i];
if i != 0 {
qs[2 * i + (2 * xs.len())] = (qs[2 * i] - qs[2 * i - 1]) / (zs[2 * i] - zs[2 * i - 1]);
}
}
for i in 2..2 * xs.len() {
for j in 2..=i {
qs[i + j * (2 * xs.len())] = (qs[i + (j - 1) * (2 * xs.len())]
- qs[i - 1 + (j - 1) * (2 * xs.len())])
/ (zs[i] - zs[i - j]);
}
}
let mut hermite = polynomial![N::zero()];
for i in (1..2 * xs.len()).rev() {
hermite += qs[i + i * (2 * xs.len())];
hermite *= polynomial![N::one(), -xs[(i - 1) / 2]]
}
hermite += qs[0];
for i in 0..=hermite.order() {
if hermite.get_coefficient(i).abs() < tol {
hermite.purge_coefficient(i);
}
}
hermite.purge_leading();
Ok(hermite)
}