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..267541b 100644
--- a/src/geomath.rs
+++ b/src/geomath.rs
@@ -136,66 +136,30 @@ 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
-    };
+    let (mut r, q) = libm::remquo(x, 90.0);
 
-    // 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
-    };
-
-    // 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();
-
-    // 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 (mut sinx, mut cosx) = r.sin_cos();
 
-    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) };
+    // special values from F.10.1.12
+    cosx += 0.0;
 
-    (s, c)
+    // special values from F.10.1.13
+    if sinx == 0.0 {
+        sinx = sinx.copysign(x);
+    }
+    (sinx, cosx)
 }
 
 // Compute atan2(y, x) with result in degrees
@@ -532,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);
+        }
+    }
 }