Skip to main content

steel_math/
trig.rs

1//! Port of Mojang `net.minecraft.util.Mth`'s float trig lookup tables.
2//!
3//! Carvers, structures, entity physics, and a lot of other gameplay code call
4//! `Mth.sin(float)` / `Mth.cos(float)`, not `Math.sin(double)`. Matching
5//! block-level determinism with vanilla requires using the same 65536-entry
6//! `float` table instead of the precise FPU sine.
7//!
8//! # Vanilla reference
9//! ```java
10//! private static final float[] SIN = new float[65536];
11//! static {
12//!     for (int i = 0; i < 65536; i++)
13//!         SIN[i] = (float) Math.sin(i * (Math.PI * 2.0 / 65536.0));
14//! }
15//! public static float sin(double i)  { return SIN[(int)((long)(i * 10430.378350470453) & 65535L)]; }
16//! public static float cos(double i)  { return SIN[(int)((long)(i * 10430.378350470453 + 16384.0) & 65535L)]; }
17//! ```
18//! Note that the multiplier `10430.378350470453` is exactly
19//! `65536 / (2π)` rounded to a double, and the table is indexed by the
20//! low 16 bits of the scaled angle.
21
22use 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);
26/// `65536 / (2π)` — Mojang's stored constant.
27const 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
32/// The `Mth.SIN` table (one entry per 2π / 65536 radians, float-valued).
33static SIN_TABLE: LazyLock<Box<[f32; TABLE_LEN]>> = LazyLock::new(|| {
34    // Box on the heap to keep the binary's .bss section small.
35    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
45/// Arcsine angles in radians for inputs from 0 to 1 in steps of 1/256.
46static 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
53/// Cosines of the angles in [`ASIN_TAB`]
54static COS_TAB: LazyLock<[f64; ATAN_TABLE_LEN]> =
55    LazyLock::new(|| array::from_fn(|index| ASIN_TAB[index].cos()));
56
57/// Vanilla `Mth.sin(double)` — returns the float from the 65536-entry sine
58/// table indexed by the angle's low-16-bit bucket.
59#[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/// Vanilla `Mth.cos(double)` — like [`sin`] but phase-shifted by +π/2 using
68/// the same table.
69#[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
77/// Approximates `1 / sqrt(value)` using a fast inverse square root algorithm.
78const 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/// Returns vanilla's approximate angle in radians for the vector `(x, y)`.
88#[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        // Index 0 is sin(0), should be 0.
149        // Index 16384 is sin(π/2), should be 1.
150        // Verify a few well-known table entries.
151        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        // cos(x) == sin(x + π/2) via the same table.
160        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}