diff --git a/Cargo.toml b/Cargo.toml
index d123fc89..83819aed 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -12,10 +12,12 @@ readme = "README.md"
 license = "Apache-2.0"
 
 [dependencies]
-libm = "0.2.2"
 serde = {version = "1.0.137", optional = true}
 regex = {version = "1.5.5", optional = true}
 serde_derive = {version = "1.0.137", optional = true}
+[dependencies.num-traits]
+version = "0.2"
+default-features = false
 
 [dev-dependencies]
 criterion = "0.3.5"
diff --git a/src/duration.rs b/src/duration.rs
index e478bfbf..dea25862 100644
--- a/src/duration.rs
+++ b/src/duration.rs
@@ -1,5 +1,3 @@
-use libm::{fabs, floor, pow};
-
 #[cfg(feature = "std")]
 use crate::ParsingErrors;
 use crate::{
@@ -493,14 +491,14 @@ impl Mul<f64> for Duration {
         let ten: f64 = 10.0;
 
         loop {
-            if fabs(floor(new_val) - new_val) < f64::EPSILON {
+            if (new_val.floor() - new_val).abs() < f64::EPSILON {
                 // Yay, we've found the precision of this number
                 break;
             }
             // Multiply by the precision
             // https://play.rust-lang.org/?version=stable&mode=debug&edition=2018&gist=b760579f103b7192c20413ebbe167b90
             p += 1;
-            new_val = q * pow(ten, f64::from(p));
+            new_val = q * ten.powi(p);
         }
 
         Duration::from_total_nanoseconds(
@@ -551,7 +549,7 @@ impl Mul<f64> for Unit {
             Unit::Microsecond => q * (NANOSECONDS_PER_MICROSECOND as f64),
             Unit::Nanosecond => q,
         };
-        if fabs(total_ns) < (i64::MAX as f64) {
+        if total_ns.abs() < (i64::MAX as f64) {
             Duration::from_truncated_nanoseconds(total_ns as i64)
         } else {
             Duration::from_total_nanoseconds(total_ns as i128)
@@ -581,7 +579,7 @@ macro_rules! impl_ops_for_type {
                     Freq::KiloHertz => NANOSECONDS_PER_MILLISECOND as f64 / (q as f64),
                     Freq::Hertz => (NANOSECONDS_PER_SECOND as f64) / (q as f64),
                 };
-                if fabs(total_ns) < (i64::MAX as f64) {
+                if total_ns.abs() < (i64::MAX as f64) {
                     Duration::from_truncated_nanoseconds(total_ns as i64)
                 } else {
                     Duration::from_total_nanoseconds(total_ns as i128)
@@ -655,7 +653,7 @@ impl fmt::LowerExp for Duration {
     // Prints the duration with appropriate units
     fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
         let seconds_f64 = self.in_seconds();
-        let seconds_f64_abs = fabs(seconds_f64);
+        let seconds_f64_abs = seconds_f64.abs();
         if seconds_f64_abs < 1e-5 {
             fmt::Display::fmt(&(seconds_f64 * 1e9), f)?;
             write!(f, " ns")
diff --git a/src/epoch.rs b/src/epoch.rs
index 44501fe6..09d5c539 100644
--- a/src/epoch.rs
+++ b/src/epoch.rs
@@ -1,5 +1,3 @@
-use libm::{fabs, sin, trunc};
-
 use crate::duration::{Duration, Unit};
 use crate::{
     Errors, TimeSystem, DAYS_GPS_TAI_OFFSET, ET_EPOCH_S, J1900_OFFSET, J2000_OFFSET, MJD_OFFSET,
@@ -303,9 +301,9 @@ impl Epoch {
         let g_rad = (PI / 180.0) * (357.528 + 35_999.050 * tt_centuries_j2k);
 
         // Decimal does not provide trig functions, so let's define the parts of the trig separately.
-        let inner = g_rad + 0.0167 * sin(g_rad);
+        let inner = g_rad + 0.0167 * g_rad.sin();
 
-        Self(tt_duration + ((ET_EPOCH_S as f64) - (0.001_658 * sin(inner))) * Unit::Second)
+        Self(tt_duration + ((ET_EPOCH_S as f64) - (0.001_658 * inner.sin())) * Unit::Second)
     }
 
     #[must_use]
@@ -815,7 +813,7 @@ impl Epoch {
         let inner = self.inner_g_rad();
 
         self.as_tt_duration() - (ET_EPOCH_S * Unit::Second)
-            + (0.001_658 * sin(inner)) * Unit::Second
+            + (0.001_658 * inner.sin()) * Unit::Second
     }
 
     #[must_use]
@@ -823,7 +821,7 @@ impl Epoch {
     pub fn as_tdb_seconds(&self) -> f64 {
         // Note that we redo the calculation of as_tdb_duration to save computational cost
         let inner = self.inner_g_rad();
-        self.as_tt_seconds() - (ET_EPOCH_S as f64) + (0.001_658 * sin(inner))
+        self.as_tt_seconds() - (ET_EPOCH_S as f64) + (0.001_658 * inner.sin())
     }
 
     /// For TDB computation, we're using f64 only because BigDecimal is far too slow for Nyx (uses FromStr).
@@ -831,7 +829,7 @@ impl Epoch {
         use core::f64::consts::PI;
         let g_rad = (PI / 180.0) * (357.528 + 35_999.050 * self.as_tt_centuries_j2k());
 
-        g_rad + 0.0167 * sin(g_rad)
+        g_rad + 0.0167 * g_rad.sin()
     }
 
     #[must_use]
@@ -853,7 +851,7 @@ impl Epoch {
     #[must_use]
     pub fn as_jde_tdb_duration(&self) -> Duration {
         let inner = self.inner_g_rad();
-        let tdb_delta = (0.001_658 * sin(inner)) * Unit::Second;
+        let tdb_delta = (0.001_658 * inner.sin()) * Unit::Second;
         self.as_jde_tt_duration() + tdb_delta
     }
 
@@ -1456,7 +1454,7 @@ fn div_rem_f64(me: f64, rhs: f64) -> (i32, f64) {
 }
 
 fn div_euclid_f64(lhs: f64, rhs: f64) -> f64 {
-    let q = trunc(lhs / rhs);
+    let q = (lhs / rhs).trunc();
     if lhs % rhs < 0.0 {
         return if rhs > 0.0 { q - 1.0 } else { q + 1.0 };
     }
@@ -1466,7 +1464,7 @@ fn div_euclid_f64(lhs: f64, rhs: f64) -> f64 {
 fn rem_euclid_f64(lhs: f64, rhs: f64) -> f64 {
     let r = lhs % rhs;
     if r < 0.0 {
-        r + fabs(rhs)
+        r + rhs.abs()
     } else {
         r
     }
diff --git a/src/lib.rs b/src/lib.rs
index 711c84b6..7faa04bf 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -71,7 +71,7 @@ pub mod prelude {
     pub use {Duration, Epoch, Freq, Frequencies, TimeSeries, TimeUnits, Unit};
 }
 
-extern crate libm;
+extern crate num_traits;
 
 #[cfg(feature = "std")]
 extern crate serde;