use libm remquo and copysign for a little speedup

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

@@ -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"

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.

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