diff --git a/crates/bevy_math/clippy.toml b/crates/bevy_math/clippy.toml
new file mode 100644
index 0000000000..bad8801c01
--- /dev/null
+++ b/crates/bevy_math/clippy.toml
@@ -0,0 +1,29 @@
+disallowed-methods = [
+  { path = "f32::powi", reason = "use ops::FloatPow::squared, ops::FloatPow::cubed, or ops::powf instead for libm determinism" },
+  { path = "f32::log", reason = "use ops::ln, ops::log2, or ops::log10 instead for libm determinism" },
+  { path = "f32::abs_sub", reason = "deprecated and deeply confusing method" },
+  { path = "f32::powf", reason = "use ops::powf instead for libm determinism" },
+  { path = "f32::exp", reason = "use ops::exp instead for libm determinism" },
+  { path = "f32::exp2", reason = "use ops::exp2 instead for libm determinism" },
+  { path = "f32::ln", reason = "use ops::ln instead for libm determinism" },
+  { path = "f32::log2", reason = "use ops::log2 instead for libm determinism" },
+  { path = "f32::log10", reason = "use ops::log10 instead for libm determinism" },
+  { path = "f32::cbrt", reason = "use ops::cbrt instead for libm determinism" },
+  { path = "f32::hypot", reason = "use ops::hypot instead for libm determinism" },
+  { path = "f32::sin", reason = "use ops::sin instead for libm determinism" },
+  { path = "f32::cos", reason = "use ops::cos instead for libm determinism" },
+  { path = "f32::tan", reason = "use ops::tan instead for libm determinism" },
+  { path = "f32::asin", reason = "use ops::asin instead for libm determinism" },
+  { path = "f32::acos", reason = "use ops::acos instead for libm determinism" },
+  { path = "f32::atan", reason = "use ops::atan instead for libm determinism" },
+  { path = "f32::atan2", reason = "use ops::atan2 instead for libm determinism" },
+  { path = "f32::sin_cos", reason = "use ops::sin_cos instead for libm determinism" },
+  { path = "f32::exp_m1", reason = "use ops::exp_m1 instead for libm determinism" },
+  { path = "f32::ln_1p", reason = "use ops::ln_1p instead for libm determinism" },
+  { path = "f32::sinh", reason = "use ops::sinh instead for libm determinism" },
+  { path = "f32::cosh", reason = "use ops::cosh instead for libm determinism" },
+  { path = "f32::tanh", reason = "use ops::tanh instead for libm determinism" },
+  { path = "f32::asinh", reason = "use ops::asinh instead for libm determinism" },
+  { path = "f32::acosh", reason = "use ops::acosh instead for libm determinism" },
+  { path = "f32::atanh", reason = "use ops::atanh instead for libm determinism" },
+]
diff --git a/crates/bevy_math/src/bounding/bounded2d/mod.rs b/crates/bevy_math/src/bounding/bounded2d/mod.rs
index 631e5559d8..1fdf6c75b0 100644
--- a/crates/bevy_math/src/bounding/bounded2d/mod.rs
+++ b/crates/bevy_math/src/bounding/bounded2d/mod.rs
@@ -2,6 +2,7 @@ mod primitive_impls;
 
 use super::{BoundingVolume, IntersectsVolume};
 use crate::{
+    ops::FloatPow,
     prelude::{Mat2, Rot2, Vec2},
     Isometry2d,
 };
@@ -251,7 +252,7 @@ impl IntersectsVolume<BoundingCircle> for Aabb2d {
     fn intersects(&self, circle: &BoundingCircle) -> bool {
         let closest_point = self.closest_point(circle.center);
         let distance_squared = circle.center.distance_squared(closest_point);
-        let radius_squared = circle.radius().powi(2);
+        let radius_squared = circle.radius().squared();
         distance_squared <= radius_squared
     }
 }
