From 8b62a2bd8c1e44db401ed4ec04d4f9fbf7674e63 Mon Sep 17 00:00:00 2001 From: Michael Kirk Date: Fri, 26 Jan 2024 16:50:12 -0800 Subject: [PATCH 1/3] use libm remquo and copysign for a little speedup MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This makes our implement of sincosd look more like the cpp version. $ cargo bench --bench="*" -- --baseline=4891fc4 Finished bench [optimized] target(s) in 0.02s Running benches/geodesic_benchmark.rs (target/release/deps/geodesic_benchmark-cda172d0984b8505) direct (c wrapper)/default time: [24.046 µs 24.071 µs 24.099 µs] change: [-0.2131% +0.0210% +0.2593%] (p = 0.86 > 0.05) No change in performance detected. Found 9 outliers among 100 measurements (9.00%) 3 (3.00%) low severe 2 (2.00%) low mild 1 (1.00%) high mild 3 (3.00%) high severe direct (rust impl)/default time: [26.129 µs 26.168 µs 26.211 µs] change: [-5.5845% -5.3729% -5.1792%] (p = 0.00 < 0.05) Performance has improved. Found 9 outliers among 100 measurements (9.00%) 4 (4.00%) high mild 5 (5.00%) high severe inverse (c wrapper)/default time: [45.061 µs 45.141 µs 45.227 µs] change: [+0.3738% +0.6326% +0.8919%] (p = 0.00 < 0.05) Change within noise threshold. Found 7 outliers among 100 measurements (7.00%) 3 (3.00%) high mild 4 (4.00%) high severe inverse (rust impl)/default time: [67.739 µs 67.796 µs 67.865 µs] change: [-3.2442% -3.0580% -2.8479%] (p = 0.00 < 0.05) Performance has improved. Found 11 outliers among 100 measurements (11.00%) 1 (1.00%) low severe 5 (5.00%) low mild 2 (2.00%) high mild 3 (3.00%) high severe --- Cargo.toml | 1 + README.md | 10 ++++---- src/geomath.rs | 64 ++++++++++---------------------------------------- 3 files changed, 18 insertions(+), 57 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index e765674..375eb17 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -25,6 +25,7 @@ default = ["accurate"] [dependencies] accurate = { version = "0.3", optional = true, default-features = false } +libm = { version = "0.2.8", default-features = false } [dev-dependencies] approx = "0.5.1" diff --git a/README.md b/README.md index bdc9019..68f636e 100644 --- a/README.md +++ b/README.md @@ -87,16 +87,16 @@ Which produces output like: ```text direct (c wrapper)/default - time: [24.055 µs 24.085 µs 24.117 µs] + time: [24.046 µs 24.071 µs 24.099 µs] direct (rust impl)/default - time: [27.760 µs 27.810 µs 27.867 µs] + time: [26.129 µs 26.168 µs 26.211 µs] inverse (c wrapper)/default - time: [46.461 µs 47.435 µs 48.557 µs] + time: [45.061 µs 45.141 µs 45.227 µs] inverse (rust impl)/default - time: [70.488 µs 70.841 µs 71.356 µs] + time: [67.739 µs 67.796 µs 67.865 µs] ``` -Showing that, at least in this benchmark, the Rust implementation is 16-52% slower than the c bindings. +Showing that, at least in this benchmark, the Rust implementation is 10-50% slower than the c bindings. diff --git a/src/geomath.rs b/src/geomath.rs index 12bf174..4482f06 100644 --- a/src/geomath.rs +++ b/src/geomath.rs @@ -136,66 +136,26 @@ pub fn ang_diff(x: f64, y: f64) -> (f64, f64) { } } -pub fn fmod(x: f64, y: f64) -> f64 { - x % y -} - /// Compute sine and cosine of x in degrees pub fn sincosd(x: f64) -> (f64, f64) { - // r = math.fmod(x, 360) if Math.isfinite(x) else Math.nan - let mut r = if x.is_finite() { - fmod(x, 360.0) - } else { - std::f64::NAN - }; - - // q = 0 if Math.isnan(r) else int(round(r / 90)) - let mut q = if r.is_nan() { - 0 - } else { - (r / 90.0).round() as i32 - }; + let (mut r, q) = libm::remquo(x, 90.0); - // r -= 90 * q; r = math.radians(r) - r -= 90.0 * q as f64; r = r.to_radians(); - // s = math.sin(r); c = math.cos(r) - let s = r.sin(); - let c = r.cos(); + let (mut sinx, mut cosx) = r.sin_cos(); - // q = q % 4 - q %= 4; - - // if q == 1: - // s, c = c, -s - // elif q == 2: - // s, c = -s, -c - // elif q == 3: - // s, c = -c, s - - let q = if q < 0 { q + 4 } else { q }; - - let (s, c) = if q == 1 { - (c, -s) - } else if q == 2 { - (-s, -c) - } else if q == 3 { - (-c, s) - } else { - debug_assert_eq!(q, 0); - (s, c) + (sinx, cosx) = match q as u32 & 3 { + 0 => (sinx, cosx), + 1 => (cosx, -sinx), + 2 => (-sinx, -cosx), + 3 => (-cosx, sinx), + _ => unreachable!(), }; - // # Remove the minus sign on -0.0 except for sin(-0.0). - // # On Windows 32-bit with python 2.7, math.fmod(-0.0, 360) = +0.0 - // # (x, c) here fixes this bug. See also Math::sincosd in the C++ library. - // # AngNormalize has a similar fix. - // s, c = (x, c) if x == 0 else (0.0+s, 0.0+c) - // return s, c - let (s, c) = if x == 0.0 { (x, c) } else { (0.0 + s, 0.0 + c) }; - - (s, c) + if sinx == 0.0 { + sinx = libm::copysign(sinx, x); // special values from F.10.1.13 + } + (sinx, cosx) } // Compute atan2(y, x) with result in degrees From ac380b1bfedcd92e682d3e277761499d7552c66f Mon Sep 17 00:00:00 2001 From: Michael Kirk Date: Tue, 6 Feb 2024 13:52:58 -0800 Subject: [PATCH 2/3] use built in `copysign` --- src/geomath.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/geomath.rs b/src/geomath.rs index 4482f06..2358ef8 100644 --- a/src/geomath.rs +++ b/src/geomath.rs @@ -153,7 +153,7 @@ pub fn sincosd(x: f64) -> (f64, f64) { }; if sinx == 0.0 { - sinx = libm::copysign(sinx, x); // special values from F.10.1.13 + sinx = sinx.copysign(x); // special values from F.10.1.13 } (sinx, cosx) } From a2c65f4e7c9933792dff843dc4ee28db92395aab Mon Sep 17 00:00:00 2001 From: Michael Kirk Date: Tue, 6 Feb 2024 14:14:20 -0800 Subject: [PATCH 3/3] fixup! use libm remquo and copysign for a little speedup Without the change, we see sign failures like: cosx(-810) = 0, but got -0 --- src/geomath.rs | 61 +++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 60 insertions(+), 1 deletion(-) diff --git a/src/geomath.rs b/src/geomath.rs index 2358ef8..267541b 100644 --- a/src/geomath.rs +++ b/src/geomath.rs @@ -152,8 +152,12 @@ pub fn sincosd(x: f64) -> (f64, f64) { _ => unreachable!(), }; + // special values from F.10.1.12 + cosx += 0.0; + + // special values from F.10.1.13 if sinx == 0.0 { - sinx = sinx.copysign(x); // special values from F.10.1.13 + sinx = sinx.copysign(x); } (sinx, cosx) } @@ -492,4 +496,59 @@ mod tests { 4.124513511893872e-05 ); } + + // corresponding to tests/signtest.cpp + mod sign_test { + use super::*; + fn is_equiv(x: f64, y: f64) -> bool { + (x.is_nan() && y.is_nan()) || (x == y && x.is_sign_positive() == y.is_sign_positive()) + } + + macro_rules! check_sincosd { + ($x: expr, $expected_sin: expr, $expected_cos: expr) => { + let (sinx, cosx) = sincosd($x); + assert!( + is_equiv(sinx, $expected_sin), + "sinx({}) = {}, but got {}", + $x, + $expected_sin, + sinx + ); + assert!( + is_equiv(cosx, $expected_cos), + "cosx({}) = {}, but got {}", + $x, + $expected_cos, + cosx + ); + }; + } + + #[test] + fn sin_cosd() { + check_sincosd!(f64::NEG_INFINITY, f64::NAN, f64::NAN); + check_sincosd!(-810.0, -1.0, 0.0); + check_sincosd!(-720.0, -0.0, 1.0); + check_sincosd!(-630.0, 1.0, 0.0); + check_sincosd!(-540.0, -0.0, -1.0); + check_sincosd!(-450.0, -1.0, 0.0); + check_sincosd!(-360.0, -0.0, 1.0); + check_sincosd!(-270.0, 1.0, 0.0); + check_sincosd!(-180.0, -0.0, -1.0); + check_sincosd!(-90.0, -1.0, 0.0); + check_sincosd!(-0.0, -0.0, 1.0); + check_sincosd!(0.0, 0.0, 1.0); + check_sincosd!(90.0, 1.0, 0.0); + check_sincosd!(180.0, 0.0, -1.0); + check_sincosd!(270.0, -1.0, 0.0); + check_sincosd!(360.0, 0.0, 1.0); + check_sincosd!(450.0, 1.0, 0.0); + check_sincosd!(540.0, 0.0, -1.0); + check_sincosd!(630.0, -1.0, 0.0); + check_sincosd!(720.0, 0.0, 1.0); + check_sincosd!(810.0, 1.0, 0.0); + check_sincosd!(f64::INFINITY, f64::NAN, f64::NAN); + check_sincosd!(f64::NAN, f64::NAN, f64::NAN); + } + } }