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
#[inline]
pub fn interpolate(a: f64, b: f64, fraction: f64) -> f64 {
a * (1.0 - fraction) + b * fraction
}
pub struct InterpolationTable<F: Fn(f64) -> f64> {
inner: Vec<f64>,
func: F,
offset: usize,
min_x: f64,
max_x: f64,
shift: f64,
}
impl<F: Fn(f64) -> f64> InterpolationTable<F> {
pub fn new(min_x: f64, max_x: f64, frac_digits: i32, func: F) -> Self {
let shift = 10.0_f64.powi(frac_digits);
let offset = (min_x * shift) as usize;
let mut table = InterpolationTable {
inner: Vec::new(),
func,
min_x,
max_x,
shift,
offset,
};
let mut i = table.index(min_x);
while i < table.index(max_x) {
let x = i as f64 / shift;
table.inner.push((table.func)(x));
i += 1;
}
table
}
#[inline]
fn index(&self, x: f64) -> usize {
(x * self.shift) as usize - self.offset
}
pub fn get(&self, x: f64) -> f64 {
if x < self.min_x || x >= self.max_x {
(self.func)(x)
} else {
let i = self.index(x);
let fraction = (x * self.shift - i as f64) / self.shift;
interpolate(self.inner[i], self.inner[i + 1], fraction)
}
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_interpolation_table() {
let table = InterpolationTable::new(0.0, 10.0, 5, |x| x.ln_1p());
for &x in &[0.02, 0.04, 0.45678686, 0.23875, 1.45345e-6] {
assert_relative_eq!(table.get(x), x.ln_1p(), epsilon = 0.00001);
}
}
}