@@ -261,7 +262,7 @@ mod aabb2d_tests {
     use super::Aabb2d;
     use crate::{
         bounding::{BoundingCircle, BoundingVolume, IntersectsVolume},
-        Vec2,
+        ops, Vec2,
     };
 
     #[test]
@@ -385,7 +386,7 @@ mod aabb2d_tests {
             max: Vec2::new(2.0, 2.0),
         };
         let transformed = a.transformed_by(Vec2::new(2.0, -2.0), std::f32::consts::FRAC_PI_4);
-        let half_length = 2_f32.hypot(2.0);
+        let half_length = ops::hypot(2.0, 2.0);
         assert_eq!(
             transformed.min,
             Vec2::new(2.0 - half_length, -half_length - 2.0)
@@ -535,7 +536,7 @@ impl BoundingVolume for BoundingCircle {
     #[inline(always)]
     fn contains(&self, other: &Self) -> bool {
         let diff = self.radius() - other.radius();
-        self.center.distance_squared(other.center) <= diff.powi(2).copysign(diff)
+        self.center.distance_squared(other.center) <= diff.squared().copysign(diff)
     }
 
     #[inline(always)]
@@ -593,7 +594,7 @@ impl IntersectsVolume<Self> for BoundingCircle {
     #[inline(always)]
     fn intersects(&self, other: &Self) -> bool {
         let center_distance_squared = self.center.distance_squared(other.center);
-        let radius_sum_squared = (self.radius() + other.radius()).powi(2);
+        let radius_sum_squared = (self.radius() + other.radius()).squared();
         center_distance_squared <= radius_sum_squared
     }
 }
diff --git a/crates/bevy_math/src/bounding/bounded2d/primitive_impls.rs b/crates/bevy_math/src/bounding/bounded2d/primitive_impls.rs
index 4452b87c69..5ed361982b 100644
--- a/crates/bevy_math/src/bounding/bounded2d/primitive_impls.rs
+++ b/crates/bevy_math/src/bounding/bounded2d/primitive_impls.rs
@@ -1,6 +1,7 @@
 //! Contains [`Bounded2d`] implementations for [geometric primitives](crate::primitives).
 
 use crate::{
+    ops,
     primitives::{
         Annulus, Arc2d, BoxedPolygon, BoxedPolyline2d, Capsule2d, Circle, CircularSector,
         CircularSegment, Ellipse, Line2d, Plane2d, Polygon, Polyline2d, Rectangle, RegularPolygon,
@@ -154,7 +155,7 @@ impl Bounded2d for Ellipse {
         let (ux, uy) = (hw * alpha_cos, hw * alpha_sin);
         let (vx, vy) = (hh * beta_cos, hh * beta_sin);
 
-        let half_size = Vec2::new(ux.hypot(vx), uy.hypot(vy));
+        let half_size = Vec2::new(ops::hypot(ux, vx), ops::hypot(uy, vy));
 
         Aabb2d::new(isometry.translation, half_size)
     }
@@ -403,6 +404,7 @@ mod tests {
 
     use crate::{
         bounding::Bounded2d,
+        ops::{self, FloatPow},
         primitives::{
             Annulus, Arc2d, Capsule2d, Circle, CircularSector, CircularSegment, Ellipse, Line2d,
             Plane2d, Polygon, Polyline2d, Rectangle, RegularPolygon, Rhombus, Segment2d,
@@ -505,7 +507,7 @@ mod tests {
                 // The exact coordinates here are not obvious, but can be computed by constructing
                 // an altitude from the midpoint of the chord to the y-axis and using the right triangle
                 // similarity theorem.
-                bounding_circle_center: Vec2::new(-apothem / 2.0, apothem.powi(2)),
+                bounding_circle_center: Vec2::new(-apothem / 2.0, apothem.squared()),
                 bounding_circle_radius: 0.5,
             },
             // Test case: handling of axis-aligned extrema
@@ -844,7 +846,7 @@ mod tests {
 
         let bounding_circle = segment.bounding_circle(isometry);
         assert_eq!(bounding_circle.center, translation);
-        assert_eq!(bounding_circle.radius(), 1.0_f32.hypot(0.5));
+        assert_eq!(bounding_circle.radius(), ops::hypot(1.0, 0.5));
     }
 
     #[test]
@@ -920,7 +922,7 @@ mod tests {
 
         let bounding_circle = rectangle.bounding_circle(Isometry2d::from_translation(translation));
         assert_eq!(bounding_circle.center, translation);
-        assert_eq!(bounding_circle.radius(), 1.0_f32.hypot(0.5));
+        assert_eq!(bounding_circle.radius(), ops::hypot(1.0, 0.5));
     }
 
     #[test]
diff --git a/crates/bevy_math/src/bounding/bounded3d/extrusion.rs b/crates/bevy_math/src/bounding/bounded3d/extrusion.rs
index b990b490f3..62171c0404 100644
--- a/crates/bevy_math/src/bounding/bounded3d/extrusion.rs
+++ b/crates/bevy_math/src/bounding/bounded3d/extrusion.rs
@@ -7,7 +7,7 @@ use crate::primitives::{
     BoxedPolygon, BoxedPolyline2d, Capsule2d, Cuboid, Cylinder, Ellipse, Extrusion, Line2d,
     Polygon, Polyline2d, Primitive2d, Rectangle, RegularPolygon, Segment2d, Triangle2d,
 };
-use crate::{Isometry2d, Isometry3d, Quat, Rot2};
+use crate::{ops, Isometry2d, Isometry3d, Quat, Rot2};
 
 use crate::{bounding::Bounded2d, primitives::Circle};
 
@@ -231,8 +231,8 @@ pub trait BoundedExtrusion: Primitive2d + Bounded2d {
             center,
             circle: Circle { radius },
         } = self.bounding_circle(Isometry2d::IDENTITY);
-        let radius = radius.hypot(half_depth);
-        let center = isometry.translation + isometry.rotation * Vec3A::from(center.extend(0.));
+        let radius = ops::hypot(radius, half_depth);
+        let center = isometry * Vec3A::from(center.extend(0.));
 
         BoundingSphere::new(center, radius)
     }
@@ -246,6 +246,7 @@ mod tests {
 
     use crate::{
         bounding::{Bounded3d, BoundingVolume},
+        ops,
         primitives::{
             Capsule2d, Circle, Ellipse, Extrusion, Line2d, Polygon, Polyline2d, Rectangle,
             RegularPolygon, Segment2d, Triangle2d,
@@ -265,7 +266,7 @@ mod tests {
 
         let bounding_sphere = cylinder.bounding_sphere(isometry);
         assert_eq!(bounding_sphere.center, translation.into());
-        assert_eq!(bounding_sphere.radius(), 1f32.hypot(0.5));
+        assert_eq!(bounding_sphere.radius(), ops::hypot(1.0, 0.5));
     }
 
     #[test]
diff --git a/crates/bevy_math/src/bounding/bounded3d/mod.rs b/crates/bevy_math/src/bounding/bounded3d/mod.rs
index fa92c3f256..60822fa3fe 100644
--- a/crates/bevy_math/src/bounding/bounded3d/mod.rs
+++ b/crates/bevy_math/src/bounding/bounded3d/mod.rs
@@ -4,7 +4,7 @@ mod primitive_impls;
 use glam::Mat3;
 
 use super::{BoundingVolume, IntersectsVolume};
-use crate::{Isometry3d, Quat, Vec3A};
+use crate::{ops::FloatPow, Isometry3d, Quat, Vec3A};
 
 #[cfg(feature = "bevy_reflect")]
 use bevy_reflect::Reflect;
@@ -253,7 +253,7 @@ impl IntersectsVolume<BoundingSphere> for Aabb3d {
     fn intersects(&self, sphere: &BoundingSphere) -> bool {
         let closest_point = self.closest_point(sphere.center);
         let distance_squared = sphere.center.distance_squared(closest_point);
-        let radius_squared = sphere.radius().powi(2);
+        let radius_squared = sphere.radius().squared();
         distance_squared <= radius_squared
     }
 }
@@ -263,7 +263,7 @@ mod aabb3d_tests {
     use super::Aabb3d;
     use crate::{
         bounding::{BoundingSphere, BoundingVolume, IntersectsVolume},
-        Quat, Vec3, Vec3A,
+        ops, Quat, Vec3, Vec3A,
     };
 
     #[test]
@@ -389,7 +389,7 @@ mod aabb3d_tests {
             Vec3A::new(2.0, -2.0, 4.0),
             Quat::from_rotation_z(std::f32::consts::FRAC_PI_4),
         );
-        let half_length = 2_f32.hypot(2.0);
+        let half_length = ops::hypot(2.0, 2.0);
         assert_eq!(
             transformed.min,
             Vec3A::new(2.0 - half_length, -half_length - 2.0, 2.0)
@@ -518,7 +518,7 @@ impl BoundingSphere {
         let radius = self.radius();
         let distance_squared = (point - self.center).length_squared();
 
-        if distance_squared <= radius.powi(2) {
+        if distance_squared <= radius.squared() {
             // The point is inside the sphere.
             point
         } else {
@@ -553,7 +553,7 @@ impl BoundingVolume for BoundingSphere {
     #[inline(always)]
     fn contains(&self, other: &Self) -> bool {
         let diff = self.radius() - other.radius();
-        self.center.distance_squared(other.center) <= diff.powi(2).copysign(diff)
+        self.center.distance_squared(other.center) <= diff.squared().copysign(diff)
     }
 
     #[inline(always)]
@@ -621,7 +621,7 @@ impl IntersectsVolume<Self> for BoundingSphere {
     #[inline(always)]
     fn intersects(&self, other: &Self) -> bool {
         let center_distance_squared = self.center.distance_squared(other.center);
-        let radius_sum_squared = (self.radius() + other.radius()).powi(2);
+        let radius_sum_squared = (self.radius() + other.radius()).squared();
         center_distance_squared <= radius_sum_squared
     }
 }
diff --git a/crates/bevy_math/src/bounding/bounded3d/primitive_impls.rs b/crates/bevy_math/src/bounding/bounded3d/primitive_impls.rs
index 0d984657ea..8827b0f253 100644
--- a/crates/bevy_math/src/bounding/bounded3d/primitive_impls.rs
+++ b/crates/bevy_math/src/bounding/bounded3d/primitive_impls.rs
@@ -4,6 +4,7 @@ use glam::Vec3A;
 
 use crate::{
     bounding::{Bounded2d, BoundingCircle},
+    ops,
     primitives::{
         BoxedPolyline3d, Capsule3d, Cone, ConicalFrustum, Cuboid, Cylinder, InfinitePlane3d,
         Line3d, Polyline3d, Segment3d, Sphere, Torus, Triangle2d, Triangle3d,
@@ -137,7 +138,7 @@ impl Bounded3d for Cylinder {
     }
 
     fn bounding_sphere(&self, isometry: Isometry3d) -> BoundingSphere {
-        let radius = self.radius.hypot(self.half_height);
+        let radius = ops::hypot(self.radius, self.half_height);
         BoundingSphere::new(isometry.translation, radius)
     }
 }
@@ -345,7 +346,7 @@ impl Bounded3d for Triangle3d {
 
 #[cfg(test)]
 mod tests {
-    use crate::{bounding::BoundingVolume, Isometry3d};
+    use crate::{bounding::BoundingVolume, ops, Isometry3d};
     use glam::{Quat, Vec3, Vec3A};
 
     use crate::{
@@ -441,7 +442,7 @@ mod tests {
 
         let bounding_sphere = segment.bounding_sphere(isometry);
         assert_eq!(bounding_sphere.center, translation.into());
-        assert_eq!(bounding_sphere.radius(), 1.0_f32.hypot(0.5));
+        assert_eq!(bounding_sphere.radius(), ops::hypot(1.0, 0.5));
     }
 
     #[test]
@@ -461,7 +462,10 @@ mod tests {
 
         let bounding_sphere = polyline.bounding_sphere(isometry);
         assert_eq!(bounding_sphere.center, translation.into());
-        assert_eq!(bounding_sphere.radius(), 1.0_f32.hypot(1.0).hypot(1.0));
+        assert_eq!(
+            bounding_sphere.radius(),
+            ops::hypot(ops::hypot(1.0, 1.0), 1.0)
+        );
     }
 
     #[test]
@@ -479,7 +483,10 @@ mod tests {
 
         let bounding_sphere = cuboid.bounding_sphere(Isometry3d::from_translation(translation));
         assert_eq!(bounding_sphere.center, translation.into());
-        assert_eq!(bounding_sphere.radius(), 1.0_f32.hypot(0.5).hypot(0.5));
+        assert_eq!(
+            bounding_sphere.radius(),
+            ops::hypot(ops::hypot(1.0, 0.5), 0.5)
+        );
     }
 
     #[test]
@@ -500,7 +507,7 @@ mod tests {
 
         let bounding_sphere = cylinder.bounding_sphere(isometry);
         assert_eq!(bounding_sphere.center, translation.into());
-        assert_eq!(bounding_sphere.radius(), 1.0_f32.hypot(0.5));
+        assert_eq!(bounding_sphere.radius(), ops::hypot(1.0, 0.5));
     }
 
     #[test]
diff --git a/crates/bevy_math/src/bounding/raycast2d.rs b/crates/bevy_math/src/bounding/raycast2d.rs
index e3a4764725..745e53de82 100644
--- a/crates/bevy_math/src/bounding/raycast2d.rs
+++ b/crates/bevy_math/src/bounding/raycast2d.rs
@@ -1,5 +1,5 @@
 use super::{Aabb2d, BoundingCircle, IntersectsVolume};
-use crate::{Dir2, Ray2d, Vec2};
+use crate::{ops::FloatPow, Dir2, Ray2d, Vec2};
 
 #[cfg(feature = "bevy_reflect")]
 use bevy_reflect::Reflect;
@@ -76,8 +76,8 @@ impl RayCast2d {
         let offset = self.ray.origin - circle.center;
         let projected = offset.dot(*self.ray.direction);
         let closest_point = offset - projected * *self.ray.direction;
-        let distance_squared = circle.radius().powi(2) - closest_point.length_squared();
-        if distance_squared < 0. || projected.powi(2).copysign(-projected) < -distance_squared {
+        let distance_squared = circle.radius().squared() - closest_point.length_squared();
+        if distance_squared < 0. || projected.squared().copysign(-projected) < -distance_squared {
             None
         } else {
             let toi = -projected - distance_squared.sqrt();
diff --git a/crates/bevy_math/src/bounding/raycast3d.rs b/crates/bevy_math/src/bounding/raycast3d.rs
index fa36a438d6..3a369d9fe1 100644
--- a/crates/bevy_math/src/bounding/raycast3d.rs
+++ b/crates/bevy_math/src/bounding/raycast3d.rs
@@ -1,5 +1,5 @@
 use super::{Aabb3d, BoundingSphere, IntersectsVolume};
-use crate::{Dir3A, Ray3d, Vec3A};
+use crate::{ops::FloatPow, Dir3A, Ray3d, Vec3A};
 
 #[cfg(feature = "bevy_reflect")]
 use bevy_reflect::Reflect;
@@ -73,8 +73,8 @@ impl RayCast3d {
         let offset = self.origin - sphere.center;
         let projected = offset.dot(*self.direction);
         let closest_point = offset - projected * *self.direction;
-        let distance_squared = sphere.radius().powi(2) - closest_point.length_squared();
-        if distance_squared < 0. || projected.powi(2).copysign(-projected) < -distance_squared {
+        let distance_squared = sphere.radius().squared() - closest_point.length_squared();
+        if distance_squared < 0. || projected.squared().copysign(-projected) < -distance_squared {
             None
         } else {
             let toi = -projected - distance_squared.sqrt();
diff --git a/crates/bevy_math/src/common_traits.rs b/crates/bevy_math/src/common_traits.rs
index cc8869f627..525d8e4035 100644
--- a/crates/bevy_math/src/common_traits.rs
+++ b/crates/bevy_math/src/common_traits.rs
@@ -1,6 +1,6 @@
 //! This module contains abstract mathematical traits shared by types used in `bevy_math`.
 
-use crate::{Dir2, Dir3, Dir3A, Quat, Rot2, Vec2, Vec3, Vec3A, Vec4};
+use crate::{ops, Dir2, Dir3, Dir3A, Quat, Rot2, Vec2, Vec3, Vec3A, Vec4};
 use std::fmt::Debug;
 use std::ops::{Add, Div, Mul, Neg, Sub};
 
@@ -256,7 +256,7 @@ pub trait StableInterpolate: Clone {
     /// object_position.smooth_nudge(&target_position, decay_rate, delta_time);
     /// ```
     fn smooth_nudge(&mut self, target: &Self, decay_rate: f32, delta: f32) {
-        self.interpolate_stable_assign(target, 1.0 - f32::exp(-decay_rate * delta));
+        self.interpolate_stable_assign(target, 1.0 - ops::exp(-decay_rate * delta));
     }
 }
 
diff --git a/crates/bevy_math/src/cubic_splines.rs b/crates/bevy_math/src/cubic_splines.rs
index 3e4bcc0e9a..b99c035c8e 100644
--- a/crates/bevy_math/src/cubic_splines.rs
+++ b/crates/bevy_math/src/cubic_splines.rs
@@ -2,7 +2,7 @@
 
 use std::{fmt::Debug, iter::once};
 
-use crate::{Vec2, VectorSpace};
+use crate::{ops::FloatPow, Vec2, VectorSpace};
 
 use itertools::Itertools;
 use thiserror::Error;
@@ -731,12 +731,12 @@ impl<P: VectorSpace> CubicNurbs<P> {
         // t[5] := t_i+2
         // t[6] := t_i+3
 
-        let m00 = (t[4] - t[3]).powi(2) / ((t[4] - t[2]) * (t[4] - t[1]));
-        let m02 = (t[3] - t[2]).powi(2) / ((t[5] - t[2]) * (t[4] - t[2]));
+        let m00 = (t[4] - t[3]).squared() / ((t[4] - t[2]) * (t[4] - t[1]));
+        let m02 = (t[3] - t[2]).squared() / ((t[5] - t[2]) * (t[4] - t[2]));
         let m12 = (3.0 * (t[4] - t[3]) * (t[3] - t[2])) / ((t[5] - t[2]) * (t[4] - t[2]));
-        let m22 = 3.0 * (t[4] - t[3]).powi(2) / ((t[5] - t[2]) * (t[4] - t[2]));
-        let m33 = (t[4] - t[3]).powi(2) / ((t[6] - t[3]) * (t[5] - t[3]));
-        let m32 = -m22 / 3.0 - m33 - (t[4] - t[3]).powi(2) / ((t[5] - t[3]) * (t[5] - t[2]));
+        let m22 = 3.0 * (t[4] - t[3]).squared() / ((t[5] - t[2]) * (t[4] - t[2]));
+        let m33 = (t[4] - t[3]).squared() / ((t[6] - t[3]) * (t[5] - t[3]));
+        let m32 = -m22 / 3.0 - m33 - (t[4] - t[3]).squared() / ((t[5] - t[3]) * (t[5] - t[2]));
         [
             [m00, 1.0 - m00 - m02, m02, 0.0],
             [-3.0 * m00, 3.0 * m00 - m12, m12, 0.0],
@@ -1254,7 +1254,7 @@ impl<P: VectorSpace> RationalSegment<P> {
         // Position = N/D therefore
         // Velocity = (N/D)' = N'/D - N * D'/D^2 = (N' * D - N * D')/D^2
         numerator_derivative / denominator
-            - numerator * (denominator_derivative / denominator.powi(2))
+            - numerator * (denominator_derivative / denominator.squared())
     }
 
     /// Instantaneous acceleration of a point at parametric value `t` in `[0, knot_span)`.
@@ -1288,10 +1288,10 @@ impl<P: VectorSpace> RationalSegment<P> {
         // Velocity = (N/D)' = N'/D - N * D'/D^2 = (N' * D - N * D')/D^2
         // Acceleration = (N/D)'' = ((N' * D - N * D')/D^2)' = N''/D + N' * (-2D'/D^2) + N * (-D''/D^2 + 2D'^2/D^3)
         numerator_second_derivative / denominator
-            + numerator_derivative * (-2.0 * denominator_derivative / denominator.powi(2))
+            + numerator_derivative * (-2.0 * denominator_derivative / denominator.squared())
             + numerator
-                * (-denominator_second_derivative / denominator.powi(2)
-                    + 2.0 * denominator_derivative.powi(2) / denominator.powi(3))
+                * (-denominator_second_derivative / denominator.squared()
+                    + 2.0 * denominator_derivative.squared() / denominator.cubed())
     }
 
     /// Calculate polynomial coefficients for the cubic polynomials using a characteristic matrix.
@@ -1512,9 +1512,12 @@ impl<P: VectorSpace> From<CubicCurve<P>> for RationalCurve<P> {
 mod tests {
     use glam::{vec2, Vec2};
 
-    use crate::cubic_splines::{
-        CubicBSpline, CubicBezier, CubicGenerator, CubicNurbs, CubicSegment, RationalCurve,
-        RationalGenerator,
+    use crate::{
+        cubic_splines::{
+            CubicBSpline, CubicBezier, CubicGenerator, CubicNurbs, CubicSegment, RationalCurve,
+            RationalGenerator,
+        },
+        ops::{self, FloatPow},
     };
 
     /// How close two floats can be and still be considered equal
@@ -1540,10 +1543,10 @@ mod tests {
     /// Manual, hardcoded function for computing the position along a cubic bezier.
     fn cubic_manual(t: f32, points: [Vec2; 4]) -> Vec2 {
         let p = points;
-        p[0] * (1.0 - t).powi(3)
-            + 3.0 * p[1] * t * (1.0 - t).powi(2)
-            + 3.0 * p[2] * t.powi(2) * (1.0 - t)
-            + p[3] * t.powi(3)
+        p[0] * (1.0 - t).cubed()
+            + 3.0 * p[1] * t * (1.0 - t).squared()
+            + 3.0 * p[2] * t.squared() * (1.0 - t)
+            + p[3] * t.cubed()
     }
 
     /// Basic cubic Bezier easing test to verify the shape of the curve.
@@ -1674,8 +1677,8 @@ mod tests {
         // subjecting ones self to a lot of tedious matrix algebra.
 
         let alpha = FRAC_PI_2;
-        let leg = 2.0 * f32::sin(alpha / 2.0) / (1.0 + 2.0 * f32::cos(alpha / 2.0));
-        let weight = (1.0 + 2.0 * f32::cos(alpha / 2.0)) / 3.0;
+        let leg = 2.0 * ops::sin(alpha / 2.0) / (1.0 + 2.0 * ops::cos(alpha / 2.0));
+        let weight = (1.0 + 2.0 * ops::cos(alpha / 2.0)) / 3.0;
         let points = [
             vec2(1.0, 0.0),
             vec2(1.0, leg),
diff --git a/crates/bevy_math/src/curve/mod.rs b/crates/bevy_math/src/curve/mod.rs
index 30d5fee9b3..54d4360a73 100644
--- a/crates/bevy_math/src/curve/mod.rs
+++ b/crates/bevy_math/src/curve/mod.rs
@@ -596,7 +596,7 @@ where
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::Quat;
+    use crate::{ops, Quat};
     use approx::{assert_abs_diff_eq, AbsDiffEq};
     use std::f32::consts::TAU;
 
@@ -616,8 +616,8 @@ mod tests {
         assert!(curve.sample_unchecked(2.0).abs_diff_eq(&4.0, f32::EPSILON));
         assert!(curve.sample_unchecked(-3.0).abs_diff_eq(&9.0, f32::EPSILON));
 
-        let curve = function_curve(interval(0.0, f32::INFINITY).unwrap(), f32::log2);
-        assert_eq!(curve.sample_unchecked(3.5), f32::log2(3.5));
+        let curve = function_curve(interval(0.0, f32::INFINITY).unwrap(), ops::log2);
+        assert_eq!(curve.sample_unchecked(3.5), ops::log2(3.5));
         assert!(curve.sample_unchecked(-1.0).is_nan());
         assert!(curve.sample(-1.0).is_none());
     }
@@ -642,10 +642,10 @@ mod tests {
 
     #[test]
     fn reparametrization() {
-        let curve = function_curve(interval(1.0, f32::INFINITY).unwrap(), f32::log2);
+        let curve = function_curve(interval(1.0, f32::INFINITY).unwrap(), ops::log2);
         let reparametrized_curve = curve
             .by_ref()
-            .reparametrize(interval(0.0, f32::INFINITY).unwrap(), f32::exp2);
+            .reparametrize(interval(0.0, f32::INFINITY).unwrap(), ops::exp2);
         assert_abs_diff_eq!(reparametrized_curve.sample_unchecked(3.5), 3.5);
         assert_abs_diff_eq!(reparametrized_curve.sample_unchecked(100.0), 100.0);
         assert_eq!(
@@ -664,8 +664,8 @@ mod tests {
     #[test]
     fn multiple_maps() {
         // Make sure these actually happen in the right order.
-        let curve = function_curve(interval(0.0, 1.0).unwrap(), f32::exp2);
-        let first_mapped = curve.map(f32::log2);
+        let curve = function_curve(interval(0.0, 1.0).unwrap(), ops::exp2);
+        let first_mapped = curve.map(ops::log2);
         let second_mapped = first_mapped.map(|x| x * -2.0);
         assert_abs_diff_eq!(second_mapped.sample_unchecked(0.0), 0.0);
         assert_abs_diff_eq!(second_mapped.sample_unchecked(0.5), -1.0);
@@ -675,8 +675,8 @@ mod tests {
     #[test]
     fn multiple_reparams() {
         // Make sure these happen in the right order too.
-        let curve = function_curve(interval(0.0, 1.0).unwrap(), f32::exp2);
-        let first_reparam = curve.reparametrize(interval(1.0, 2.0).unwrap(), f32::log2);
+        let curve = function_curve(interval(0.0, 1.0).unwrap(), ops::exp2);
+        let first_reparam = curve.reparametrize(interval(1.0, 2.0).unwrap(), ops::log2);
         let second_reparam = first_reparam.reparametrize(interval(0.0, 1.0).unwrap(), |t| t + 1.0);
         assert_abs_diff_eq!(second_reparam.sample_unchecked(0.0), 1.0);
         assert_abs_diff_eq!(second_reparam.sample_unchecked(0.5), 1.5);
diff --git a/crates/bevy_math/src/direction.rs b/crates/bevy_math/src/direction.rs
index 2ecb33ebf2..8f4cf7d9ea 100644
--- a/crates/bevy_math/src/direction.rs
+++ b/crates/bevy_math/src/direction.rs
@@ -860,6 +860,8 @@ impl approx::UlpsEq for Dir3A {
 
 #[cfg(test)]
 mod tests {
+    use crate::ops;
+
     use super::*;
     use approx::assert_relative_eq;
 
@@ -916,7 +918,7 @@ mod tests {
     #[test]
     fn dir2_renorm() {
         // Evil denormalized Rot2
-        let (sin, cos) = 1.0_f32.sin_cos();
+        let (sin, cos) = ops::sin_cos(1.0_f32);
         let rot2 = Rot2::from_sin_cos(sin * (1.0 + 1e-5), cos * (1.0 + 1e-5));
         let mut dir_a = Dir2::X;
         let mut dir_b = Dir2::X;
diff --git a/crates/bevy_math/src/lib.rs b/crates/bevy_math/src/lib.rs
index f34e452fce..8a61bd28e4 100644
--- a/crates/bevy_math/src/lib.rs
+++ b/crates/bevy_math/src/lib.rs
@@ -21,6 +21,7 @@ pub mod curve;
 mod direction;
 mod float_ord;
 mod isometry;
+mod ops;
 pub mod primitives;
 mod ray;
 mod rects;
diff --git a/crates/bevy_math/src/ops.rs b/crates/bevy_math/src/ops.rs
new file mode 100644
index 0000000000..798625d374
--- /dev/null
+++ b/crates/bevy_math/src/ops.rs
@@ -0,0 +1,290 @@
+//! This mod re-exports the correct versions of floating-point operations with
+//! unspecified precision in the standard library depending on whether the `libm`
+//! crate feature is enabled.
+//!
+//! All the functions here are named according to their versions in the standard
+//! library.
+
+#![allow(dead_code)]
+#![allow(clippy::disallowed_methods)]
+
+// Note: There are some Rust methods with unspecified precision without a `libm`
+// equivalent:
+// - `f32::powi` (integer powers)
+// - `f32::log` (logarithm with specified base)
+// - `f32::abs_sub` (actually unsure if `libm` has this, but don't use it regardless)
+//
+// Additionally, the following nightly API functions are not presently integrated
+// into this, but they would be candidates once standardized:
+// - `f32::gamma`
+// - `f32::ln_gamma`
+
+#[cfg(not(feature = "libm"))]
+mod std_ops {
+    #[inline(always)]
+    pub(crate) fn powf(x: f32, y: f32) -> f32 {
+        f32::powf(x, y)
+    }
+
+    #[inline(always)]
+    pub(crate) fn exp(x: f32) -> f32 {
+        f32::exp(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn exp2(x: f32) -> f32 {
+        f32::exp2(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn ln(x: f32) -> f32 {
+        f32::ln(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn log2(x: f32) -> f32 {
+        f32::log2(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn log10(x: f32) -> f32 {
+        f32::log10(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn cbrt(x: f32) -> f32 {
+        f32::cbrt(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn hypot(x: f32, y: f32) -> f32 {
+        f32::hypot(x, y)
+    }
+
+    #[inline(always)]
+    pub(crate) fn sin(x: f32) -> f32 {
+        f32::sin(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn cos(x: f32) -> f32 {
+        f32::cos(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn tan(x: f32) -> f32 {
+        f32::tan(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn asin(x: f32) -> f32 {
+        f32::asin(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn acos(x: f32) -> f32 {
+        f32::acos(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn atan(x: f32) -> f32 {
+        f32::atan(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn atan2(x: f32, y: f32) -> f32 {
+        f32::atan2(x, y)
+    }
+
+    #[inline(always)]
+    pub(crate) fn sin_cos(x: f32) -> (f32, f32) {
+        f32::sin_cos(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn exp_m1(x: f32) -> f32 {
+        f32::exp_m1(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn ln_1p(x: f32) -> f32 {
+        f32::ln_1p(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn sinh(x: f32) -> f32 {
+        f32::sinh(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn cosh(x: f32) -> f32 {
+        f32::cosh(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn tanh(x: f32) -> f32 {
+        f32::tanh(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn asinh(x: f32) -> f32 {
+        f32::asinh(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn acosh(x: f32) -> f32 {
+        f32::acosh(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn atanh(x: f32) -> f32 {
+        f32::atanh(x)
+    }
+}
+
+#[cfg(feature = "libm")]
+mod libm_ops {
+    #[inline(always)]
+    pub(crate) fn powf(x: f32, y: f32) -> f32 {
+        libm::powf(x, y)
+    }
+
+    #[inline(always)]
+    pub(crate) fn exp(x: f32) -> f32 {
+        libm::expf(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn exp2(x: f32) -> f32 {
+        libm::exp2f(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn ln(x: f32) -> f32 {
+        // This isn't documented in `libm` but this is actually the base e logarithm.
+        libm::logf(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn log2(x: f32) -> f32 {
+        libm::log2f(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn log10(x: f32) -> f32 {
+        libm::log10f(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn cbrt(x: f32) -> f32 {
+        libm::cbrtf(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn hypot(x: f32, y: f32) -> f32 {
+        libm::hypotf(x, y)
+    }
+
+    #[inline(always)]
+    pub(crate) fn sin(x: f32) -> f32 {
+        libm::sinf(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn cos(x: f32) -> f32 {
+        libm::cosf(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn tan(x: f32) -> f32 {
+        libm::tanf(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn asin(x: f32) -> f32 {
+        libm::asinf(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn acos(x: f32) -> f32 {
+        libm::acosf(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn atan(x: f32) -> f32 {
+        libm::atanf(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn atan2(x: f32, y: f32) -> f32 {
+        libm::atan2f(x, y)
+    }
+
+    #[inline(always)]
+    pub(crate) fn sin_cos(x: f32) -> (f32, f32) {
+        libm::sincosf(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn exp_m1(x: f32) -> f32 {
+        libm::expm1f(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn ln_1p(x: f32) -> f32 {
+        libm::log1pf(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn sinh(x: f32) -> f32 {
+        libm::sinhf(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn cosh(x: f32) -> f32 {
+        libm::coshf(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn tanh(x: f32) -> f32 {
+        libm::tanhf(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn asinh(x: f32) -> f32 {
+        libm::asinhf(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn acosh(x: f32) -> f32 {
+        libm::acoshf(x)
+    }
+
+    #[inline(always)]
+    pub(crate) fn atanh(x: f32) -> f32 {
+        libm::atanhf(x)
+    }
+}
+
+#[cfg(feature = "libm")]
+pub(crate) use libm_ops::*;
+#[cfg(not(feature = "libm"))]
+pub(crate) use std_ops::*;
+
+/// This extension trait covers shortfall in determinacy from the lack of a `libm` counterpart
+/// to `f32::powi`. Use this for the common small exponents.
+pub(crate) trait FloatPow {
+    fn squared(self) -> Self;
+    fn cubed(self) -> Self;
+}
+
+impl FloatPow for f32 {
+    #[inline]
+    fn squared(self) -> Self {
+        self * self
+    }
+    #[inline]
+    fn cubed(self) -> Self {
+        self * self * self
+    }
+}
diff --git a/crates/bevy_math/src/primitives/dim2.rs b/crates/bevy_math/src/primitives/dim2.rs
index e961603678..68e4cd3d7a 100644
--- a/crates/bevy_math/src/primitives/dim2.rs
+++ b/crates/bevy_math/src/primitives/dim2.rs
@@ -1,7 +1,10 @@
-use std::f32::consts::{FRAC_PI_2, FRAC_PI_3, PI};
+use std::f32::consts::{FRAC_1_SQRT_2, FRAC_PI_2, FRAC_PI_3, PI};
 
 use super::{Measured2d, Primitive2d, WindingOrder};
-use crate::{Dir2, Vec2};
+use crate::{
+    ops::{self, FloatPow},
+    Dir2, Vec2,
+};
 
 #[cfg(feature = "bevy_reflect")]
 use bevy_reflect::{std_traits::ReflectDefault, Reflect};
@@ -54,7 +57,7 @@ impl Circle {
     pub fn closest_point(&self, point: Vec2) -> Vec2 {
         let distance_squared = point.length_squared();
 
-        if distance_squared <= self.radius.powi(2) {
+        if distance_squared <= self.radius.squared() {
             // The point is inside the circle.
             point
         } else {
@@ -70,7 +73,7 @@ impl Measured2d for Circle {
     /// Get the area of the circle
     #[inline(always)]
     fn area(&self) -> f32 {
-        PI * self.radius.powi(2)
+        PI * self.radius.squared()
     }
 
     /// Get the perimeter or circumference of the circle
@@ -200,7 +203,7 @@ impl Arc2d {
     /// Get half the distance between the endpoints (half the length of the chord)
     #[inline(always)]
     pub fn half_chord_length(&self) -> f32 {
-        self.radius * f32::sin(self.half_angle)
+        self.radius * ops::sin(self.half_angle)
     }
 
     /// Get the distance between the endpoints (the length of the chord)
@@ -226,7 +229,7 @@ impl Arc2d {
     // used by Wolfram MathWorld, which is the distance rather than the segment.
     pub fn apothem(&self) -> f32 {
         let sign = if self.is_minor() { 1.0 } else { -1.0 };
-        sign * f32::sqrt(self.radius.powi(2) - self.half_chord_length().powi(2))
+        sign * f32::sqrt(self.radius.squared() - self.half_chord_length().squared())
     }
 
     /// Get the length of the sagitta of this arc, that is,
@@ -388,7 +391,7 @@ impl CircularSector {
     /// Returns the area of this sector
     #[inline(always)]
     pub fn area(&self) -> f32 {
-        self.arc.radius.powi(2) * self.arc.half_angle
+        self.arc.radius.squared() * self.arc.half_angle
     }
 }
 
@@ -527,7 +530,7 @@ impl CircularSegment {
     /// Returns the area of this segment
     #[inline(always)]
     pub fn area(&self) -> f32 {
-        0.5 * self.arc.radius.powi(2) * (self.arc.angle() - self.arc.angle().sin())
+        0.5 * self.arc.radius.squared() * (self.arc.angle() - ops::sin(self.arc.angle()))
     }
 }
 
@@ -882,11 +885,11 @@ impl Measured2d for Ellipse {
         // The algorithm used here is the Gauss-Kummer infinite series expansion of the elliptic integral expression for the perimeter of ellipses
         // For more information see https://www.wolframalpha.com/input/?i=gauss-kummer+series
         // We only use the terms up to `i == 20` for this approximation
-        let h = ((a - b) / (a + b)).powi(2);
+        let h = ((a - b) / (a + b)).squared();
 
         PI * (a + b)
             * (0..=20)
-                .map(|i| BINOMIAL_COEFFICIENTS[i] * h.powi(i as i32))
+                .map(|i| BINOMIAL_COEFFICIENTS[i] * ops::powf(h, i as f32))
                 .sum::<f32>()
     }
 }
@@ -953,8 +956,8 @@ impl Annulus {
     pub fn closest_point(&self, point: Vec2) -> Vec2 {
         let distance_squared = point.length_squared();
 
-        if self.inner_circle.radius.powi(2) <= distance_squared {
-            if distance_squared <= self.outer_circle.radius.powi(2) {
+        if self.inner_circle.radius.squared() <= distance_squared {
+            if distance_squared <= self.outer_circle.radius.squared() {
                 // The point is inside the annulus.
                 point
             } else {
@@ -976,7 +979,7 @@ impl Measured2d for Annulus {
     /// Get the area of the annulus
     #[inline(always)]
     fn area(&self) -> f32 {
-        PI * (self.outer_circle.radius.powi(2) - self.inner_circle.radius.powi(2))
+        PI * (self.outer_circle.radius.squared() - self.inner_circle.radius.squared())
     }
 
     /// Get the perimeter or circumference of the annulus,
@@ -1031,7 +1034,7 @@ impl Rhombus {
     #[inline(always)]
     pub fn from_side(side: f32) -> Self {
         Self {
-            half_diagonals: Vec2::splat(side.hypot(side) / 2.0),
+            half_diagonals: Vec2::splat(side * FRAC_1_SQRT_2),
         }
     }
 
@@ -1694,13 +1697,13 @@ impl RegularPolygon {
     #[inline(always)]
     #[doc(alias = "apothem")]
     pub fn inradius(&self) -> f32 {
-        self.circumradius() * (PI / self.sides as f32).cos()
+        self.circumradius() * ops::cos(PI / self.sides as f32)
     }
 
     /// Get the length of one side of the regular polygon
     #[inline(always)]
     pub fn side_length(&self) -> f32 {
-        2.0 * self.circumradius() * (PI / self.sides as f32).sin()
+        2.0 * self.circumradius() * ops::sin(PI / self.sides as f32)
     }
 
     /// Get the internal angle of the regular polygon in degrees.
@@ -1750,7 +1753,7 @@ impl RegularPolygon {
 
         (0..self.sides).map(move |i| {
             let theta = start_angle + i as f32 * step;
-            let (sin, cos) = theta.sin_cos();
+            let (sin, cos) = ops::sin_cos(theta);
             Vec2::new(cos, sin) * self.circumcircle.radius
         })
     }
@@ -1761,7 +1764,7 @@ impl Measured2d for RegularPolygon {
     #[inline(always)]
     fn area(&self) -> f32 {
         let angle: f32 = 2.0 * PI / (self.sides as f32);
-        (self.sides as f32) * self.circumradius().powi(2) * angle.sin() / 2.0
+        (self.sides as f32) * self.circumradius().squared() * ops::sin(angle) / 2.0
     }
 
     /// Get the perimeter of the regular polygon.
@@ -1821,7 +1824,7 @@ mod tests {
     // Reference values were computed by hand and/or with external tools
 
     use super::*;
-    use approx::assert_relative_eq;
+    use approx::{assert_abs_diff_eq, assert_relative_eq};
 
     #[test]
     fn rectangle_closest_point() {
@@ -1913,10 +1916,10 @@ mod tests {
         assert_eq!(rhombus.inradius(), 0.0, "incorrect inradius");
         assert_eq!(rhombus.circumradius(), 0.0, "incorrect circumradius");
         let rhombus = Rhombus::from_side(std::f32::consts::SQRT_2);
-        assert_eq!(rhombus, Rhombus::new(2.0, 2.0));
-        assert_eq!(
-            rhombus,
-            Rhombus::from_inradius(std::f32::consts::FRAC_1_SQRT_2)
+        assert_abs_diff_eq!(rhombus.half_diagonals, Vec2::new(1.0, 1.0));
+        assert_abs_diff_eq!(
+            rhombus.half_diagonals,
+            Rhombus::from_inradius(std::f32::consts::FRAC_1_SQRT_2).half_diagonals
         );
     }
 
diff --git a/crates/bevy_math/src/primitives/dim3.rs b/crates/bevy_math/src/primitives/dim3.rs
index 9c21b9ee4a..ed692d7ead 100644
--- a/crates/bevy_math/src/primitives/dim3.rs
+++ b/crates/bevy_math/src/primitives/dim3.rs
@@ -1,7 +1,7 @@
 use std::f32::consts::{FRAC_PI_3, PI};
 
 use super::{Circle, Measured2d, Measured3d, Primitive2d, Primitive3d};
-use crate::{Dir3, InvalidDirectionError, Mat3, Vec2, Vec3};
+use crate::{ops, ops::FloatPow, Dir3, InvalidDirectionError, Mat3, Vec2, Vec3};
 
 #[cfg(feature = "bevy_reflect")]
 use bevy_reflect::{std_traits::ReflectDefault, Reflect};
@@ -54,7 +54,7 @@ impl Sphere {
     pub fn closest_point(&self, point: Vec3) -> Vec3 {
         let distance_squared = point.length_squared();
 
-        if distance_squared <= self.radius.powi(2) {
+        if distance_squared <= self.radius.squared() {
             // The point is inside the sphere.
             point
         } else {
@@ -70,13 +70,13 @@ impl Measured3d for Sphere {
     /// Get the surface area of the sphere
     #[inline(always)]
     fn area(&self) -> f32 {
-        4.0 * PI * self.radius.powi(2)
+        4.0 * PI * self.radius.squared()
     }
 
     /// Get the volume of the sphere
     #[inline(always)]
     fn volume(&self) -> f32 {
-        4.0 * FRAC_PI_3 * self.radius.powi(3)
+        4.0 * FRAC_PI_3 * self.radius.cubed()
     }
 }
 
@@ -505,7 +505,7 @@ impl Cylinder {
     /// Get the surface area of one base of the cylinder
     #[inline(always)]
     pub fn base_area(&self) -> f32 {
-        PI * self.radius.powi(2)
+        PI * self.radius.squared()
     }
 }
 
@@ -642,7 +642,7 @@ impl Cone {
     #[inline(always)]
     #[doc(alias = "side_length")]
     pub fn slant_height(&self) -> f32 {
-        self.radius.hypot(self.height)
+        ops::hypot(self.radius, self.height)
     }
 
     /// Get the surface area of the side of the cone,
@@ -656,7 +656,7 @@ impl Cone {
     /// Get the surface area of the base of the cone
     #[inline(always)]
     pub fn base_area(&self) -> f32 {
-        PI * self.radius.powi(2)
+        PI * self.radius.squared()
     }
 }
 
@@ -828,14 +828,14 @@ impl Measured3d for Torus {
     /// the expected result when the torus has a ring and isn't self-intersecting
     #[inline(always)]
     fn area(&self) -> f32 {
-        4.0 * PI.powi(2) * self.major_radius * self.minor_radius
+        4.0 * PI.squared() * self.major_radius * self.minor_radius
     }
 
     /// Get the volume of the torus. Note that this only produces
     /// the expected result when the torus has a ring and isn't self-intersecting
     #[inline(always)]
     fn volume(&self) -> f32 {
-        2.0 * PI.powi(2) * self.major_radius * self.minor_radius.powi(2)
+        2.0 * PI.squared() * self.major_radius * self.minor_radius.squared()
     }
 }
 
diff --git a/crates/bevy_math/src/rotation2d.rs b/crates/bevy_math/src/rotation2d.rs
index 5f75ab804c..e9fd1f7f13 100644
--- a/crates/bevy_math/src/rotation2d.rs
+++ b/crates/bevy_math/src/rotation2d.rs
@@ -1,6 +1,9 @@
 use glam::FloatExt;
 
-use crate::prelude::{Mat2, Vec2};
+use crate::{
+    ops,
+    prelude::{Mat2, Vec2},
+};
 
 #[cfg(feature = "bevy_reflect")]
 use bevy_reflect::{std_traits::ReflectDefault, Reflect};
@@ -99,14 +102,7 @@ impl Rot2 {
     /// Creates a [`Rot2`] from a counterclockwise angle in radians.
     #[inline]
     pub fn radians(radians: f32) -> Self {
-        #[cfg(feature = "libm")]
-        let (sin, cos) = (
-            libm::sin(radians as f64) as f32,
-            libm::cos(radians as f64) as f32,
-        );
-        #[cfg(not(feature = "libm"))]
-        let (sin, cos) = radians.sin_cos();
-
+        let (sin, cos) = ops::sin_cos(radians);
         Self::from_sin_cos(sin, cos)
     }
 
@@ -136,14 +132,7 @@ impl Rot2 {
     /// Returns the rotation in radians in the `(-pi, pi]` range.
     #[inline]
     pub fn as_radians(self) -> f32 {
-        #[cfg(feature = "libm")]
-        {
-            libm::atan2(self.sin as f64, self.cos as f64) as f32
-        }
-        #[cfg(not(feature = "libm"))]
-        {
-            f32::atan2(self.sin, self.cos)
-        }
+        ops::atan2(self.sin, self.cos)
     }
 
     /// Returns the rotation in degrees in the `(-180, 180]` range.
diff --git a/crates/bevy_math/src/sampling/shape_sampling.rs b/crates/bevy_math/src/sampling/shape_sampling.rs
index 411c4b99f8..10f9e17f60 100644
--- a/crates/bevy_math/src/sampling/shape_sampling.rs
+++ b/crates/bevy_math/src/sampling/shape_sampling.rs
@@ -40,7 +40,7 @@
 
 use std::f32::consts::{PI, TAU};
 
-use crate::{primitives::*, NormedVectorSpace, Vec2, Vec3};
+use crate::{ops, primitives::*, NormedVectorSpace, Vec2, Vec3};
 use rand::{
     distributions::{Distribution, WeightedIndex},
     Rng,
@@ -155,12 +155,14 @@ impl ShapeSample for Circle {
         let theta = rng.gen_range(0.0..TAU);
         let r_squared = rng.gen_range(0.0..=(self.radius * self.radius));
         let r = r_squared.sqrt();
-        Vec2::new(r * theta.cos(), r * theta.sin())
+        let (sin, cos) = ops::sin_cos(theta);
+        Vec2::new(r * cos, r * sin)
     }
 
     fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
         let theta = rng.gen_range(0.0..TAU);
-        Vec2::new(self.radius * theta.cos(), self.radius * theta.sin())
+        let (sin, cos) = ops::sin_cos(theta);
+        Vec2::new(self.radius * cos, self.radius * sin)
     }
 }
 
@@ -168,7 +170,7 @@ impl ShapeSample for Circle {
 #[inline]
 fn sample_unit_sphere_boundary<R: Rng + ?Sized>(rng: &mut R) -> Vec3 {
     let z = rng.gen_range(-1f32..=1f32);
-    let (a_sin, a_cos) = rng.gen_range(-PI..=PI).sin_cos();
+    let (a_sin, a_cos) = ops::sin_cos(rng.gen_range(-PI..=PI));
     let c = (1f32 - z * z).sqrt();
     let x = a_sin * c;
     let y = a_cos * c;
@@ -181,7 +183,7 @@ impl ShapeSample for Sphere {
 
     fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
         let r_cubed = rng.gen_range(0.0..=(self.radius * self.radius * self.radius));
-        let r = r_cubed.cbrt();
+        let r = ops::cbrt(r_cubed);
 
         r * sample_unit_sphere_boundary(rng)
     }
@@ -202,8 +204,9 @@ impl ShapeSample for Annulus {
         let r_squared = rng.gen_range((inner_radius * inner_radius)..(outer_radius * outer_radius));
         let r = r_squared.sqrt();
         let theta = rng.gen_range(0.0..TAU);
+        let (sin, cos) = ops::sin_cos(theta);
 
-        Vec2::new(r * theta.cos(), r * theta.sin())
+        Vec2::new(r * cos, r * sin)
     }
 
     fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
@@ -624,7 +627,7 @@ mod tests {
         for _ in 0..5000 {
             let point = circle.sample_boundary(&mut rng);
 
-            let angle = f32::atan(point.y / point.x) + PI / 2.0;
+            let angle = ops::atan(point.y / point.x) + PI / 2.0;
             let wedge = (angle * 8.0 / PI).floor() as usize;
             wedge_hits[wedge] += 1;
         }