diff --git a/src/math/ceil.rs b/src/math/ceil.rs
index eda28b9a..22d89297 100644
--- a/src/math/ceil.rs
+++ b/src/math/ceil.rs
@@ -1,3 +1,4 @@
+#![allow(unreachable_code)]
 use core::f64;
 
 const TOINT: f64 = 1. / f64::EPSILON;
@@ -15,6 +16,24 @@ pub fn ceil(x: f64) -> f64 {
             return unsafe { ::core::intrinsics::ceilf64(x) }
         }
     }
+    #[cfg(all(target_arch = "x86", not(target_feature = "sse2")))]
+    {
+        //use an alternative implementation on x86, because the
+        //main implementation fails with the x87 FPU used by
+        //debian i386, probablly due to excess precision issues.
+        //basic implementation taken from https://github.com/rust-lang/libm/issues/219
+        use super::fabs;
+        if fabs(x).to_bits() < 4503599627370496.0_f64.to_bits() {
+            let truncated = x as i64 as f64;
+            if truncated < x {
+                return truncated + 1.0;
+            } else {
+                return truncated;
+            }
+        } else {
+            return x;
+        }
+    }
     let u: u64 = x.to_bits();
     let e: i64 = (u >> 52 & 0x7ff) as i64;
     let y: f64;
diff --git a/src/math/floor.rs b/src/math/floor.rs
index b2b76057..d09f9a1a 100644
--- a/src/math/floor.rs
+++ b/src/math/floor.rs
@@ -1,3 +1,4 @@
+#![allow(unreachable_code)]
 use core::f64;
 
 const TOINT: f64 = 1. / f64::EPSILON;
@@ -15,6 +16,24 @@ pub fn floor(x: f64) -> f64 {
             return unsafe { ::core::intrinsics::floorf64(x) }
         }
     }
+    #[cfg(all(target_arch = "x86", not(target_feature = "sse2")))]
+    {
+        //use an alternative implementation on x86, because the
+        //main implementation fails with the x87 FPU used by
+        //debian i386, probablly due to excess precision issues.
+        //basic implementation taken from https://github.com/rust-lang/libm/issues/219
+        use super::fabs;
+        if fabs(x).to_bits() < 4503599627370496.0_f64.to_bits() {
+            let truncated = x as i64 as f64;
+            if truncated > x {
+                return truncated - 1.0;
+            } else {
+                return truncated;
+            }
+        } else {
+            return x;
+        }
+    }
     let ui = x.to_bits();
     let e = ((ui >> 52) & 0x7ff) as i32;
 
diff --git a/src/math/j1f.rs b/src/math/j1f.rs
index 5095894d..17586f14 100644
--- a/src/math/j1f.rs
+++ b/src/math/j1f.rs
@@ -367,6 +367,10 @@ mod tests {
     }
     #[test]
     fn test_y1f_2002() {
-        assert_eq!(y1f(2.0000002_f32), -0.10703229_f32);
+        //allow slightly different result on x87
+        let res = y1f(2.0000002_f32);
+        if res != -0.10703231_f32 {
+            assert_eq!(res, -0.10703229_f32);
+        }
     }
 }
diff --git a/src/math/rem_pio2f.rs b/src/math/rem_pio2f.rs
index 5d392ba2..4d4499b9 100644
--- a/src/math/rem_pio2f.rs
+++ b/src/math/rem_pio2f.rs
@@ -43,7 +43,9 @@ pub(crate) fn rem_pio2f(x: f32) -> (i32, f64) {
     if ix < 0x4dc90fdb {
         /* |x| ~< 2^28*(pi/2), medium size */
         /* Use a specialized rint() to get fn.  Assume round-to-nearest. */
-        let f_n = x64 * INV_PIO2 + TOINT - TOINT;
+        // use to_bits and from_bits to force rounding to storage format on
+        // x87.
+        let f_n = f64::from_bits((x64 * INV_PIO2 + TOINT).to_bits()) - TOINT;
         return (f_n as i32, x64 - f_n * PIO2_1 - f_n * PIO2_1T);
     }
     if ix >= 0x7f800000 {
diff --git a/src/math/round.rs b/src/math/round.rs
index bf72f5b9..dbb36852 100644
--- a/src/math/round.rs
+++ b/src/math/round.rs
@@ -19,7 +19,9 @@ pub fn round(mut x: f64) -> f64 {
     if i >> 63 != 0 {
         x = -x;
     }
-    y = x + TOINT - TOINT - x;
+    // use to_bits and from_bits to force rounding to storage format on
+    // x87.
+    y = f64::from_bits((x + TOINT).to_bits()) - TOINT - x;
     if y > 0.5 {
         y = y + x - 1.0;
     } else if y <= -0.5 {
diff --git a/src/math/roundf.rs b/src/math/roundf.rs
index 497e88d6..1b1c0333 100644
--- a/src/math/roundf.rs
+++ b/src/math/roundf.rs
@@ -18,7 +18,9 @@ pub fn roundf(mut x: f32) -> f32 {
     if i >> 31 != 0 {
         x = -x;
     }
-    y = x + TOINT - TOINT - x;
+    // use to_bits and from_bits to force rounding to storage format on
+    // x87.
+    y = f32::from_bits((x + TOINT).to_bits()) - TOINT - x;
     if y > 0.5f32 {
         y = y + x - 1.0;
     } else if y <= -0.5f32 {
diff --git a/src/math/sincosf.rs b/src/math/sincosf.rs
index 2725caad..4e863236 100644
--- a/src/math/sincosf.rs
+++ b/src/math/sincosf.rs
@@ -145,10 +145,38 @@ mod tests {
             let (s_minus, c_minus) = sincosf(theta - 2. * PI);
 
             const TOLERANCE: f32 = 1e-6;
-            assert!((s - s_plus).abs() < TOLERANCE);
-            assert!((s - s_minus).abs() < TOLERANCE);
-            assert!((c - c_plus).abs() < TOLERANCE);
-            assert!((c - c_minus).abs() < TOLERANCE);
+            assert!(
+                (s - s_plus).abs() < TOLERANCE,
+                "|{} - {}| = {} >= {}",
+                s,
+                s_plus,
+                (s - s_plus).abs(),
+                TOLERANCE
+            );
+            assert!(
+                (s - s_minus).abs() < TOLERANCE,
+                "|{} - {}| = {} >= {}",
+                s,
+                s_minus,
+                (s - s_minus).abs(),
+                TOLERANCE
+            );
+            assert!(
+                (c - c_plus).abs() < TOLERANCE,
+                "|{} - {}| = {} >= {}",
+                c,
+                c_plus,
+                (c - c_plus).abs(),
+                TOLERANCE
+            );
+            assert!(
+                (c - c_minus).abs() < TOLERANCE,
+                "|{} - {}| = {} >= {}",
+                c,
+                c_minus,
+                (c - c_minus).abs(),
+                TOLERANCE
+            );
         }
     }
 }