1use std::{array, f64::consts::PI, mem, sync::LazyLock};
23
24const ATAN_TABLE_LEN: usize = 257;
25const FRAC_BIAS: f64 = f64::from_bits(4_805_340_802_404_319_232);
26const INDEX_SCALE: f64 = 10_430.378_350_470_453;
28const ONE_SIXTH: f64 = 0.166_666_666_666_666_66;
29const TABLE_LEN: usize = 65_536;
30const TABLE_MASK: i64 = 0xFFFF;
31
32static SIN_TABLE: LazyLock<Box<[f32; TABLE_LEN]>> = LazyLock::new(|| {
34 let mut table: Box<[f32; TABLE_LEN]> = vec![0.0_f32; TABLE_LEN]
36 .into_boxed_slice()
37 .try_into()
38 .expect("65536-element vec");
39 for (i, slot) in table.iter_mut().enumerate() {
40 *slot = (i as f64 / INDEX_SCALE).sin() as f32;
41 }
42 table
43});
44
45static ASIN_TAB: LazyLock<[f64; ATAN_TABLE_LEN]> = LazyLock::new(|| {
47 array::from_fn(|index| {
48 let value = index as f64 / 256.0;
49 value.asin()
50 })
51});
52
53static COS_TAB: LazyLock<[f64; ATAN_TABLE_LEN]> =
55 LazyLock::new(|| array::from_fn(|index| ASIN_TAB[index].cos()));
56
57#[inline]
60#[must_use]
61pub fn sin(angle: impl Into<f64>) -> f32 {
62 let angle = angle.into();
63 let idx = (((angle * INDEX_SCALE) as i64) & TABLE_MASK) as usize;
64 SIN_TABLE[idx]
65}
66
67#[inline]
70#[must_use]
71pub fn cos(angle: impl Into<f64>) -> f32 {
72 let angle = angle.into();
73 let idx = (((angle * INDEX_SCALE + 16_384.0) as i64) & TABLE_MASK) as usize;
74 SIN_TABLE[idx]
75}
76
77const fn fast_inv_sqrt(value: f64) -> f64 {
79 let half = 0.5 * value;
80 let bits = value.to_bits() as i64;
81 let bits = 6_910_469_410_427_058_090_i64.wrapping_sub(bits >> 1);
82 let estimate = f64::from_bits(bits as u64);
83
84 estimate * (1.5 - half * estimate * estimate)
85}
86
87#[must_use]
89pub fn atan2(mut y: f64, mut x: f64) -> f64 {
90 let d2 = x * x + y * y;
91 if d2.is_nan() {
92 return f64::NAN;
93 }
94
95 let neg_y = y < 0.0;
96 if neg_y {
97 y = -y;
98 }
99
100 let neg_x = x < 0.0;
101 if neg_x {
102 x = -x;
103 }
104
105 let steep = y > x;
106 if steep {
107 mem::swap(&mut x, &mut y);
108 }
109
110 let rinv = fast_inv_sqrt(d2);
111 x *= rinv;
112 y *= rinv;
113
114 let yp = FRAC_BIAS + y;
115 let index = (yp.to_bits() as i32) as usize;
116 let phi = ASIN_TAB[index];
117 let c_phi = COS_TAB[index];
118 let s_phi = yp - FRAC_BIAS;
119 let sd = y * c_phi - x * s_phi;
120 let d = (6.0 + sd * sd) * sd * ONE_SIXTH;
121 let mut theta = phi + d;
122
123 if steep {
124 theta = (PI / 2.0) - theta;
125 }
126 if neg_x {
127 theta = PI - theta;
128 }
129 if neg_y {
130 theta = -theta;
131 }
132
133 theta
134}
135
136#[cfg(test)]
137mod test {
138 use super::*;
139
140 #[test]
141 fn sin_at_zero_is_zero() {
142 assert_eq!(sin(0.0), 0.0);
143 assert_eq!(cos(0.0), 1.0);
144 }
145
146 #[test]
147 fn sin_matches_table_entries() {
148 for &i in &[0usize, 1000, 16_384, 32_768, 49_152] {
152 let expected = (i as f64 / INDEX_SCALE).sin() as f32;
153 assert_eq!(SIN_TABLE[i], expected);
154 }
155 }
156
157 #[test]
158 fn cos_is_sin_shifted() {
159 let angle = 0.7;
161 let expected = {
162 let idx = (((angle * INDEX_SCALE + 16_384.0) as i64) & TABLE_MASK) as usize;
163 SIN_TABLE[idx]
164 };
165 assert_eq!(cos(angle), expected);
166 }
167}