phyz_math/
fp.rs

1//! Reproducible floating-point primitives.
2//!
3//! `+`, `-`, `*`, `/` and `sqrt` are correctly rounded by IEEE-754, so every
4//! conforming target computes them to the same bits. The transcendental
5//! functions are not: `sin`, `cos`, `atan2`, `pow` and `tanh` come from the
6//! platform's libm, and Apple's, glibc's and musl's implementations disagree
7//! in the last ulp for many inputs. That is enough. A 1-ulp difference in one
8//! `sin` inside `integrate_configuration` is amplified by a contact-rich
9//! rollout until the trajectories are visibly different, which is exactly the
10//! failure mode where you cannot tell chaos from a bug.
11//!
12//! So the rollout path does not call the platform's libm. It calls this
13//! module, which forwards to the [`libm`] crate — a pure-Rust port of musl's
14//! libm that computes the same bits from the same inputs on every target,
15//! because it is the same Rust source compiled with the same IEEE semantics
16//! (Rust never contracts to FMA or reassociates FP without explicit opt-in).
17//!
18//! # The `portable-fp` feature
19//!
20//! On by default. Turning it off routes every function here straight to `std`,
21//! which is faster (the platform libm is usually vectorized and tuned) and
22//! gives up cross-platform bit equality. Do that only if you have measured the
23//! difference and do not need to compare numbers across machines.
24//!
25//! # What this does *not* fix
26//!
27//! Nothing here makes a chaotic rollout stop being chaotic. It makes the
28//! *same* rollout on two machines take the same path through that chaos. See
29//! `docs/determinism.md` for the full contract.
30
31/// Sine.
32#[inline]
33pub fn sin(x: f64) -> f64 {
34    #[cfg(feature = "portable-fp")]
35    {
36        libm::sin(x)
37    }
38    #[cfg(not(feature = "portable-fp"))]
39    {
40        x.sin()
41    }
42}
43
44/// Cosine.
45#[inline]
46pub fn cos(x: f64) -> f64 {
47    #[cfg(feature = "portable-fp")]
48    {
49        libm::cos(x)
50    }
51    #[cfg(not(feature = "portable-fp"))]
52    {
53        x.cos()
54    }
55}
56
57/// `(sin(x), cos(x))`, matching [`f64::sin_cos`]'s tuple order.
58#[inline]
59pub fn sin_cos(x: f64) -> (f64, f64) {
60    #[cfg(feature = "portable-fp")]
61    {
62        libm::sincos(x)
63    }
64    #[cfg(not(feature = "portable-fp"))]
65    {
66        x.sin_cos()
67    }
68}
69
70/// Two-argument arctangent.
71#[inline]
72pub fn atan2(y: f64, x: f64) -> f64 {
73    #[cfg(feature = "portable-fp")]
74    {
75        libm::atan2(y, x)
76    }
77    #[cfg(not(feature = "portable-fp"))]
78    {
79        y.atan2(x)
80    }
81}
82
83/// Arccosine.
84#[inline]
85pub fn acos(x: f64) -> f64 {
86    #[cfg(feature = "portable-fp")]
87    {
88        libm::acos(x)
89    }
90    #[cfg(not(feature = "portable-fp"))]
91    {
92        x.acos()
93    }
94}
95
96/// `x^y` for real `y`.
97#[inline]
98pub fn pow(x: f64, y: f64) -> f64 {
99    #[cfg(feature = "portable-fp")]
100    {
101        libm::pow(x, y)
102    }
103    #[cfg(not(feature = "portable-fp"))]
104    {
105        x.powf(y)
106    }
107}
108
109/// Hyperbolic tangent.
110#[inline]
111pub fn tanh(x: f64) -> f64 {
112    #[cfg(feature = "portable-fp")]
113    {
114        libm::tanh(x)
115    }
116    #[cfg(not(feature = "portable-fp"))]
117    {
118        x.tanh()
119    }
120}
121
122/// Natural exponential.
123#[inline]
124pub fn exp(x: f64) -> f64 {
125    #[cfg(feature = "portable-fp")]
126    {
127        libm::exp(x)
128    }
129    #[cfg(not(feature = "portable-fp"))]
130    {
131        x.exp()
132    }
133}
134
135/// Natural logarithm.
136#[inline]
137pub fn ln(x: f64) -> f64 {
138    #[cfg(feature = "portable-fp")]
139    {
140        libm::log(x)
141    }
142    #[cfg(not(feature = "portable-fp"))]
143    {
144        x.ln()
145    }
146}
147
148#[cfg(test)]
149mod tests {
150    use super::*;
151
152    /// The portable path must still be *correct*, not merely consistent: musl
153    /// and the platform libm agree to within an ulp or so, and a bug in the
154    /// wiring (a swapped `sincos` tuple, say) would show up as a gross
155    /// mismatch rather than a last-digit one.
156    #[test]
157    fn agrees_with_std_to_a_few_ulp() {
158        let xs: [f64; 12] = [
159            -3.0,
160            -1.25,
161            -0.5,
162            -1e-9,
163            0.0,
164            1e-9,
165            0.5,
166            1.0,
167            1.25,
168            2.0,
169            3.0,
170            std::f64::consts::TAU,
171        ];
172        for x in xs {
173            let tol = 8.0 * f64::EPSILON * (1.0 + x.abs());
174            assert!((sin(x) - x.sin()).abs() <= tol, "sin({x})");
175            assert!((cos(x) - x.cos()).abs() <= tol, "cos({x})");
176            assert!((tanh(x) - x.tanh()).abs() <= tol, "tanh({x})");
177            let (s, c) = sin_cos(x);
178            assert_eq!((s, c), (sin(x), cos(x)), "sin_cos({x}) tuple order");
179            assert!((atan2(x, 1.5) - x.atan2(1.5)).abs() <= tol, "atan2({x})");
180            if x > 0.0 {
181                assert!((pow(x, 1.5) - x.powf(1.5)).abs() <= 8.0 * tol, "pow({x})");
182                assert!((ln(x) - x.ln()).abs() <= 8.0 * tol, "ln({x})");
183            }
184        }
185        assert!((acos(0.25) - 0.25f64.acos()).abs() <= 8.0 * f64::EPSILON);
186    }
187}
phyz_math/fp.rs

phyz_math/
fp.rs
