fixed_trig/
lookup.rs

1//! Implementations based on lookup tables.
2
3use fixed::types::I48F16;
4use std::convert::TryInto;
5use std::include_bytes;
6
7const SIN_TABLE: &[u8] = include_bytes!("tables/sin.table");
8const SIN_TABLE_BITS: u8 = SIN_TABLE[0];
9const SIN_ANGLE_BITS: u8 = 16;
10const QUADRANT_HIGH_MASK: i64 = 1 << (SIN_ANGLE_BITS - 1);
11const QUADRANT_LOW_MASK: i64 = 1 << (SIN_ANGLE_BITS - 2);
12const SIN_TABLE_OFFSET: u8 = SIN_ANGLE_BITS - 2;
13const SIN_INTERP_OFFSET: u8 = SIN_TABLE_OFFSET - SIN_TABLE_BITS;
14const SIN_TABLE_MASK: i64 = ((1 << SIN_TABLE_BITS) - 1) << SIN_INTERP_OFFSET;
15const SIN_INTERP_MASK: i64 = (1 << SIN_INTERP_OFFSET) - 1;
16
17fn lookup_table(table: &[u8], index: i64) -> i64 {
18    let i = (index * 2 + 1) as usize;
19    i64::from(u16::from_le_bytes(table[i..(i + 2)].try_into().unwrap()))
20}
21
22pub fn sin(angle: I48F16) -> I48F16 {
23    let two_pi = I48F16::from_num(2_f64 * std::f64::consts::PI);
24    let a = angle.rem_euclid(two_pi);
25    let an = (i64::from(a.to_bits()) * (1 << SIN_ANGLE_BITS)) / two_pi.to_bits();
26    let is_negative_quadrant: bool = (an & QUADRANT_HIGH_MASK) != 0;
27    let is_odd_quadrant: bool = (an & QUADRANT_LOW_MASK) == 0;
28    let mut index = (an & SIN_TABLE_MASK) >> SIN_INTERP_OFFSET;
29    if !is_odd_quadrant {
30        index = (1 << SIN_TABLE_BITS) - 1 - index;
31    }
32    let x1 = lookup_table(SIN_TABLE, index);
33    let x2 = lookup_table(SIN_TABLE, index + 1);
34    let interp = an & SIN_INTERP_MASK;
35    let approx = ((x2 - x1) * interp) >> SIN_INTERP_OFFSET;
36    let mut sine = if is_odd_quadrant {
37        x1 + approx
38    } else {
39        x2 - approx
40    };
41    if is_negative_quadrant {
42        sine = -sine;
43    }
44    I48F16::from_bits(sine)
45}
46
47pub fn cos(val: I48F16) -> I48F16 {
48    let cos_offset = I48F16::from_num((5_f64 * std::f64::consts::PI) / 2_f64);
49    sin(val + cos_offset)
50}
51
52const ASIN_TABLE: &[u8] = include_bytes!("tables/asin.table");
53const ASIN_TABLE_BITS: u8 = ASIN_TABLE[0];
54const ASIN_VAL_BITS: u8 = 16;
55const ASIN_TABLE_OFFSET: u8 = ASIN_VAL_BITS;
56const ASIN_INTERP_OFFSET: u8 = ASIN_TABLE_OFFSET - ASIN_TABLE_BITS;
57const ASIN_TABLE_MASK: i64 = ((1 << ASIN_TABLE_BITS) - 1) << ASIN_INTERP_OFFSET;
58const ASIN_INTERP_MASK: i64 = (1 << ASIN_INTERP_OFFSET) - 1;
59
60pub fn asin(mut val: I48F16) -> I48F16 {
61    let pi_2 = I48F16::from_num(std::f64::consts::FRAC_PI_2);
62    let is_neg = val < 0;
63    if is_neg {
64        val = -val;
65    }
66    if val > 1 {
67        panic!("invalid asin value");
68    } else if val == 1 {
69        return if is_neg { -pi_2 } else { pi_2 };
70    }
71    let val_bits = val.to_bits();
72    let index = (val_bits & ASIN_TABLE_MASK) >> ASIN_INTERP_OFFSET;
73    let x1 = lookup_table(ASIN_TABLE, index);
74    let x2 = lookup_table(ASIN_TABLE, index + 1);
75    let interp = val_bits & ASIN_INTERP_MASK;
76    let approx = ((x2 - x1) * interp) >> ASIN_INTERP_OFFSET;
77    let mut asine = x1 + approx;
78    if is_neg {
79        asine = -asine;
80    }
81    I48F16::from_bits(asine) * pi_2
82}
83
84pub fn acos(val: I48F16) -> I48F16 {
85    asin(-val) + I48F16::from_num(std::f64::consts::FRAC_PI_2)
86}
87
88#[cfg(test)]
89mod tests {
90    use super::*;
91
92    #[test]
93    fn test_sin() {
94        for i in -100..100 {
95            let fx = f64::from(i) * 0.1_f64;
96            let x: I48F16 = I48F16::from_num(fx);
97            let sin_x: f64 = sin(x).to_num();
98            let sin_fx = fx.sin();
99            let dsx = (sin_x - sin_fx).abs();
100            if dsx > 0.001 {
101                panic!("mismatch sin({}): {} {}", x, sin_x, sin_fx);
102            }
103        }
104    }
105
106    #[test]
107    fn test_sin_comprehensive() {
108        for i in 0..(1 << 20) {
109            let x = I48F16::from_bits(i);
110            let fx: f64 = x.to_num();
111            let sin_x: f64 = sin(x).to_num();
112            let sin_fx = fx.sin();
113            let dsx = (sin_x - sin_fx).abs();
114            if dsx > 0.001 {
115                panic!("mismatch sin({}): {} {}", x, sin_x, sin_fx);
116            }
117        }
118    }
119
120    #[test]
121    fn test_cos() {
122        for i in -100..100 {
123            let fx = f64::from(i) * 0.1_f64;
124            let x: I48F16 = I48F16::from_num(fx);
125            let cos_x: f64 = cos(x).to_num();
126            let cos_fx = fx.cos();
127            let dsx = (cos_x - cos_fx).abs();
128            if dsx > 0.001 {
129                panic!("mismatch cos({}): {} {}", x, cos_x, cos_fx);
130            }
131        }
132    }
133
134    #[test]
135    fn test_cos_comprehensive() {
136        for i in 0..(1 << 20) {
137            let x = I48F16::from_bits(i);
138            let fx: f64 = x.to_num();
139            let cos_x: f64 = cos(x).to_num();
140            let cos_fx = fx.cos();
141            let dsx = (cos_x - cos_fx).abs();
142            if dsx > 0.001 {
143                panic!("mismatch cos({}): {} {}", x, cos_x, cos_fx);
144            }
145        }
146    }
147
148    #[test]
149    fn test_asin_comprehensive() {
150        let one = 1 << 16;
151        for i in -one..(one + 1) {
152            let x = I48F16::from_bits(i);
153            let fx: f64 = x.to_num();
154            let asin_x: f64 = asin(x).to_num();
155            let asin_fx = fx.asin();
156            let dsx = (asin_x - asin_fx).abs();
157            // does not do well at values very close to 1 or -1
158            if dsx > 0.001 && ((fx > -0.999 && fx < 0.999) || dsx > 0.02) {
159                panic!("mismatch asin({}): {} {}", x, asin_x, asin_fx);
160            }
161        }
162    }
163
164    #[test]
165    fn test_acos_comprehensive() {
166        let one = 1 << 16;
167        for i in -one..(one + 1) {
168            let x = I48F16::from_bits(i);
169            let fx: f64 = x.to_num();
170            let acos_x: f64 = acos(x).to_num();
171            let acos_fx = fx.acos();
172            let dsx = (acos_x - acos_fx).abs();
173            // does not do well at values very close to 1 or -1
174            if dsx > 0.001 && ((fx > -0.999 && fx < 0.999) || dsx > 0.02) {
175                panic!("mismatch acos({}): {} {}", x, acos_x, acos_fx);
176            }
177        }
178    }
179}
fixed_trig/lookup.rs

fixed_trig/
lookup.rs
