Repository: bitshifter/glam-rs
Pull request #389 (2023-04-02 11:02:54)

Make enabling libm override std math functions


Patch:
diff --git a/Cargo.toml b/Cargo.toml
index a5d784a..82febfb 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -28,9 +28,6 @@ glam-assert = []
 # this is primarily for testing the fallback implementation
 scalar-math = []
 
-# libm is required when building no_std
-libm = ["num-traits/libm"]
-
 # align types to match CUDA requirements
 cuda = []
 
@@ -48,11 +45,11 @@ core-simd = []
 approx = { version = "0.5", optional = true, default-features = false }
 bytemuck = { version = "1.5", optional = true, default-features = false }
 mint = { version = "0.5.8", optional = true, default-features = false }
-num-traits = { version = "0.2.14", optional = true, default-features = false }
 rand = { version = "0.8", optional = true, default-features = false }
 serde = { version = "1.0", optional = true, default-features = false }
 rkyv = { version = "0.7", optional = true }
 bytecheck = { version = "0.6", optional = true, default-features = false}
+libm = { version = "0.2", optional = true, default-features = false}
 
 [dev-dependencies]
 # rand_xoshiro is required for tests if rand is enabled
diff --git a/README.md b/README.md
index dde7a4c..8baff6d 100644
--- a/README.md
+++ b/README.md
@@ -95,25 +95,23 @@ glam = { version = "0.23", default-features = false }
 
 * [`approx`] - traits and macros for approximate float comparisons
 * [`bytemuck`] - for casting into slices of bytes
-* [`libm`] - required to compile with `no_std`
+* [`libm`] - uses `libm` math functions instead of `std`, required to compile
+  with `no_std`
 * [`mint`] - for interoperating with other 3D math libraries
-* [`num-traits`] - required to compile `no_std`, will be included when enabling
-  the `libm` feature
 * [`rand`] - implementations of `Distribution` trait for all `glam` types.
 * [`serde`] - implementations of `Serialize` and `Deserialize` for all `glam`
   types. Note that serialization should work between builds of `glam` with and
   without SIMD enabled
-* [`rkyv`] - implementations of `Archive`, `Serialize` and `Deserialize` for all
-  `glam` types. Note that serialization is not interoperable with and without the
-  `scalar-math` feature. It should work between all other builds of `glam`.
-  Endian conversion is currently not supported
+* [`rkyv`] - implementations of `Archive`, `Serialize` and `Deserialize` for
+  all `glam` types. Note that serialization is not interoperable with and
+  without the `scalar-math` feature. It should work between all other builds of
+  `glam`.  Endian conversion is currently not supported
 * [`bytecheck`] - to perform archive validation when using the `rkyv` feature
 
 [`approx`]: https://docs.rs/approx
 [`bytemuck`]: https://docs.rs/bytemuck
 [`libm`]: https://github.com/rust-lang/libm
 [`mint`]: https://github.com/kvark/mint
-[`num-traits`]: https://github.com/rust-num/num-traits
 [`rand`]: https://github.com/rust-random/rand
 [`serde`]: https://serde.rs
 [`rkyv`]: https://github.com/rkyv/rkyv
diff --git a/codegen/templates/affine.rs.tera b/codegen/templates/affine.rs.tera
index a9c3b42..edb5405 100644
--- a/codegen/templates/affine.rs.tera
+++ b/codegen/templates/affine.rs.tera
@@ -454,15 +454,12 @@ impl {{ self_t }} {
     /// vector contains any zero elements when `glam_assert` is enabled.
     #[inline]
     pub fn to_scale_rotation_translation(&self) -> ({{ vec3_t }}, {{ quat_t }}, {{ vec3_t }}) {
-        #[cfg(feature = "libm")]
-        #[allow(unused_imports)]
-        use num_traits::Float;
-
+        use crate::{{ scalar_t }}::math;
         let det = self.matrix3.determinant();
         glam_assert!(det != 0.0);
 
         let scale = {{ vec3_t }}::new(
-            self.matrix3.x_axis.length() * det.signum(),
+            self.matrix3.x_axis.length() * math::signum(det),
             self.matrix3.y_axis.length(),
             self.matrix3.z_axis.length(),
         );
diff --git a/codegen/templates/mat.rs.tera b/codegen/templates/mat.rs.tera
index a187d90..33ab8d8 100644
--- a/codegen/templates/mat.rs.tera
+++ b/codegen/templates/mat.rs.tera
@@ -91,6 +91,7 @@ use crate::{
         coresimd::*,
     {% endif %}
 {% endif %}
+    {{ scalar_t }}::math,
     swizzles::*,
 };
 #[cfg(not(target_arch = "spirv"))]
@@ -109,10 +110,6 @@ use core::arch::wasm32::*;
 use core::simd::{Which::*, *};
 {% endif %}
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 {% if self_t == "Mat2" and is_sse2 %}
 union UnionCast {
     a: [f32; 4],
@@ -398,14 +395,14 @@ impl {{ self_t }} {
     /// `angle` (in radians).
     #[inline]
     pub fn from_scale_angle(scale: {{ col_t }}, angle: {{ scalar_t }}) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::new(cos * scale.x, sin * scale.x, -sin * scale.y, cos * scale.y)
     }
 
     /// Creates a {{ nxn }} matrix containing a rotation of `angle` (in radians).
     #[inline]
     pub fn from_angle(angle: {{ scalar_t }}) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::new(cos, sin, -sin, cos)
     }
 
@@ -482,7 +479,7 @@ impl {{ self_t }} {
         {# TODO: make common with dim == 4 #}
         glam_assert!(axis.is_normalized());
 
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         let (xsin, ysin, zsin) = axis.mul(sin).into();
         let (x, y, z) = axis.into();
         let (x2, y2, z2) = axis.mul(axis).into();
@@ -508,7 +505,7 @@ impl {{ self_t }} {
     /// Creates a 3D rotation matrix from `angle` (in radians) around the x axis.
     #[inline]
     pub fn from_rotation_x(angle: {{ scalar_t }}) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             {{ col_t }}::X,
             {{ col_t }}::new(0.0, cosa, sina),
@@ -519,7 +516,7 @@ impl {{ self_t }} {
     /// Creates a 3D rotation matrix from `angle` (in radians) around the y axis.
     #[inline]
     pub fn from_rotation_y(angle: {{ scalar_t }}) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             {{ col_t }}::new(cosa, 0.0, -sina),
             {{ col_t }}::Y,
@@ -530,7 +527,7 @@ impl {{ self_t }} {
     /// Creates a 3D rotation matrix from `angle` (in radians) around the z axis.
     #[inline]
     pub fn from_rotation_z(angle: {{ scalar_t }}) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             {{ col_t }}::new(cosa, sina, 0.0),
             {{ col_t }}::new(-sina, cosa, 0.0),
@@ -557,7 +554,7 @@ impl {{ self_t }} {
     /// [`Self::transform_point2()`] and [`Self::transform_vector2()`].
     #[inline]
     pub fn from_angle(angle: {{ scalar_t }}) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::from_cols(
             {{ col_t }}::new(cos, sin, 0.0),
             {{ col_t }}::new(-sin, cos, 0.0),
@@ -572,7 +569,7 @@ impl {{ self_t }} {
     /// [`Self::transform_point2()`] and [`Self::transform_vector2()`].
     #[inline]
     pub fn from_scale_angle_translation(scale: {{ vec2_t }}, angle: {{ scalar_t }}, translation: {{ vec2_t }}) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::from_cols(
             {{ col_t }}::new(cos * scale.x, sin * scale.x, 0.0),
             {{ col_t }}::new(-sin * scale.y, cos * scale.y, 0.0),
@@ -685,7 +682,7 @@ impl {{ self_t }} {
         glam_assert!(det != 0.0);
 
         let scale = {{ vec3_t }}::new(
-            self.x_axis.length() * det.signum(),
+            self.x_axis.length() * math::signum(det),
             self.y_axis.length(),
             self.z_axis.length(),
         );
@@ -779,7 +776,7 @@ impl {{ self_t }} {
         {# TODO: make common with dim == 3 #}
         glam_assert!(axis.is_normalized());
 
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         let axis_sin = axis.mul(sin);
         let axis_sq = axis.mul(axis);
         let omc = 1.0 - cos;
@@ -827,7 +824,7 @@ impl {{ self_t }} {
     /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
     #[inline]
     pub fn from_rotation_x(angle: {{ scalar_t }}) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             {{ col_t }}::X,
             {{ col_t }}::new(0.0, cosa, sina, 0.0),
@@ -843,7 +840,7 @@ impl {{ self_t }} {
     /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
     #[inline]
     pub fn from_rotation_y(angle: {{ scalar_t }}) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             {{ col_t }}::new(cosa, 0.0, -sina, 0.0),
             {{ col_t }}::Y,
@@ -859,7 +856,7 @@ impl {{ self_t }} {
     /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
     #[inline]
     pub fn from_rotation_z(angle: {{ scalar_t }}) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             {{ col_t }}::new(cosa, sina, 0.0, 0.0),
             {{ col_t }}::new(-sina, cosa, 0.0, 0.0),
@@ -1509,7 +1506,7 @@ impl {{ self_t }} {
         z_far: {{ scalar_t }},
     ) -> Self {
         let inv_length = 1.0 / (z_near - z_far);
-        let f = 1.0 / (0.5 * fov_y_radians).tan();
+        let f = 1.0 / math::tan(0.5 * fov_y_radians);
         let a = f / aspect_ratio;
         let b = (z_near + z_far) * inv_length;
         let c = (2.0 * z_near * z_far) * inv_length;
@@ -1530,7 +1527,7 @@ impl {{ self_t }} {
     #[inline]
     pub fn perspective_lh(fov_y_radians: {{ scalar_t }}, aspect_ratio: {{ scalar_t }}, z_near: {{ scalar_t }}, z_far: {{ scalar_t }}) -> Self {
         glam_assert!(z_near > 0.0 && z_far > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         let r = z_far / (z_far - z_near);
@@ -1551,7 +1548,7 @@ impl {{ self_t }} {
     #[inline]
     pub fn perspective_rh(fov_y_radians: {{ scalar_t }}, aspect_ratio: {{ scalar_t }}, z_near: {{ scalar_t }}, z_far: {{ scalar_t }}) -> Self {
         glam_assert!(z_near > 0.0 && z_far > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         let r = z_far / (z_near - z_far);
@@ -1571,7 +1568,7 @@ impl {{ self_t }} {
     #[inline]
     pub fn perspective_infinite_lh(fov_y_radians: {{ scalar_t }}, aspect_ratio: {{ scalar_t }}, z_near: {{ scalar_t }}) -> Self {
         glam_assert!(z_near > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         Self::from_cols(
@@ -1594,7 +1591,7 @@ impl {{ self_t }} {
         z_near: {{ scalar_t }},
     ) -> Self {
         glam_assert!(z_near > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         Self::from_cols(
@@ -1610,7 +1607,7 @@ impl {{ self_t }} {
     #[inline]
     pub fn perspective_infinite_rh(fov_y_radians: {{ scalar_t }}, aspect_ratio: {{ scalar_t }}, z_near: {{ scalar_t }}) -> Self {
         glam_assert!(z_near > 0.0);
-        let f = 1.0 / (0.5 * fov_y_radians).tan();
+        let f = 1.0 / math::tan(0.5 * fov_y_radians);
         Self::from_cols(
             {{ col_t }}::new(f / aspect_ratio, 0.0, 0.0, 0.0),
             {{ col_t }}::new(0.0, f, 0.0, 0.0),
@@ -1628,7 +1625,7 @@ impl {{ self_t }} {
         z_near: {{ scalar_t }},
     ) -> Self {
         glam_assert!(z_near > 0.0);
-        let f = 1.0 / (0.5 * fov_y_radians).tan();
+        let f = 1.0 / math::tan(0.5 * fov_y_radians);
         Self::from_cols(
             {{ col_t }}::new(f / aspect_ratio, 0.0, 0.0, 0.0),
             {{ col_t }}::new(0.0, f, 0.0, 0.0),
diff --git a/codegen/templates/quat.rs.tera b/codegen/templates/quat.rs.tera
index 1686c04..d5071e5 100644
--- a/codegen/templates/quat.rs.tera
+++ b/codegen/templates/quat.rs.tera
@@ -30,7 +30,7 @@
 {% endif %}
 
 use crate::{
-    FloatEx,
+    {{ scalar_t }}::math,
     euler::{EulerFromQuaternion, EulerRot, EulerToQuaternion},
     {% if scalar_t == "f32" %}
         DQuat, Mat3, Mat3A, Mat4, Vec2, Vec3, Vec3A, Vec4,
@@ -46,10 +46,6 @@ use crate::{
     {% endif %}
 };
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 {% if is_sse2 %}
 #[cfg(target_arch = "x86")]
 use core::arch::x86::*;
@@ -227,7 +223,7 @@ impl {{ self_t }} {
     #[inline]
     pub fn from_axis_angle(axis: {{ vec3_t }}, angle: {{ scalar_t }}) -> Self {
         glam_assert!(axis.is_normalized());
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         let v = axis * s;
         Self::from_xyzw(v.x, v.y, v.z, c)
     }
@@ -248,21 +244,21 @@ impl {{ self_t }} {
     /// Creates a quaternion from the `angle` (in radians) around the x axis.
     #[inline]
     pub fn from_rotation_x(angle: {{ scalar_t }}) -> Self {
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         Self::from_xyzw(s, 0.0, 0.0, c)
     }
 
     /// Creates a quaternion from the `angle` (in radians) around the y axis.
     #[inline]
     pub fn from_rotation_y(angle: {{ scalar_t }}) -> Self {
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         Self::from_xyzw(0.0, s, 0.0, c)
     }
 
     /// Creates a quaternion from the `angle` (in radians) around the z axis.
     #[inline]
     pub fn from_rotation_z(angle: {{ scalar_t }}) -> Self {
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         Self::from_xyzw(0.0, 0.0, s, c)
     }
 
@@ -286,7 +282,7 @@ impl {{ self_t }} {
             if dif10 <= 0.0 {
                 // x^2 >= y^2
                 let four_xsq = omm22 - dif10;
-                let inv4x = 0.5 / four_xsq.sqrt();
+                let inv4x = 0.5 / math::sqrt(four_xsq);
                 Self::from_xyzw(
                     four_xsq * inv4x,
                     (m01 + m10) * inv4x,
@@ -296,7 +292,7 @@ impl {{ self_t }} {
             } else {
                 // y^2 >= x^2
                 let four_ysq = omm22 + dif10;
-                let inv4y = 0.5 / four_ysq.sqrt();
+                let inv4y = 0.5 / math::sqrt(four_ysq);
                 Self::from_xyzw(
                     (m01 + m10) * inv4y,
                     four_ysq * inv4y,
@@ -311,7 +307,7 @@ impl {{ self_t }} {
             if sum10 <= 0.0 {
                 // z^2 >= w^2
                 let four_zsq = opm22 - sum10;
-                let inv4z = 0.5 / four_zsq.sqrt();
+                let inv4z = 0.5 / math::sqrt(four_zsq);
                 Self::from_xyzw(
                     (m02 + m20) * inv4z,
                     (m12 + m21) * inv4z,
@@ -321,7 +317,7 @@ impl {{ self_t }} {
             } else {
                 // w^2 >= z^2
                 let four_wsq = opm22 + sum10;
-                let inv4w = 0.5 / four_wsq.sqrt();
+                let inv4w = 0.5 / math::sqrt(four_wsq);
                 Self::from_xyzw(
                     (m12 - m21) * inv4w,
                     (m20 - m02) * inv4w,
@@ -451,7 +447,7 @@ impl {{ self_t }} {
             let z = from.x * to.y - to.x * from.y;
             let w = 1.0 + dot;
             // calculate length with x=0 and y=0 to normalize
-            let len_rcp = 1.0 / (z * z + w * w).sqrt();
+            let len_rcp = 1.0 / math::sqrt(z * z + w * w);
             Self::from_xyzw(0.0, 0.0, z * len_rcp, w * len_rcp)
         }
     }
@@ -463,7 +459,7 @@ impl {{ self_t }} {
         let v = {{ vec3_t }}::new(self.x, self.y, self.z);
         let length = v.length();
         if length >= EPSILON {
-            let angle = 2.0 * {{ scalar_t }}::atan2(length, self.w);
+            let angle = 2.0 * math::atan2(length, self.w);
             let axis = v / length;
             (axis, angle)
         } else {
@@ -613,7 +609,7 @@ impl {{ self_t }} {
         // If the quat.w is close to -1.0, the angle will be near 2*PI which is close to
         // a negative 0 rotation. By forcing quat.w to be positive, we'll end up with
         // the shortest path.
-        let positive_w_angle = self.w.abs().acos_approx() * 2.0;
+        let positive_w_angle = math::acos_approx(math::abs(self.w)) * 2.0;
         positive_w_angle < threshold_angle
     }
 
@@ -628,7 +624,7 @@ impl {{ self_t }} {
     #[inline]
     pub fn angle_between(self, rhs: Self) -> {{ scalar_t }} {
         glam_assert!(self.is_normalized() && rhs.is_normalized());
-        self.dot(rhs).abs().acos_approx() * 2.0
+        math::acos_approx(math::abs(self.dot(rhs))) * 2.0
     }
 
     /// Returns true if the absolute difference of all elements between `self` and `rhs`
@@ -740,15 +736,15 @@ impl {{ self_t }} {
             // assumes lerp returns a normalized quaternion
             self.lerp(end, s)
         } else {
-            let theta = dot.acos_approx();
+            let theta = math::acos_approx(dot);
             {% if is_scalar %}
-                let scale1 = (theta * (1.0 - s)).sin();
-                let scale2 = (theta * s).sin();
-                let theta_sin = theta.sin();
+                let scale1 = math::sin(theta * (1.0 - s));
+                let scale2 = math::sin(theta * s);
+                let theta_sin = math::sin(theta);
 
                 self.mul(scale1)
                     .add(end.mul(scale2))
-                    .mul(theta_sin.recip())
+                    .mul(1.0 / theta_sin)
             {% elif is_sse2 %}
                 let x = 1.0 - s;
                 let y = s;
@@ -775,9 +771,9 @@ impl {{ self_t }} {
                 // let w = 0.0;
                 // let tmp = f32x4_mul(f32x4_splat(theta), f32x4(x, y, z, w));
                 // let tmp = v128_sin(tmp);
-                let x = (theta * (1.0 - s)).sin();
-                let y = (theta * s).sin();
-                let z = theta.sin();
+                let x = math::sin(theta * (1.0 - s));
+                let y = math::sin(theta * s);
+                let z = math::sin(theta);
                 let w = 0.0;
                 let tmp = f32x4(x, y, z, w);
 
@@ -790,9 +786,9 @@ impl {{ self_t }} {
                     theta_sin,
                 ))
             {% elif is_coresimd %}
-                let x = (theta * (1.0 - s)).sin();
-                let y = (theta * s).sin();
-                let z = theta.sin();
+                let x = math::sin(theta * (1.0 - s));
+                let y = math::sin(theta * s);
+                let z = math::sin(theta);
                 let w = 0.0;
                 let tmp = f32x4::from_array([x, y, z, w]);
 
diff --git a/codegen/templates/vec.rs.tera b/codegen/templates/vec.rs.tera
index cfb0f58..bc05278 100644
--- a/codegen/templates/vec.rs.tera
+++ b/codegen/templates/vec.rs.tera
@@ -104,6 +104,9 @@ use crate::{
     {% elif is_coresimd %}
         coresimd::*,
     {% endif %}
+    {% if is_float %}
+        {{ scalar_t }}::math,
+    {% endif %}
 };
 
 #[cfg(not(target_arch = "spirv"))]
@@ -123,12 +126,6 @@ use core::simd::*;
 use std::simd::StdFloat;
 {% endif %}
 
-{% if is_float %}
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-{% endif %}
-
 {% if is_sse2 or is_wasm32 %}
 union UnionCast {
     a: [f32; 4],
@@ -804,7 +801,11 @@ impl {{ self_t }} {
         {% if is_scalar %}
             Self {
                 {% for c in components %}
-                    {{ c }}: self.{{ c }}.abs(),
+                    {%- if is_float %}
+                        {{ c }}: math::abs(self.{{ c }}),
+                    {%- else %}
+                        {{ c }}: self.{{ c }}.abs(),
+                    {%- endif %}
                 {%- endfor %}
             }
         {% elif is_sse2 %}
@@ -832,7 +833,11 @@ impl {{ self_t }} {
         {% if is_scalar %}
             Self {
                 {% for c in components %}
-                    {{ c }}: self.{{ c }}.signum(),
+                    {%- if is_float %}
+                        {{ c }}: math::signum(self.{{ c }}),
+                    {%- else %}
+                        {{ c }}: self.{{ c }}.signum(),
+                    {%- endif %}
                 {%- endfor %}
             }
         {% elif is_coresimd %}
@@ -858,7 +863,7 @@ impl {{ self_t }} {
         {% if is_scalar and is_float %}
             Self {
                 {% for c in components %}
-                    {{ c }}: self.{{ c }}.copysign(rhs.{{ c }}),
+                    {{ c }}: math::copysign(self.{{ c }}, rhs.{{ c }}),
                 {%- endfor %}
             }
         {% elif is_scalar %}
@@ -977,7 +982,7 @@ impl {{ self_t }} {
     #[inline]
     pub fn length(self) -> {{ scalar_t }} {
         {% if is_scalar %}
-            self.dot(self).sqrt()
+            math::sqrt(self.dot(self))
         {% elif is_sse2 %}
             unsafe {
                 let dot = dot{{ dim }}_in_x(self.0, self.0);
@@ -1122,7 +1127,7 @@ impl {{ self_t }} {
     #[inline]
     pub fn is_normalized(self) -> bool {
         // TODO: do something with epsilon
-        (self.length_squared() - 1.0).abs() <= 1e-4
+        math::abs(self.length_squared() - 1.0) <= 1e-4
     }
 
     /// Returns the vector projection of `self` onto `rhs`.
@@ -1193,7 +1198,7 @@ impl {{ self_t }} {
         {% if is_scalar %}
             Self {
                 {% for c in components %}
-                    {{ c }}: self.{{ c }}.round(),
+                    {{ c }}: math::round(self.{{ c }}),
                 {%- endfor %}
             }
         {% elif is_sse2 %}
@@ -1212,7 +1217,7 @@ impl {{ self_t }} {
         {% if is_scalar %}
             Self {
                 {% for c in components %}
-                    {{ c }}: self.{{ c }}.floor(),
+                    {{ c }}: math::floor(self.{{ c }}),
                 {%- endfor %}
             }
         {% elif is_sse2 %}
@@ -1231,7 +1236,7 @@ impl {{ self_t }} {
         {% if is_scalar %}
             Self {
                 {% for c in components %}
-                    {{ c }}: self.{{ c }}.ceil(),
+                    {{ c }}: math::ceil(self.{{ c }}),
                 {%- endfor %}
             }
         {% elif is_sse2 %}
@@ -1258,7 +1263,7 @@ impl {{ self_t }} {
     pub fn exp(self) -> Self {
         Self::new(
             {% for c in components %}
-                self.{{ c }}.exp(),
+                math::exp(self.{{ c }}),
             {%- endfor %}
         )
     }
@@ -1268,7 +1273,7 @@ impl {{ self_t }} {
     pub fn powf(self, n: {{ scalar_t }}) -> Self {
         Self::new(
             {% for c in components %}
-                self.{{ c }}.powf(n),
+                math::powf(self.{{ c }}, n),
             {%- endfor %}
         )
     }
@@ -1279,7 +1284,7 @@ impl {{ self_t }} {
         {% if is_scalar %}
             Self {
                 {% for c in components %}
-                    {{ c }}: self.{{ c }}.recip(),
+                    {{ c }}: 1.0 / self.{{ c }},
                 {%- endfor %}
             }
         {% elif is_sse2 %}
@@ -1326,9 +1331,9 @@ impl {{ self_t }} {
         glam_assert!(min <= max);
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -1338,7 +1343,7 @@ impl {{ self_t }} {
     pub fn clamp_length_max(self, max: {{ scalar_t }}) -> Self {
         let length_sq = self.length_squared();
         if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -1348,7 +1353,7 @@ impl {{ self_t }} {
     pub fn clamp_length_min(self, min: {{ scalar_t }}) -> Self {
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -1373,7 +1378,7 @@ impl {{ self_t }} {
         {% else %}
             Self::new(
                 {% for c in components %}
-                    self.{{ c }}.mul_add(a.{{ c }}, b.{{ c }}),
+                    math::mul_add(self.{{ c }}, a.{{ c }}, b.{{ c }}),
                 {%- endfor %}
             )
         {% endif %}
@@ -1385,7 +1390,7 @@ impl {{ self_t }} {
     /// create the vector [-1, 0] and rotate `Vec2::Y` around it returning `-Vec2::Y`.
     #[inline]
     pub fn from_angle(angle: {{ scalar_t }}) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self {
             x: cos,
             y: sin,
@@ -1397,11 +1402,10 @@ impl {{ self_t }} {
     /// The input vectors do not need to be unit length however they must be non-zero.
     #[inline]
     pub fn angle_between(self, rhs: Self) -> {{ scalar_t }} {
-        use crate::FloatEx;
-        let angle = (self.dot(rhs) / (self.length_squared() * rhs.length_squared()).sqrt())
-            .acos_approx();
+        let angle = math::acos_approx(
+            self.dot(rhs) / math::sqrt(self.length_squared() * rhs.length_squared()));
 
-        angle * self.perp_dot(rhs).signum()
+        angle * math::signum(self.perp_dot(rhs))
     }
 {% elif dim == 3 %}
     /// Returns the angle (in radians) between two vectors.
@@ -1409,10 +1413,9 @@ impl {{ self_t }} {
     /// The input vectors do not need to be unit length however they must be non-zero.
     #[inline]
     pub fn angle_between(self, rhs: Self) -> {{ scalar_t }} {
-        use crate::FloatEx;
-        self.dot(rhs)
-            .div(self.length_squared().mul(rhs.length_squared()).sqrt())
-            .acos_approx()
+        math::acos_approx(
+            self.dot(rhs).div(
+                math::sqrt(self.length_squared().mul(rhs.length_squared()))))
     }
 
     /// Returns some vector that is orthogonal to the given one.
@@ -1424,7 +1427,7 @@ impl {{ self_t }} {
     #[inline]
     pub fn any_orthogonal_vector(&self) -> Self {
         // This can probably be optimized
-        if self.x.abs() > self.y.abs() {
+        if math::abs(self.x) > math::abs(self.y) {
             Self::new(-self.z, 0.0, self.x) // self.cross(Self::Y)
         } else {
             Self::new(0.0, self.z, -self.y) // self.cross(Self::X)
@@ -1441,10 +1444,7 @@ impl {{ self_t }} {
     pub fn any_orthonormal_vector(&self) -> Self {
         glam_assert!(self.is_normalized());
         // From https://graphics.pixar.com/library/OrthonormalB/paper.pdf
-        #[cfg(feature = "std")]
-        let sign = (1.0_{{ scalar_t }}).copysign(self.z);
-        #[cfg(not(feature = "std"))]
-        let sign = self.z.signum();
+        let sign = math::signum(self.z);
         let a = -1.0 / (sign + self.z);
         let b = self.x * self.y * a;
         Self::new(b, sign + self.y * self.y * a, -self.y)
@@ -1460,10 +1460,7 @@ impl {{ self_t }} {
     pub fn any_orthonormal_pair(&self) -> (Self, Self) {
         glam_assert!(self.is_normalized());
         // From https://graphics.pixar.com/library/OrthonormalB/paper.pdf
-        #[cfg(feature = "std")]
-        let sign = (1.0_{{ scalar_t }}).copysign(self.z);
-        #[cfg(not(feature = "std"))]
-        let sign = self.z.signum();
+        let sign = math::signum(self.z);
         let a = -1.0 / (sign + self.z);
         let b = self.x * self.y * a;
         (
diff --git a/src/euler.rs b/src/euler.rs
index cd658bf..bb9236b 100644
--- a/src/euler.rs
+++ b/src/euler.rs
@@ -4,11 +4,7 @@ Conversion from quaternions to Euler rotation sequences.
 From: http://bediyap.com/programming/convert-quaternion-to-euler-rotations/
 */
 
-use super::{DQuat, Quat};
-
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
+use crate::{DQuat, Quat};
 
 /// Euler rotation sequences.
 ///
@@ -62,79 +58,82 @@ pub(crate) trait EulerToQuaternion<T>: Copy {
     fn new_quat(self, u: T, v: T, w: T) -> Self::Output;
 }
 
-/// Adds a atan2 that handles the negative zero case.
-/// Basically forces positive zero in the x-argument for atan2.
-pub(crate) trait Atan2Fixed<T = Self> {
-    fn atan2_fixed(self, other: T) -> T;
-}
-
-impl Atan2Fixed for f32 {
-    fn atan2_fixed(self, other: f32) -> f32 {
-        self.atan2(if other == 0.0f32 { 0.0f32 } else { other })
-    }
-}
-impl Atan2Fixed for f64 {
-    fn atan2_fixed(self, other: f64) -> f64 {
-        self.atan2(if other == 0.0f64 { 0.0f64 } else { other })
-    }
-}
-
 macro_rules! impl_from_quat {
-    ($t:ty, $quat:ident) => {
+    ($t:ident, $quat:ident) => {
         impl EulerFromQuaternion<$quat> for EulerRot {
             type Output = $t;
             fn first(self, q: $quat) -> $t {
+                use crate::$t::math;
                 use EulerRot::*;
                 match self {
-                    ZYX => (2.0 * (q.x * q.y + q.w * q.z))
-                        .atan2(q.w * q.w + q.x * q.x - q.y * q.y - q.z * q.z),
-                    ZXY => (-2.0 * (q.x * q.y - q.w * q.z))
-                        .atan2(q.w * q.w - q.x * q.x + q.y * q.y - q.z * q.z),
-                    YXZ => (2.0 * (q.x * q.z + q.w * q.y))
-                        .atan2(q.w * q.w - q.x * q.x - q.y * q.y + q.z * q.z),
-                    YZX => (-2.0 * (q.x * q.z - q.w * q.y))
-                        .atan2(q.w * q.w + q.x * q.x - q.y * q.y - q.z * q.z),
-                    XYZ => (-2.0 * (q.y * q.z - q.w * q.x))
-                        .atan2(q.w * q.w - q.x * q.x - q.y * q.y + q.z * q.z),
-                    XZY => (2.0 * (q.y * q.z + q.w * q.x))
-                        .atan2(q.w * q.w - q.x * q.x + q.y * q.y - q.z * q.z),
+                    ZYX => math::atan2(
+                        2.0 * (q.x * q.y + q.w * q.z),
+                        q.w * q.w + q.x * q.x - q.y * q.y - q.z * q.z,
+                    ),
+                    ZXY => math::atan2(
+                        -2.0 * (q.x * q.y - q.w * q.z),
+                        q.w * q.w - q.x * q.x + q.y * q.y - q.z * q.z,
+                    ),
+                    YXZ => math::atan2(
+                        2.0 * (q.x * q.z + q.w * q.y),
+                        q.w * q.w - q.x * q.x - q.y * q.y + q.z * q.z,
+                    ),
+                    YZX => math::atan2(
+                        -2.0 * (q.x * q.z - q.w * q.y),
+                        q.w * q.w + q.x * q.x - q.y * q.y - q.z * q.z,
+                    ),
+                    XYZ => math::atan2(
+                        -2.0 * (q.y * q.z - q.w * q.x),
+                        q.w * q.w - q.x * q.x - q.y * q.y + q.z * q.z,
+                    ),
+                    XZY => math::atan2(
+                        2.0 * (q.y * q.z + q.w * q.x),
+                        q.w * q.w - q.x * q.x + q.y * q.y - q.z * q.z,
+                    ),
                 }
             }
 
             fn second(self, q: $quat) -> $t {
+                use crate::$t::math;
                 use EulerRot::*;
-
-                /// Clamp number to range [-1,1](-1,1) for asin() and acos(), else NaN is possible.
-                #[inline(always)]
-                fn arc_clamp(val: $t) -> $t {
-                    <$t>::min(<$t>::max(val, -1.0), 1.0)
-                }
-
                 match self {
-                    ZYX => arc_clamp(-2.0 * (q.x * q.z - q.w * q.y)).asin(),
-                    ZXY => arc_clamp(2.0 * (q.y * q.z + q.w * q.x)).asin(),
-                    YXZ => arc_clamp(-2.0 * (q.y * q.z - q.w * q.x)).asin(),
-                    YZX => arc_clamp(2.0 * (q.x * q.y + q.w * q.z)).asin(),
-                    XYZ => arc_clamp(2.0 * (q.x * q.z + q.w * q.y)).asin(),
-                    XZY => arc_clamp(-2.0 * (q.x * q.y - q.w * q.z)).asin(),
+                    ZYX => math::asin_clamped(-2.0 * (q.x * q.z - q.w * q.y)),
+                    ZXY => math::asin_clamped(2.0 * (q.y * q.z + q.w * q.x)),
+                    YXZ => math::asin_clamped(-2.0 * (q.y * q.z - q.w * q.x)),
+                    YZX => math::asin_clamped(2.0 * (q.x * q.y + q.w * q.z)),
+                    XYZ => math::asin_clamped(2.0 * (q.x * q.z + q.w * q.y)),
+                    XZY => math::asin_clamped(-2.0 * (q.x * q.y - q.w * q.z)),
                 }
             }
 
             fn third(self, q: $quat) -> $t {
+                use crate::$t::math;
                 use EulerRot::*;
                 match self {
-                    ZYX => (2.0 * (q.y * q.z + q.w * q.x))
-                        .atan2(q.w * q.w - q.x * q.x - q.y * q.y + q.z * q.z),
-                    ZXY => (-2.0 * (q.x * q.z - q.w * q.y))
-                        .atan2(q.w * q.w - q.x * q.x - q.y * q.y + q.z * q.z),
-                    YXZ => (2.0 * (q.x * q.y + q.w * q.z))
-                        .atan2(q.w * q.w - q.x * q.x + q.y * q.y - q.z * q.z),
-                    YZX => (-2.0 * (q.y * q.z - q.w * q.x))
-                        .atan2(q.w * q.w - q.x * q.x + q.y * q.y - q.z * q.z),
-                    XYZ => (-2.0 * (q.x * q.y - q.w * q.z))
-                        .atan2(q.w * q.w + q.x * q.x - q.y * q.y - q.z * q.z),
-                    XZY => (2.0 * (q.x * q.z + q.w * q.y))
-                        .atan2(q.w * q.w + q.x * q.x - q.y * q.y - q.z * q.z),
+                    ZYX => math::atan2(
+                        2.0 * (q.y * q.z + q.w * q.x),
+                        q.w * q.w - q.x * q.x - q.y * q.y + q.z * q.z,
+                    ),
+                    ZXY => math::atan2(
+                        -2.0 * (q.x * q.z - q.w * q.y),
+                        q.w * q.w - q.x * q.x - q.y * q.y + q.z * q.z,
+                    ),
+                    YXZ => math::atan2(
+                        2.0 * (q.x * q.y + q.w * q.z),
+                        q.w * q.w - q.x * q.x + q.y * q.y - q.z * q.z,
+                    ),
+                    YZX => math::atan2(
+                        -2.0 * (q.y * q.z - q.w * q.x),
+                        q.w * q.w - q.x * q.x + q.y * q.y - q.z * q.z,
+                    ),
+                    XYZ => math::atan2(
+                        -2.0 * (q.x * q.y - q.w * q.z),
+                        q.w * q.w + q.x * q.x - q.y * q.y - q.z * q.z,
+                    ),
+                    XZY => math::atan2(
+                        2.0 * (q.x * q.z + q.w * q.y),
+                        q.w * q.w + q.x * q.x - q.y * q.y - q.z * q.z,
+                    ),
                 }
             }
         }
diff --git a/src/f32.rs b/src/f32.rs
index 2076bd4..58fc417 100644
--- a/src/f32.rs
+++ b/src/f32.rs
@@ -1,6 +1,7 @@
 mod affine2;
 mod affine3a;
 mod mat3;
+pub(crate) mod math;
 mod vec2;
 mod vec3;
 
diff --git a/src/f32/affine3a.rs b/src/f32/affine3a.rs
index 916b5f7..0a1f8db 100644
--- a/src/f32/affine3a.rs
+++ b/src/f32/affine3a.rs
@@ -262,15 +262,12 @@ impl Affine3A {
     /// vector contains any zero elements when `glam_assert` is enabled.
     #[inline]
     pub fn to_scale_rotation_translation(&self) -> (Vec3, Quat, Vec3) {
-        #[cfg(feature = "libm")]
-        #[allow(unused_imports)]
-        use num_traits::Float;
-
+        use crate::f32::math;
         let det = self.matrix3.determinant();
         glam_assert!(det != 0.0);
 
         let scale = Vec3::new(
-            self.matrix3.x_axis.length() * det.signum(),
+            self.matrix3.x_axis.length() * math::signum(det),
             self.matrix3.y_axis.length(),
             self.matrix3.z_axis.length(),
         );
diff --git a/src/f32/coresimd/mat2.rs b/src/f32/coresimd/mat2.rs
index 68aee47..ea9169b 100644
--- a/src/f32/coresimd/mat2.rs
+++ b/src/f32/coresimd/mat2.rs
@@ -1,6 +1,6 @@
 // Generated from mat.rs.tera template. Edit the template, not the generated file.
 
-use crate::{swizzles::*, DMat2, Mat3, Mat3A, Vec2};
+use crate::{f32::math, swizzles::*, DMat2, Mat3, Mat3A, Vec2};
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
@@ -8,10 +8,6 @@ use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign};
 
 use core::simd::{Which::*, *};
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 2x2 matrix from column vectors.
 #[inline(always)]
 pub const fn mat2(x_axis: Vec2, y_axis: Vec2) -> Mat2 {
@@ -90,14 +86,14 @@ impl Mat2 {
     /// `angle` (in radians).
     #[inline]
     pub fn from_scale_angle(scale: Vec2, angle: f32) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::new(cos * scale.x, sin * scale.x, -sin * scale.y, cos * scale.y)
     }
 
     /// Creates a 2x2 matrix containing a rotation of `angle` (in radians).
     #[inline]
     pub fn from_angle(angle: f32) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::new(cos, sin, -sin, cos)
     }
 
diff --git a/src/f32/coresimd/mat3a.rs b/src/f32/coresimd/mat3a.rs
index bda6d27..97a83cf 100644
--- a/src/f32/coresimd/mat3a.rs
+++ b/src/f32/coresimd/mat3a.rs
@@ -1,6 +1,6 @@
 // Generated from mat.rs.tera template. Edit the template, not the generated file.
 
-use crate::{swizzles::*, DMat3, EulerRot, Mat2, Mat3, Mat4, Quat, Vec2, Vec3, Vec3A};
+use crate::{f32::math, swizzles::*, DMat3, EulerRot, Mat2, Mat3, Mat4, Quat, Vec2, Vec3, Vec3A};
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
@@ -8,10 +8,6 @@ use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign};
 
 use core::simd::{Which::*, *};
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 3x3 matrix from column vectors.
 #[inline(always)]
 pub const fn mat3a(x_axis: Vec3A, y_axis: Vec3A, z_axis: Vec3A) -> Mat3A {
@@ -188,7 +184,7 @@ impl Mat3A {
     pub fn from_axis_angle(axis: Vec3, angle: f32) -> Self {
         glam_assert!(axis.is_normalized());
 
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         let (xsin, ysin, zsin) = axis.mul(sin).into();
         let (x, y, z) = axis.into();
         let (x2, y2, z2) = axis.mul(axis).into();
@@ -214,7 +210,7 @@ impl Mat3A {
     /// Creates a 3D rotation matrix from `angle` (in radians) around the x axis.
     #[inline]
     pub fn from_rotation_x(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec3A::X,
             Vec3A::new(0.0, cosa, sina),
@@ -225,7 +221,7 @@ impl Mat3A {
     /// Creates a 3D rotation matrix from `angle` (in radians) around the y axis.
     #[inline]
     pub fn from_rotation_y(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec3A::new(cosa, 0.0, -sina),
             Vec3A::Y,
@@ -236,7 +232,7 @@ impl Mat3A {
     /// Creates a 3D rotation matrix from `angle` (in radians) around the z axis.
     #[inline]
     pub fn from_rotation_z(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec3A::new(cosa, sina, 0.0),
             Vec3A::new(-sina, cosa, 0.0),
@@ -264,7 +260,7 @@ impl Mat3A {
     /// [`Self::transform_point2()`] and [`Self::transform_vector2()`].
     #[inline]
     pub fn from_angle(angle: f32) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::from_cols(
             Vec3A::new(cos, sin, 0.0),
             Vec3A::new(-sin, cos, 0.0),
@@ -279,7 +275,7 @@ impl Mat3A {
     /// [`Self::transform_point2()`] and [`Self::transform_vector2()`].
     #[inline]
     pub fn from_scale_angle_translation(scale: Vec2, angle: f32, translation: Vec2) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::from_cols(
             Vec3A::new(cos * scale.x, sin * scale.x, 0.0),
             Vec3A::new(-sin * scale.y, cos * scale.y, 0.0),
diff --git a/src/f32/coresimd/mat4.rs b/src/f32/coresimd/mat4.rs
index 26239ba..9fe5445 100644
--- a/src/f32/coresimd/mat4.rs
+++ b/src/f32/coresimd/mat4.rs
@@ -1,6 +1,8 @@
 // Generated from mat.rs.tera template. Edit the template, not the generated file.
 
-use crate::{coresimd::*, swizzles::*, DMat4, EulerRot, Mat3, Mat3A, Quat, Vec3, Vec3A, Vec4};
+use crate::{
+    coresimd::*, f32::math, swizzles::*, DMat4, EulerRot, Mat3, Mat3A, Quat, Vec3, Vec3A, Vec4,
+};
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
@@ -8,10 +10,6 @@ use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign};
 
 use core::simd::{Which::*, *};
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 4x4 matrix from column vectors.
 #[inline(always)]
 pub const fn mat4(x_axis: Vec4, y_axis: Vec4, z_axis: Vec4, w_axis: Vec4) -> Mat4 {
@@ -238,7 +236,7 @@ impl Mat4 {
         glam_assert!(det != 0.0);
 
         let scale = Vec3::new(
-            self.x_axis.length() * det.signum(),
+            self.x_axis.length() * math::signum(det),
             self.y_axis.length(),
             self.z_axis.length(),
         );
@@ -329,7 +327,7 @@ impl Mat4 {
     pub fn from_axis_angle(axis: Vec3, angle: f32) -> Self {
         glam_assert!(axis.is_normalized());
 
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         let axis_sin = axis.mul(sin);
         let axis_sq = axis.mul(axis);
         let omc = 1.0 - cos;
@@ -377,7 +375,7 @@ impl Mat4 {
     /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
     #[inline]
     pub fn from_rotation_x(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec4::X,
             Vec4::new(0.0, cosa, sina, 0.0),
@@ -393,7 +391,7 @@ impl Mat4 {
     /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
     #[inline]
     pub fn from_rotation_y(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec4::new(cosa, 0.0, -sina, 0.0),
             Vec4::Y,
@@ -409,7 +407,7 @@ impl Mat4 {
     /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
     #[inline]
     pub fn from_rotation_z(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec4::new(cosa, sina, 0.0, 0.0),
             Vec4::new(-sina, cosa, 0.0, 0.0),
@@ -952,7 +950,7 @@ impl Mat4 {
         z_far: f32,
     ) -> Self {
         let inv_length = 1.0 / (z_near - z_far);
-        let f = 1.0 / (0.5 * fov_y_radians).tan();
+        let f = 1.0 / math::tan(0.5 * fov_y_radians);
         let a = f / aspect_ratio;
         let b = (z_near + z_far) * inv_length;
         let c = (2.0 * z_near * z_far) * inv_length;
@@ -973,7 +971,7 @@ impl Mat4 {
     #[inline]
     pub fn perspective_lh(fov_y_radians: f32, aspect_ratio: f32, z_near: f32, z_far: f32) -> Self {
         glam_assert!(z_near > 0.0 && z_far > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         let r = z_far / (z_far - z_near);
@@ -994,7 +992,7 @@ impl Mat4 {
     #[inline]
     pub fn perspective_rh(fov_y_radians: f32, aspect_ratio: f32, z_near: f32, z_far: f32) -> Self {
         glam_assert!(z_near > 0.0 && z_far > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         let r = z_far / (z_near - z_far);
@@ -1014,7 +1012,7 @@ impl Mat4 {
     #[inline]
     pub fn perspective_infinite_lh(fov_y_radians: f32, aspect_ratio: f32, z_near: f32) -> Self {
         glam_assert!(z_near > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         Self::from_cols(
@@ -1037,7 +1035,7 @@ impl Mat4 {
         z_near: f32,
     ) -> Self {
         glam_assert!(z_near > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         Self::from_cols(
@@ -1053,7 +1051,7 @@ impl Mat4 {
     #[inline]
     pub fn perspective_infinite_rh(fov_y_radians: f32, aspect_ratio: f32, z_near: f32) -> Self {
         glam_assert!(z_near > 0.0);
-        let f = 1.0 / (0.5 * fov_y_radians).tan();
+        let f = 1.0 / math::tan(0.5 * fov_y_radians);
         Self::from_cols(
             Vec4::new(f / aspect_ratio, 0.0, 0.0, 0.0),
             Vec4::new(0.0, f, 0.0, 0.0),
@@ -1071,7 +1069,7 @@ impl Mat4 {
         z_near: f32,
     ) -> Self {
         glam_assert!(z_near > 0.0);
-        let f = 1.0 / (0.5 * fov_y_radians).tan();
+        let f = 1.0 / math::tan(0.5 * fov_y_radians);
         Self::from_cols(
             Vec4::new(f / aspect_ratio, 0.0, 0.0, 0.0),
             Vec4::new(0.0, f, 0.0, 0.0),
diff --git a/src/f32/coresimd/quat.rs b/src/f32/coresimd/quat.rs
index 65454fd..6f5b67c 100644
--- a/src/f32/coresimd/quat.rs
+++ b/src/f32/coresimd/quat.rs
@@ -3,13 +3,10 @@
 use crate::{
     coresimd::*,
     euler::{EulerFromQuaternion, EulerRot, EulerToQuaternion},
-    DQuat, FloatEx, Mat3, Mat3A, Mat4, Vec2, Vec3, Vec3A, Vec4,
+    f32::math,
+    DQuat, Mat3, Mat3A, Mat4, Vec2, Vec3, Vec3A, Vec4,
 };
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 use core::simd::*;
 
 #[cfg(not(target_arch = "spirv"))]
@@ -124,7 +121,7 @@ impl Quat {
     #[inline]
     pub fn from_axis_angle(axis: Vec3, angle: f32) -> Self {
         glam_assert!(axis.is_normalized());
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         let v = axis * s;
         Self::from_xyzw(v.x, v.y, v.z, c)
     }
@@ -145,21 +142,21 @@ impl Quat {
     /// Creates a quaternion from the `angle` (in radians) around the x axis.
     #[inline]
     pub fn from_rotation_x(angle: f32) -> Self {
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         Self::from_xyzw(s, 0.0, 0.0, c)
     }
 
     /// Creates a quaternion from the `angle` (in radians) around the y axis.
     #[inline]
     pub fn from_rotation_y(angle: f32) -> Self {
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         Self::from_xyzw(0.0, s, 0.0, c)
     }
 
     /// Creates a quaternion from the `angle` (in radians) around the z axis.
     #[inline]
     pub fn from_rotation_z(angle: f32) -> Self {
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         Self::from_xyzw(0.0, 0.0, s, c)
     }
 
@@ -183,7 +180,7 @@ impl Quat {
             if dif10 <= 0.0 {
                 // x^2 >= y^2
                 let four_xsq = omm22 - dif10;
-                let inv4x = 0.5 / four_xsq.sqrt();
+                let inv4x = 0.5 / math::sqrt(four_xsq);
                 Self::from_xyzw(
                     four_xsq * inv4x,
                     (m01 + m10) * inv4x,
@@ -193,7 +190,7 @@ impl Quat {
             } else {
                 // y^2 >= x^2
                 let four_ysq = omm22 + dif10;
-                let inv4y = 0.5 / four_ysq.sqrt();
+                let inv4y = 0.5 / math::sqrt(four_ysq);
                 Self::from_xyzw(
                     (m01 + m10) * inv4y,
                     four_ysq * inv4y,
@@ -208,7 +205,7 @@ impl Quat {
             if sum10 <= 0.0 {
                 // z^2 >= w^2
                 let four_zsq = opm22 - sum10;
-                let inv4z = 0.5 / four_zsq.sqrt();
+                let inv4z = 0.5 / math::sqrt(four_zsq);
                 Self::from_xyzw(
                     (m02 + m20) * inv4z,
                     (m12 + m21) * inv4z,
@@ -218,7 +215,7 @@ impl Quat {
             } else {
                 // w^2 >= z^2
                 let four_wsq = opm22 + sum10;
-                let inv4w = 0.5 / four_wsq.sqrt();
+                let inv4w = 0.5 / math::sqrt(four_wsq);
                 Self::from_xyzw(
                     (m12 - m21) * inv4w,
                     (m20 - m02) * inv4w,
@@ -338,7 +335,7 @@ impl Quat {
             let z = from.x * to.y - to.x * from.y;
             let w = 1.0 + dot;
             // calculate length with x=0 and y=0 to normalize
-            let len_rcp = 1.0 / (z * z + w * w).sqrt();
+            let len_rcp = 1.0 / math::sqrt(z * z + w * w);
             Self::from_xyzw(0.0, 0.0, z * len_rcp, w * len_rcp)
         }
     }
@@ -350,7 +347,7 @@ impl Quat {
         let v = Vec3::new(self.x, self.y, self.z);
         let length = v.length();
         if length >= EPSILON {
-            let angle = 2.0 * f32::atan2(length, self.w);
+            let angle = 2.0 * math::atan2(length, self.w);
             let axis = v / length;
             (axis, angle)
         } else {
@@ -490,7 +487,7 @@ impl Quat {
         // If the quat.w is close to -1.0, the angle will be near 2*PI which is close to
         // a negative 0 rotation. By forcing quat.w to be positive, we'll end up with
         // the shortest path.
-        let positive_w_angle = self.w.abs().acos_approx() * 2.0;
+        let positive_w_angle = math::acos_approx(math::abs(self.w)) * 2.0;
         positive_w_angle < threshold_angle
     }
 
@@ -505,7 +502,7 @@ impl Quat {
     #[inline]
     pub fn angle_between(self, rhs: Self) -> f32 {
         glam_assert!(self.is_normalized() && rhs.is_normalized());
-        self.dot(rhs).abs().acos_approx() * 2.0
+        math::acos_approx(math::abs(self.dot(rhs))) * 2.0
     }
 
     /// Returns true if the absolute difference of all elements between `self` and `rhs`
@@ -581,11 +578,11 @@ impl Quat {
             // assumes lerp returns a normalized quaternion
             self.lerp(end, s)
         } else {
-            let theta = dot.acos_approx();
+            let theta = math::acos_approx(dot);
 
-            let x = (theta * (1.0 - s)).sin();
-            let y = (theta * s).sin();
-            let z = theta.sin();
+            let x = math::sin(theta * (1.0 - s));
+            let y = math::sin(theta * s);
+            let z = math::sin(theta);
             let w = 0.0;
             let tmp = f32x4::from_array([x, y, z, w]);
 
diff --git a/src/f32/coresimd/vec3a.rs b/src/f32/coresimd/vec3a.rs
index eff4464..83489a8 100644
--- a/src/f32/coresimd/vec3a.rs
+++ b/src/f32/coresimd/vec3a.rs
@@ -1,6 +1,6 @@
 // Generated from vec.rs.tera template. Edit the template, not the generated file.
 
-use crate::{coresimd::*, BVec3A, Vec2, Vec3, Vec4};
+use crate::{coresimd::*, f32::math, BVec3A, Vec2, Vec3, Vec4};
 
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
@@ -10,10 +10,6 @@ use core::{f32, ops::*};
 use core::simd::*;
 use std::simd::StdFloat;
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 3-dimensional vector.
 #[inline(always)]
 pub const fn vec3a(x: f32, y: f32, z: f32) -> Vec3A {
@@ -430,7 +426,7 @@ impl Vec3A {
     #[inline]
     pub fn is_normalized(self) -> bool {
         // TODO: do something with epsilon
-        (self.length_squared() - 1.0).abs() <= 1e-4
+        math::abs(self.length_squared() - 1.0) <= 1e-4
     }
 
     /// Returns the vector projection of `self` onto `rhs`.
@@ -528,13 +524,17 @@ impl Vec3A {
     /// `self`.
     #[inline]
     pub fn exp(self) -> Self {
-        Self::new(self.x.exp(), self.y.exp(), self.z.exp())
+        Self::new(math::exp(self.x), math::exp(self.y), math::exp(self.z))
     }
 
     /// Returns a vector containing each element of `self` raised to the power of `n`.
     #[inline]
     pub fn powf(self, n: f32) -> Self {
-        Self::new(self.x.powf(n), self.y.powf(n), self.z.powf(n))
+        Self::new(
+            math::powf(self.x, n),
+            math::powf(self.y, n),
+            math::powf(self.z, n),
+        )
     }
 
     /// Returns a vector containing the reciprocal `1.0/n` of each element of `self`.
@@ -578,9 +578,9 @@ impl Vec3A {
         glam_assert!(min <= max);
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -590,7 +590,7 @@ impl Vec3A {
     pub fn clamp_length_max(self, max: f32) -> Self {
         let length_sq = self.length_squared();
         if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -600,7 +600,7 @@ impl Vec3A {
     pub fn clamp_length_min(self, min: f32) -> Self {
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -623,10 +623,10 @@ impl Vec3A {
     /// The input vectors do not need to be unit length however they must be non-zero.
     #[inline]
     pub fn angle_between(self, rhs: Self) -> f32 {
-        use crate::FloatEx;
-        self.dot(rhs)
-            .div(self.length_squared().mul(rhs.length_squared()).sqrt())
-            .acos_approx()
+        math::acos_approx(
+            self.dot(rhs)
+                .div(math::sqrt(self.length_squared().mul(rhs.length_squared()))),
+        )
     }
 
     /// Returns some vector that is orthogonal to the given one.
@@ -638,7 +638,7 @@ impl Vec3A {
     #[inline]
     pub fn any_orthogonal_vector(&self) -> Self {
         // This can probably be optimized
-        if self.x.abs() > self.y.abs() {
+        if math::abs(self.x) > math::abs(self.y) {
             Self::new(-self.z, 0.0, self.x) // self.cross(Self::Y)
         } else {
             Self::new(0.0, self.z, -self.y) // self.cross(Self::X)
@@ -655,10 +655,7 @@ impl Vec3A {
     pub fn any_orthonormal_vector(&self) -> Self {
         glam_assert!(self.is_normalized());
         // From https://graphics.pixar.com/library/OrthonormalB/paper.pdf
-        #[cfg(feature = "std")]
-        let sign = (1.0_f32).copysign(self.z);
-        #[cfg(not(feature = "std"))]
-        let sign = self.z.signum();
+        let sign = math::signum(self.z);
         let a = -1.0 / (sign + self.z);
         let b = self.x * self.y * a;
         Self::new(b, sign + self.y * self.y * a, -self.y)
@@ -674,10 +671,7 @@ impl Vec3A {
     pub fn any_orthonormal_pair(&self) -> (Self, Self) {
         glam_assert!(self.is_normalized());
         // From https://graphics.pixar.com/library/OrthonormalB/paper.pdf
-        #[cfg(feature = "std")]
-        let sign = (1.0_f32).copysign(self.z);
-        #[cfg(not(feature = "std"))]
-        let sign = self.z.signum();
+        let sign = math::signum(self.z);
         let a = -1.0 / (sign + self.z);
         let b = self.x * self.y * a;
         (
diff --git a/src/f32/coresimd/vec4.rs b/src/f32/coresimd/vec4.rs
index 3a476fd..e78563a 100644
--- a/src/f32/coresimd/vec4.rs
+++ b/src/f32/coresimd/vec4.rs
@@ -1,6 +1,6 @@
 // Generated from vec.rs.tera template. Edit the template, not the generated file.
 
-use crate::{coresimd::*, BVec4A, Vec2, Vec3, Vec3A};
+use crate::{coresimd::*, f32::math, BVec4A, Vec2, Vec3, Vec3A};
 
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
@@ -10,10 +10,6 @@ use core::{f32, ops::*};
 use core::simd::*;
 use std::simd::StdFloat;
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 4-dimensional vector.
 #[inline(always)]
 pub const fn vec4(x: f32, y: f32, z: f32, w: f32) -> Vec4 {
@@ -403,7 +399,7 @@ impl Vec4 {
     #[inline]
     pub fn is_normalized(self) -> bool {
         // TODO: do something with epsilon
-        (self.length_squared() - 1.0).abs() <= 1e-4
+        math::abs(self.length_squared() - 1.0) <= 1e-4
     }
 
     /// Returns the vector projection of `self` onto `rhs`.
@@ -501,17 +497,22 @@ impl Vec4 {
     /// `self`.
     #[inline]
     pub fn exp(self) -> Self {
-        Self::new(self.x.exp(), self.y.exp(), self.z.exp(), self.w.exp())
+        Self::new(
+            math::exp(self.x),
+            math::exp(self.y),
+            math::exp(self.z),
+            math::exp(self.w),
+        )
     }
 
     /// Returns a vector containing each element of `self` raised to the power of `n`.
     #[inline]
     pub fn powf(self, n: f32) -> Self {
         Self::new(
-            self.x.powf(n),
-            self.y.powf(n),
-            self.z.powf(n),
-            self.w.powf(n),
+            math::powf(self.x, n),
+            math::powf(self.y, n),
+            math::powf(self.z, n),
+            math::powf(self.w, n),
         )
     }
 
@@ -556,9 +557,9 @@ impl Vec4 {
         glam_assert!(min <= max);
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -568,7 +569,7 @@ impl Vec4 {
     pub fn clamp_length_max(self, max: f32) -> Self {
         let length_sq = self.length_squared();
         if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -578,7 +579,7 @@ impl Vec4 {
     pub fn clamp_length_min(self, min: f32) -> Self {
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else {
             self
         }
diff --git a/src/f32/mat3.rs b/src/f32/mat3.rs
index 6c2ed4d..7f1fe00 100644
--- a/src/f32/mat3.rs
+++ b/src/f32/mat3.rs
@@ -1,15 +1,11 @@
 // Generated from mat.rs.tera template. Edit the template, not the generated file.
 
-use crate::{swizzles::*, DMat3, EulerRot, Mat2, Mat3A, Mat4, Quat, Vec2, Vec3, Vec3A};
+use crate::{f32::math, swizzles::*, DMat3, EulerRot, Mat2, Mat3A, Mat4, Quat, Vec2, Vec3, Vec3A};
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
 use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign};
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 3x3 matrix from column vectors.
 #[inline(always)]
 pub const fn mat3(x_axis: Vec3, y_axis: Vec3, z_axis: Vec3) -> Mat3 {
@@ -189,7 +185,7 @@ impl Mat3 {
     pub fn from_axis_angle(axis: Vec3, angle: f32) -> Self {
         glam_assert!(axis.is_normalized());
 
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         let (xsin, ysin, zsin) = axis.mul(sin).into();
         let (x, y, z) = axis.into();
         let (x2, y2, z2) = axis.mul(axis).into();
@@ -215,7 +211,7 @@ impl Mat3 {
     /// Creates a 3D rotation matrix from `angle` (in radians) around the x axis.
     #[inline]
     pub fn from_rotation_x(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec3::X,
             Vec3::new(0.0, cosa, sina),
@@ -226,7 +222,7 @@ impl Mat3 {
     /// Creates a 3D rotation matrix from `angle` (in radians) around the y axis.
     #[inline]
     pub fn from_rotation_y(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec3::new(cosa, 0.0, -sina),
             Vec3::Y,
@@ -237,7 +233,7 @@ impl Mat3 {
     /// Creates a 3D rotation matrix from `angle` (in radians) around the z axis.
     #[inline]
     pub fn from_rotation_z(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec3::new(cosa, sina, 0.0),
             Vec3::new(-sina, cosa, 0.0),
@@ -265,7 +261,7 @@ impl Mat3 {
     /// [`Self::transform_point2()`] and [`Self::transform_vector2()`].
     #[inline]
     pub fn from_angle(angle: f32) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::from_cols(Vec3::new(cos, sin, 0.0), Vec3::new(-sin, cos, 0.0), Vec3::Z)
     }
 
@@ -276,7 +272,7 @@ impl Mat3 {
     /// [`Self::transform_point2()`] and [`Self::transform_vector2()`].
     #[inline]
     pub fn from_scale_angle_translation(scale: Vec2, angle: f32, translation: Vec2) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::from_cols(
             Vec3::new(cos * scale.x, sin * scale.x, 0.0),
             Vec3::new(-sin * scale.y, cos * scale.y, 0.0),
diff --git a/src/f32/math.rs b/src/f32/math.rs
new file mode 100644
index 0000000..6a71826
--- /dev/null
+++ b/src/f32/math.rs
@@ -0,0 +1,217 @@
+/// Returns a very close approximation of `self.clamp(-1.0, 1.0).acos()`.
+#[inline]
+fn acos_approx_f32(v: f32) -> f32 {
+    // Based on https://github.com/microsoft/DirectXMath `XMScalarAcos`
+    // Clamp input to [-1,1].
+    let nonnegative = v >= 0.0;
+    let x = abs(v);
+    let mut omx = 1.0 - x;
+    if omx < 0.0 {
+        omx = 0.0;
+    }
+    let root = sqrt(omx);
+
+    // 7-degree minimax approximation
+    #[allow(clippy::approx_constant)]
+    let mut result =
+        ((((((-0.001_262_491_1 * x + 0.006_670_09) * x - 0.017_088_126) * x + 0.030_891_88) * x
+            - 0.050_174_303)
+            * x
+            + 0.088_978_99)
+            * x
+            - 0.214_598_8)
+            * x
+            + 1.570_796_3;
+    result *= root;
+
+    // acos(x) = pi - acos(-x) when x < 0
+    if nonnegative {
+        result
+    } else {
+        core::f32::consts::PI - result
+    }
+}
+
+#[cfg(feature = "libm")]
+mod libm_math {
+    #[inline(always)]
+    pub(crate) fn abs(f: f32) -> f32 {
+        if f.is_sign_positive() {
+            return f;
+        }
+        if f.is_sign_negative() {
+            return -f;
+        }
+        f32::NAN
+    }
+
+    #[inline(always)]
+    pub(crate) fn acos_approx(f: f32) -> f32 {
+        super::acos_approx_f32(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn asin_clamped(f: f32) -> f32 {
+        libm::asinf(f.clamp(-1.0, 1.0))
+    }
+
+    #[inline(always)]
+    pub(crate) fn atan2(f: f32, other: f32) -> f32 {
+        libm::atan2f(f, other)
+    }
+
+    #[allow(unused)]
+    #[inline(always)]
+    pub(crate) fn sin(f: f32) -> f32 {
+        libm::sinf(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn sin_cos(f: f32) -> (f32, f32) {
+        libm::sincosf(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn tan(f: f32) -> f32 {
+        libm::tanf(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn sqrt(f: f32) -> f32 {
+        libm::sqrtf(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn copysign(f: f32, sign: f32) -> f32 {
+        libm::copysignf(f, sign)
+    }
+
+    #[inline(always)]
+    pub(crate) fn signum(f: f32) -> f32 {
+        if f.is_nan() {
+            f32::NAN
+        } else {
+            copysign(1.0, f)
+        }
+    }
+
+    #[inline(always)]
+    pub(crate) fn round(f: f32) -> f32 {
+        libm::roundf(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn ceil(f: f32) -> f32 {
+        libm::ceilf(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn floor(f: f32) -> f32 {
+        libm::floorf(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn exp(f: f32) -> f32 {
+        libm::expf(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn powf(f: f32, n: f32) -> f32 {
+        libm::powf(f, n)
+    }
+
+    #[inline(always)]
+    pub(crate) fn mul_add(a: f32, b: f32, c: f32) -> f32 {
+        libm::fmaf(a, b, c)
+    }
+}
+
+#[cfg(not(feature = "libm"))]
+mod std_math {
+    #[inline(always)]
+    pub(crate) fn abs(f: f32) -> f32 {
+        f32::abs(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn acos_approx(f: f32) -> f32 {
+        super::acos_approx_f32(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn asin_clamped(f: f32) -> f32 {
+        f32::asin(f32::clamp(f, -1.0, 1.0))
+    }
+
+    #[inline(always)]
+    pub(crate) fn atan2(f: f32, other: f32) -> f32 {
+        f32::atan2(f, other)
+    }
+
+    #[allow(unused)]
+    #[inline(always)]
+    pub(crate) fn sin(f: f32) -> f32 {
+        f32::sin(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn sin_cos(f: f32) -> (f32, f32) {
+        f32::sin_cos(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn tan(f: f32) -> f32 {
+        f32::tan(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn sqrt(f: f32) -> f32 {
+        f32::sqrt(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn copysign(f: f32, sign: f32) -> f32 {
+        f32::copysign(f, sign)
+    }
+
+    #[inline(always)]
+    pub(crate) fn signum(f: f32) -> f32 {
+        f32::signum(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn round(f: f32) -> f32 {
+        f32::round(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn ceil(f: f32) -> f32 {
+        f32::ceil(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn floor(f: f32) -> f32 {
+        f32::floor(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn exp(f: f32) -> f32 {
+        f32::exp(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn powf(f: f32, n: f32) -> f32 {
+        f32::powf(f, n)
+    }
+
+    #[inline(always)]
+    pub(crate) fn mul_add(a: f32, b: f32, c: f32) -> f32 {
+        f32::mul_add(a, b, c)
+    }
+}
+
+#[cfg(feature = "libm")]
+pub(crate) use libm_math::*;
+
+#[cfg(not(feature = "libm"))]
+pub(crate) use std_math::*;
diff --git a/src/f32/scalar/mat2.rs b/src/f32/scalar/mat2.rs
index 9c4c253..604deda 100644
--- a/src/f32/scalar/mat2.rs
+++ b/src/f32/scalar/mat2.rs
@@ -1,15 +1,11 @@
 // Generated from mat.rs.tera template. Edit the template, not the generated file.
 
-use crate::{swizzles::*, DMat2, Mat3, Mat3A, Vec2};
+use crate::{f32::math, swizzles::*, DMat2, Mat3, Mat3A, Vec2};
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
 use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign};
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 2x2 matrix from column vectors.
 #[inline(always)]
 pub const fn mat2(x_axis: Vec2, y_axis: Vec2) -> Mat2 {
@@ -95,14 +91,14 @@ impl Mat2 {
     /// `angle` (in radians).
     #[inline]
     pub fn from_scale_angle(scale: Vec2, angle: f32) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::new(cos * scale.x, sin * scale.x, -sin * scale.y, cos * scale.y)
     }
 
     /// Creates a 2x2 matrix containing a rotation of `angle` (in radians).
     #[inline]
     pub fn from_angle(angle: f32) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::new(cos, sin, -sin, cos)
     }
 
diff --git a/src/f32/scalar/mat3a.rs b/src/f32/scalar/mat3a.rs
index ea56981..9fc8193 100644
--- a/src/f32/scalar/mat3a.rs
+++ b/src/f32/scalar/mat3a.rs
@@ -1,15 +1,11 @@
 // Generated from mat.rs.tera template. Edit the template, not the generated file.
 
-use crate::{swizzles::*, DMat3, EulerRot, Mat2, Mat3, Mat4, Quat, Vec2, Vec3, Vec3A};
+use crate::{f32::math, swizzles::*, DMat3, EulerRot, Mat2, Mat3, Mat4, Quat, Vec2, Vec3, Vec3A};
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
 use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign};
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 3x3 matrix from column vectors.
 #[inline(always)]
 pub const fn mat3a(x_axis: Vec3A, y_axis: Vec3A, z_axis: Vec3A) -> Mat3A {
@@ -189,7 +185,7 @@ impl Mat3A {
     pub fn from_axis_angle(axis: Vec3, angle: f32) -> Self {
         glam_assert!(axis.is_normalized());
 
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         let (xsin, ysin, zsin) = axis.mul(sin).into();
         let (x, y, z) = axis.into();
         let (x2, y2, z2) = axis.mul(axis).into();
@@ -215,7 +211,7 @@ impl Mat3A {
     /// Creates a 3D rotation matrix from `angle` (in radians) around the x axis.
     #[inline]
     pub fn from_rotation_x(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec3A::X,
             Vec3A::new(0.0, cosa, sina),
@@ -226,7 +222,7 @@ impl Mat3A {
     /// Creates a 3D rotation matrix from `angle` (in radians) around the y axis.
     #[inline]
     pub fn from_rotation_y(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec3A::new(cosa, 0.0, -sina),
             Vec3A::Y,
@@ -237,7 +233,7 @@ impl Mat3A {
     /// Creates a 3D rotation matrix from `angle` (in radians) around the z axis.
     #[inline]
     pub fn from_rotation_z(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec3A::new(cosa, sina, 0.0),
             Vec3A::new(-sina, cosa, 0.0),
@@ -265,7 +261,7 @@ impl Mat3A {
     /// [`Self::transform_point2()`] and [`Self::transform_vector2()`].
     #[inline]
     pub fn from_angle(angle: f32) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::from_cols(
             Vec3A::new(cos, sin, 0.0),
             Vec3A::new(-sin, cos, 0.0),
@@ -280,7 +276,7 @@ impl Mat3A {
     /// [`Self::transform_point2()`] and [`Self::transform_vector2()`].
     #[inline]
     pub fn from_scale_angle_translation(scale: Vec2, angle: f32, translation: Vec2) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::from_cols(
             Vec3A::new(cos * scale.x, sin * scale.x, 0.0),
             Vec3A::new(-sin * scale.y, cos * scale.y, 0.0),
diff --git a/src/f32/scalar/mat4.rs b/src/f32/scalar/mat4.rs
index 78d7a52..c5d66c9 100644
--- a/src/f32/scalar/mat4.rs
+++ b/src/f32/scalar/mat4.rs
@@ -1,15 +1,11 @@
 // Generated from mat.rs.tera template. Edit the template, not the generated file.
 
-use crate::{swizzles::*, DMat4, EulerRot, Mat3, Mat3A, Quat, Vec3, Vec3A, Vec4};
+use crate::{f32::math, swizzles::*, DMat4, EulerRot, Mat3, Mat3A, Quat, Vec3, Vec3A, Vec4};
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
 use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign};
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 4x4 matrix from column vectors.
 #[inline(always)]
 pub const fn mat4(x_axis: Vec4, y_axis: Vec4, z_axis: Vec4, w_axis: Vec4) -> Mat4 {
@@ -251,7 +247,7 @@ impl Mat4 {
         glam_assert!(det != 0.0);
 
         let scale = Vec3::new(
-            self.x_axis.length() * det.signum(),
+            self.x_axis.length() * math::signum(det),
             self.y_axis.length(),
             self.z_axis.length(),
         );
@@ -342,7 +338,7 @@ impl Mat4 {
     pub fn from_axis_angle(axis: Vec3, angle: f32) -> Self {
         glam_assert!(axis.is_normalized());
 
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         let axis_sin = axis.mul(sin);
         let axis_sq = axis.mul(axis);
         let omc = 1.0 - cos;
@@ -390,7 +386,7 @@ impl Mat4 {
     /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
     #[inline]
     pub fn from_rotation_x(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec4::X,
             Vec4::new(0.0, cosa, sina, 0.0),
@@ -406,7 +402,7 @@ impl Mat4 {
     /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
     #[inline]
     pub fn from_rotation_y(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec4::new(cosa, 0.0, -sina, 0.0),
             Vec4::Y,
@@ -422,7 +418,7 @@ impl Mat4 {
     /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
     #[inline]
     pub fn from_rotation_z(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec4::new(cosa, sina, 0.0, 0.0),
             Vec4::new(-sina, cosa, 0.0, 0.0),
@@ -731,7 +727,7 @@ impl Mat4 {
         z_far: f32,
     ) -> Self {
         let inv_length = 1.0 / (z_near - z_far);
-        let f = 1.0 / (0.5 * fov_y_radians).tan();
+        let f = 1.0 / math::tan(0.5 * fov_y_radians);
         let a = f / aspect_ratio;
         let b = (z_near + z_far) * inv_length;
         let c = (2.0 * z_near * z_far) * inv_length;
@@ -752,7 +748,7 @@ impl Mat4 {
     #[inline]
     pub fn perspective_lh(fov_y_radians: f32, aspect_ratio: f32, z_near: f32, z_far: f32) -> Self {
         glam_assert!(z_near > 0.0 && z_far > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         let r = z_far / (z_far - z_near);
@@ -773,7 +769,7 @@ impl Mat4 {
     #[inline]
     pub fn perspective_rh(fov_y_radians: f32, aspect_ratio: f32, z_near: f32, z_far: f32) -> Self {
         glam_assert!(z_near > 0.0 && z_far > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         let r = z_far / (z_near - z_far);
@@ -793,7 +789,7 @@ impl Mat4 {
     #[inline]
     pub fn perspective_infinite_lh(fov_y_radians: f32, aspect_ratio: f32, z_near: f32) -> Self {
         glam_assert!(z_near > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         Self::from_cols(
@@ -816,7 +812,7 @@ impl Mat4 {
         z_near: f32,
     ) -> Self {
         glam_assert!(z_near > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         Self::from_cols(
@@ -832,7 +828,7 @@ impl Mat4 {
     #[inline]
     pub fn perspective_infinite_rh(fov_y_radians: f32, aspect_ratio: f32, z_near: f32) -> Self {
         glam_assert!(z_near > 0.0);
-        let f = 1.0 / (0.5 * fov_y_radians).tan();
+        let f = 1.0 / math::tan(0.5 * fov_y_radians);
         Self::from_cols(
             Vec4::new(f / aspect_ratio, 0.0, 0.0, 0.0),
             Vec4::new(0.0, f, 0.0, 0.0),
@@ -850,7 +846,7 @@ impl Mat4 {
         z_near: f32,
     ) -> Self {
         glam_assert!(z_near > 0.0);
-        let f = 1.0 / (0.5 * fov_y_radians).tan();
+        let f = 1.0 / math::tan(0.5 * fov_y_radians);
         Self::from_cols(
             Vec4::new(f / aspect_ratio, 0.0, 0.0, 0.0),
             Vec4::new(0.0, f, 0.0, 0.0),
diff --git a/src/f32/scalar/quat.rs b/src/f32/scalar/quat.rs
index 3ec4f96..718fe6f 100644
--- a/src/f32/scalar/quat.rs
+++ b/src/f32/scalar/quat.rs
@@ -2,13 +2,10 @@
 
 use crate::{
     euler::{EulerFromQuaternion, EulerRot, EulerToQuaternion},
-    DQuat, FloatEx, Mat3, Mat3A, Mat4, Vec2, Vec3, Vec3A, Vec4,
+    f32::math,
+    DQuat, Mat3, Mat3A, Mat4, Vec2, Vec3, Vec3A, Vec4,
 };
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
@@ -132,7 +129,7 @@ impl Quat {
     #[inline]
     pub fn from_axis_angle(axis: Vec3, angle: f32) -> Self {
         glam_assert!(axis.is_normalized());
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         let v = axis * s;
         Self::from_xyzw(v.x, v.y, v.z, c)
     }
@@ -153,21 +150,21 @@ impl Quat {
     /// Creates a quaternion from the `angle` (in radians) around the x axis.
     #[inline]
     pub fn from_rotation_x(angle: f32) -> Self {
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         Self::from_xyzw(s, 0.0, 0.0, c)
     }
 
     /// Creates a quaternion from the `angle` (in radians) around the y axis.
     #[inline]
     pub fn from_rotation_y(angle: f32) -> Self {
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         Self::from_xyzw(0.0, s, 0.0, c)
     }
 
     /// Creates a quaternion from the `angle` (in radians) around the z axis.
     #[inline]
     pub fn from_rotation_z(angle: f32) -> Self {
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         Self::from_xyzw(0.0, 0.0, s, c)
     }
 
@@ -191,7 +188,7 @@ impl Quat {
             if dif10 <= 0.0 {
                 // x^2 >= y^2
                 let four_xsq = omm22 - dif10;
-                let inv4x = 0.5 / four_xsq.sqrt();
+                let inv4x = 0.5 / math::sqrt(four_xsq);
                 Self::from_xyzw(
                     four_xsq * inv4x,
                     (m01 + m10) * inv4x,
@@ -201,7 +198,7 @@ impl Quat {
             } else {
                 // y^2 >= x^2
                 let four_ysq = omm22 + dif10;
-                let inv4y = 0.5 / four_ysq.sqrt();
+                let inv4y = 0.5 / math::sqrt(four_ysq);
                 Self::from_xyzw(
                     (m01 + m10) * inv4y,
                     four_ysq * inv4y,
@@ -216,7 +213,7 @@ impl Quat {
             if sum10 <= 0.0 {
                 // z^2 >= w^2
                 let four_zsq = opm22 - sum10;
-                let inv4z = 0.5 / four_zsq.sqrt();
+                let inv4z = 0.5 / math::sqrt(four_zsq);
                 Self::from_xyzw(
                     (m02 + m20) * inv4z,
                     (m12 + m21) * inv4z,
@@ -226,7 +223,7 @@ impl Quat {
             } else {
                 // w^2 >= z^2
                 let four_wsq = opm22 + sum10;
-                let inv4w = 0.5 / four_wsq.sqrt();
+                let inv4w = 0.5 / math::sqrt(four_wsq);
                 Self::from_xyzw(
                     (m12 - m21) * inv4w,
                     (m20 - m02) * inv4w,
@@ -346,7 +343,7 @@ impl Quat {
             let z = from.x * to.y - to.x * from.y;
             let w = 1.0 + dot;
             // calculate length with x=0 and y=0 to normalize
-            let len_rcp = 1.0 / (z * z + w * w).sqrt();
+            let len_rcp = 1.0 / math::sqrt(z * z + w * w);
             Self::from_xyzw(0.0, 0.0, z * len_rcp, w * len_rcp)
         }
     }
@@ -358,7 +355,7 @@ impl Quat {
         let v = Vec3::new(self.x, self.y, self.z);
         let length = v.length();
         if length >= EPSILON {
-            let angle = 2.0 * f32::atan2(length, self.w);
+            let angle = 2.0 * math::atan2(length, self.w);
             let axis = v / length;
             (axis, angle)
         } else {
@@ -502,7 +499,7 @@ impl Quat {
         // If the quat.w is close to -1.0, the angle will be near 2*PI which is close to
         // a negative 0 rotation. By forcing quat.w to be positive, we'll end up with
         // the shortest path.
-        let positive_w_angle = self.w.abs().acos_approx() * 2.0;
+        let positive_w_angle = math::acos_approx(math::abs(self.w)) * 2.0;
         positive_w_angle < threshold_angle
     }
 
@@ -517,7 +514,7 @@ impl Quat {
     #[inline]
     pub fn angle_between(self, rhs: Self) -> f32 {
         glam_assert!(self.is_normalized() && rhs.is_normalized());
-        self.dot(rhs).abs().acos_approx() * 2.0
+        math::acos_approx(math::abs(self.dot(rhs))) * 2.0
     }
 
     /// Returns true if the absolute difference of all elements between `self` and `rhs`
@@ -589,13 +586,13 @@ impl Quat {
             // assumes lerp returns a normalized quaternion
             self.lerp(end, s)
         } else {
-            let theta = dot.acos_approx();
+            let theta = math::acos_approx(dot);
 
-            let scale1 = (theta * (1.0 - s)).sin();
-            let scale2 = (theta * s).sin();
-            let theta_sin = theta.sin();
+            let scale1 = math::sin(theta * (1.0 - s));
+            let scale2 = math::sin(theta * s);
+            let theta_sin = math::sin(theta);
 
-            self.mul(scale1).add(end.mul(scale2)).mul(theta_sin.recip())
+            self.mul(scale1).add(end.mul(scale2)).mul(1.0 / theta_sin)
         }
     }
 
diff --git a/src/f32/scalar/vec3a.rs b/src/f32/scalar/vec3a.rs
index 1fd25af..c745f2b 100644
--- a/src/f32/scalar/vec3a.rs
+++ b/src/f32/scalar/vec3a.rs
@@ -1,16 +1,12 @@
 // Generated from vec.rs.tera template. Edit the template, not the generated file.
 
-use crate::{BVec3A, Vec2, Vec3, Vec4};
+use crate::{f32::math, BVec3A, Vec2, Vec3, Vec4};
 
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
 use core::{f32, ops::*};
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 3-dimensional vector.
 #[inline(always)]
 pub const fn vec3a(x: f32, y: f32, z: f32) -> Vec3A {
@@ -295,9 +291,9 @@ impl Vec3A {
     #[inline]
     pub fn abs(self) -> Self {
         Self {
-            x: self.x.abs(),
-            y: self.y.abs(),
-            z: self.z.abs(),
+            x: math::abs(self.x),
+            y: math::abs(self.y),
+            z: math::abs(self.z),
         }
     }
 
@@ -309,9 +305,9 @@ impl Vec3A {
     #[inline]
     pub fn signum(self) -> Self {
         Self {
-            x: self.x.signum(),
-            y: self.y.signum(),
-            z: self.z.signum(),
+            x: math::signum(self.x),
+            y: math::signum(self.y),
+            z: math::signum(self.z),
         }
     }
 
@@ -319,9 +315,9 @@ impl Vec3A {
     #[inline]
     pub fn copysign(self, rhs: Self) -> Self {
         Self {
-            x: self.x.copysign(rhs.x),
-            y: self.y.copysign(rhs.y),
-            z: self.z.copysign(rhs.z),
+            x: math::copysign(self.x, rhs.x),
+            y: math::copysign(self.y, rhs.y),
+            z: math::copysign(self.z, rhs.z),
         }
     }
 
@@ -361,7 +357,7 @@ impl Vec3A {
     #[doc(alias = "magnitude")]
     #[inline]
     pub fn length(self) -> f32 {
-        self.dot(self).sqrt()
+        math::sqrt(self.dot(self))
     }
 
     /// Computes the squared length of `self`.
@@ -451,7 +447,7 @@ impl Vec3A {
     #[inline]
     pub fn is_normalized(self) -> bool {
         // TODO: do something with epsilon
-        (self.length_squared() - 1.0).abs() <= 1e-4
+        math::abs(self.length_squared() - 1.0) <= 1e-4
     }
 
     /// Returns the vector projection of `self` onto `rhs`.
@@ -520,9 +516,9 @@ impl Vec3A {
     #[inline]
     pub fn round(self) -> Self {
         Self {
-            x: self.x.round(),
-            y: self.y.round(),
-            z: self.z.round(),
+            x: math::round(self.x),
+            y: math::round(self.y),
+            z: math::round(self.z),
         }
     }
 
@@ -531,9 +527,9 @@ impl Vec3A {
     #[inline]
     pub fn floor(self) -> Self {
         Self {
-            x: self.x.floor(),
-            y: self.y.floor(),
-            z: self.z.floor(),
+            x: math::floor(self.x),
+            y: math::floor(self.y),
+            z: math::floor(self.z),
         }
     }
 
@@ -542,9 +538,9 @@ impl Vec3A {
     #[inline]
     pub fn ceil(self) -> Self {
         Self {
-            x: self.x.ceil(),
-            y: self.y.ceil(),
-            z: self.z.ceil(),
+            x: math::ceil(self.x),
+            y: math::ceil(self.y),
+            z: math::ceil(self.z),
         }
     }
 
@@ -561,22 +557,26 @@ impl Vec3A {
     /// `self`.
     #[inline]
     pub fn exp(self) -> Self {
-        Self::new(self.x.exp(), self.y.exp(), self.z.exp())
+        Self::new(math::exp(self.x), math::exp(self.y), math::exp(self.z))
     }
 
     /// Returns a vector containing each element of `self` raised to the power of `n`.
     #[inline]
     pub fn powf(self, n: f32) -> Self {
-        Self::new(self.x.powf(n), self.y.powf(n), self.z.powf(n))
+        Self::new(
+            math::powf(self.x, n),
+            math::powf(self.y, n),
+            math::powf(self.z, n),
+        )
     }
 
     /// Returns a vector containing the reciprocal `1.0/n` of each element of `self`.
     #[inline]
     pub fn recip(self) -> Self {
         Self {
-            x: self.x.recip(),
-            y: self.y.recip(),
-            z: self.z.recip(),
+            x: 1.0 / self.x,
+            y: 1.0 / self.y,
+            z: 1.0 / self.z,
         }
     }
 
@@ -615,9 +615,9 @@ impl Vec3A {
         glam_assert!(min <= max);
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -627,7 +627,7 @@ impl Vec3A {
     pub fn clamp_length_max(self, max: f32) -> Self {
         let length_sq = self.length_squared();
         if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -637,7 +637,7 @@ impl Vec3A {
     pub fn clamp_length_min(self, min: f32) -> Self {
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -653,9 +653,9 @@ impl Vec3A {
     #[inline]
     pub fn mul_add(self, a: Self, b: Self) -> Self {
         Self::new(
-            self.x.mul_add(a.x, b.x),
-            self.y.mul_add(a.y, b.y),
-            self.z.mul_add(a.z, b.z),
+            math::mul_add(self.x, a.x, b.x),
+            math::mul_add(self.y, a.y, b.y),
+            math::mul_add(self.z, a.z, b.z),
         )
     }
 
@@ -664,10 +664,10 @@ impl Vec3A {
     /// The input vectors do not need to be unit length however they must be non-zero.
     #[inline]
     pub fn angle_between(self, rhs: Self) -> f32 {
-        use crate::FloatEx;
-        self.dot(rhs)
-            .div(self.length_squared().mul(rhs.length_squared()).sqrt())
-            .acos_approx()
+        math::acos_approx(
+            self.dot(rhs)
+                .div(math::sqrt(self.length_squared().mul(rhs.length_squared()))),
+        )
     }
 
     /// Returns some vector that is orthogonal to the given one.
@@ -679,7 +679,7 @@ impl Vec3A {
     #[inline]
     pub fn any_orthogonal_vector(&self) -> Self {
         // This can probably be optimized
-        if self.x.abs() > self.y.abs() {
+        if math::abs(self.x) > math::abs(self.y) {
             Self::new(-self.z, 0.0, self.x) // self.cross(Self::Y)
         } else {
             Self::new(0.0, self.z, -self.y) // self.cross(Self::X)
@@ -696,10 +696,7 @@ impl Vec3A {
     pub fn any_orthonormal_vector(&self) -> Self {
         glam_assert!(self.is_normalized());
         // From https://graphics.pixar.com/library/OrthonormalB/paper.pdf
-        #[cfg(feature = "std")]
-        let sign = (1.0_f32).copysign(self.z);
-        #[cfg(not(feature = "std"))]
-        let sign = self.z.signum();
+        let sign = math::signum(self.z);
         let a = -1.0 / (sign + self.z);
         let b = self.x * self.y * a;
         Self::new(b, sign + self.y * self.y * a, -self.y)
@@ -715,10 +712,7 @@ impl Vec3A {
     pub fn any_orthonormal_pair(&self) -> (Self, Self) {
         glam_assert!(self.is_normalized());
         // From https://graphics.pixar.com/library/OrthonormalB/paper.pdf
-        #[cfg(feature = "std")]
-        let sign = (1.0_f32).copysign(self.z);
-        #[cfg(not(feature = "std"))]
-        let sign = self.z.signum();
+        let sign = math::signum(self.z);
         let a = -1.0 / (sign + self.z);
         let b = self.x * self.y * a;
         (
diff --git a/src/f32/scalar/vec4.rs b/src/f32/scalar/vec4.rs
index f43c3d1..d55746f 100644
--- a/src/f32/scalar/vec4.rs
+++ b/src/f32/scalar/vec4.rs
@@ -1,16 +1,12 @@
 // Generated from vec.rs.tera template. Edit the template, not the generated file.
 
-use crate::{BVec4, Vec2, Vec3, Vec3A};
+use crate::{f32::math, BVec4, Vec2, Vec3, Vec3A};
 
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
 use core::{f32, ops::*};
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 4-dimensional vector.
 #[inline(always)]
 pub const fn vec4(x: f32, y: f32, z: f32, w: f32) -> Vec4 {
@@ -317,10 +313,10 @@ impl Vec4 {
     #[inline]
     pub fn abs(self) -> Self {
         Self {
-            x: self.x.abs(),
-            y: self.y.abs(),
-            z: self.z.abs(),
-            w: self.w.abs(),
+            x: math::abs(self.x),
+            y: math::abs(self.y),
+            z: math::abs(self.z),
+            w: math::abs(self.w),
         }
     }
 
@@ -332,10 +328,10 @@ impl Vec4 {
     #[inline]
     pub fn signum(self) -> Self {
         Self {
-            x: self.x.signum(),
-            y: self.y.signum(),
-            z: self.z.signum(),
-            w: self.w.signum(),
+            x: math::signum(self.x),
+            y: math::signum(self.y),
+            z: math::signum(self.z),
+            w: math::signum(self.w),
         }
     }
 
@@ -343,10 +339,10 @@ impl Vec4 {
     #[inline]
     pub fn copysign(self, rhs: Self) -> Self {
         Self {
-            x: self.x.copysign(rhs.x),
-            y: self.y.copysign(rhs.y),
-            z: self.z.copysign(rhs.z),
-            w: self.w.copysign(rhs.w),
+            x: math::copysign(self.x, rhs.x),
+            y: math::copysign(self.y, rhs.y),
+            z: math::copysign(self.z, rhs.z),
+            w: math::copysign(self.w, rhs.w),
         }
     }
 
@@ -392,7 +388,7 @@ impl Vec4 {
     #[doc(alias = "magnitude")]
     #[inline]
     pub fn length(self) -> f32 {
-        self.dot(self).sqrt()
+        math::sqrt(self.dot(self))
     }
 
     /// Computes the squared length of `self`.
@@ -482,7 +478,7 @@ impl Vec4 {
     #[inline]
     pub fn is_normalized(self) -> bool {
         // TODO: do something with epsilon
-        (self.length_squared() - 1.0).abs() <= 1e-4
+        math::abs(self.length_squared() - 1.0) <= 1e-4
     }
 
     /// Returns the vector projection of `self` onto `rhs`.
@@ -551,10 +547,10 @@ impl Vec4 {
     #[inline]
     pub fn round(self) -> Self {
         Self {
-            x: self.x.round(),
-            y: self.y.round(),
-            z: self.z.round(),
-            w: self.w.round(),
+            x: math::round(self.x),
+            y: math::round(self.y),
+            z: math::round(self.z),
+            w: math::round(self.w),
         }
     }
 
@@ -563,10 +559,10 @@ impl Vec4 {
     #[inline]
     pub fn floor(self) -> Self {
         Self {
-            x: self.x.floor(),
-            y: self.y.floor(),
-            z: self.z.floor(),
-            w: self.w.floor(),
+            x: math::floor(self.x),
+            y: math::floor(self.y),
+            z: math::floor(self.z),
+            w: math::floor(self.w),
         }
     }
 
@@ -575,10 +571,10 @@ impl Vec4 {
     #[inline]
     pub fn ceil(self) -> Self {
         Self {
-            x: self.x.ceil(),
-            y: self.y.ceil(),
-            z: self.z.ceil(),
-            w: self.w.ceil(),
+            x: math::ceil(self.x),
+            y: math::ceil(self.y),
+            z: math::ceil(self.z),
+            w: math::ceil(self.w),
         }
     }
 
@@ -595,17 +591,22 @@ impl Vec4 {
     /// `self`.
     #[inline]
     pub fn exp(self) -> Self {
-        Self::new(self.x.exp(), self.y.exp(), self.z.exp(), self.w.exp())
+        Self::new(
+            math::exp(self.x),
+            math::exp(self.y),
+            math::exp(self.z),
+            math::exp(self.w),
+        )
     }
 
     /// Returns a vector containing each element of `self` raised to the power of `n`.
     #[inline]
     pub fn powf(self, n: f32) -> Self {
         Self::new(
-            self.x.powf(n),
-            self.y.powf(n),
-            self.z.powf(n),
-            self.w.powf(n),
+            math::powf(self.x, n),
+            math::powf(self.y, n),
+            math::powf(self.z, n),
+            math::powf(self.w, n),
         )
     }
 
@@ -613,10 +614,10 @@ impl Vec4 {
     #[inline]
     pub fn recip(self) -> Self {
         Self {
-            x: self.x.recip(),
-            y: self.y.recip(),
-            z: self.z.recip(),
-            w: self.w.recip(),
+            x: 1.0 / self.x,
+            y: 1.0 / self.y,
+            z: 1.0 / self.z,
+            w: 1.0 / self.w,
         }
     }
 
@@ -655,9 +656,9 @@ impl Vec4 {
         glam_assert!(min <= max);
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -667,7 +668,7 @@ impl Vec4 {
     pub fn clamp_length_max(self, max: f32) -> Self {
         let length_sq = self.length_squared();
         if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -677,7 +678,7 @@ impl Vec4 {
     pub fn clamp_length_min(self, min: f32) -> Self {
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -693,10 +694,10 @@ impl Vec4 {
     #[inline]
     pub fn mul_add(self, a: Self, b: Self) -> Self {
         Self::new(
-            self.x.mul_add(a.x, b.x),
-            self.y.mul_add(a.y, b.y),
-            self.z.mul_add(a.z, b.z),
-            self.w.mul_add(a.w, b.w),
+            math::mul_add(self.x, a.x, b.x),
+            math::mul_add(self.y, a.y, b.y),
+            math::mul_add(self.z, a.z, b.z),
+            math::mul_add(self.w, a.w, b.w),
         )
     }
 
diff --git a/src/f32/sse2/mat2.rs b/src/f32/sse2/mat2.rs
index 556c0cc..4b4baec 100644
--- a/src/f32/sse2/mat2.rs
+++ b/src/f32/sse2/mat2.rs
@@ -1,6 +1,6 @@
 // Generated from mat.rs.tera template. Edit the template, not the generated file.
 
-use crate::{swizzles::*, DMat2, Mat3, Mat3A, Vec2};
+use crate::{f32::math, swizzles::*, DMat2, Mat3, Mat3A, Vec2};
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
@@ -11,10 +11,6 @@ use core::arch::x86::*;
 #[cfg(target_arch = "x86_64")]
 use core::arch::x86_64::*;
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 union UnionCast {
     a: [f32; 4],
     v: Mat2,
@@ -108,14 +104,14 @@ impl Mat2 {
     /// `angle` (in radians).
     #[inline]
     pub fn from_scale_angle(scale: Vec2, angle: f32) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::new(cos * scale.x, sin * scale.x, -sin * scale.y, cos * scale.y)
     }
 
     /// Creates a 2x2 matrix containing a rotation of `angle` (in radians).
     #[inline]
     pub fn from_angle(angle: f32) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::new(cos, sin, -sin, cos)
     }
 
diff --git a/src/f32/sse2/mat3a.rs b/src/f32/sse2/mat3a.rs
index 14912dd..3ff0e5d 100644
--- a/src/f32/sse2/mat3a.rs
+++ b/src/f32/sse2/mat3a.rs
@@ -1,6 +1,6 @@
 // Generated from mat.rs.tera template. Edit the template, not the generated file.
 
-use crate::{swizzles::*, DMat3, EulerRot, Mat2, Mat3, Mat4, Quat, Vec2, Vec3, Vec3A};
+use crate::{f32::math, swizzles::*, DMat3, EulerRot, Mat2, Mat3, Mat4, Quat, Vec2, Vec3, Vec3A};
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
@@ -11,10 +11,6 @@ use core::arch::x86::*;
 #[cfg(target_arch = "x86_64")]
 use core::arch::x86_64::*;
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 3x3 matrix from column vectors.
 #[inline(always)]
 pub const fn mat3a(x_axis: Vec3A, y_axis: Vec3A, z_axis: Vec3A) -> Mat3A {
@@ -191,7 +187,7 @@ impl Mat3A {
     pub fn from_axis_angle(axis: Vec3, angle: f32) -> Self {
         glam_assert!(axis.is_normalized());
 
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         let (xsin, ysin, zsin) = axis.mul(sin).into();
         let (x, y, z) = axis.into();
         let (x2, y2, z2) = axis.mul(axis).into();
@@ -217,7 +213,7 @@ impl Mat3A {
     /// Creates a 3D rotation matrix from `angle` (in radians) around the x axis.
     #[inline]
     pub fn from_rotation_x(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec3A::X,
             Vec3A::new(0.0, cosa, sina),
@@ -228,7 +224,7 @@ impl Mat3A {
     /// Creates a 3D rotation matrix from `angle` (in radians) around the y axis.
     #[inline]
     pub fn from_rotation_y(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec3A::new(cosa, 0.0, -sina),
             Vec3A::Y,
@@ -239,7 +235,7 @@ impl Mat3A {
     /// Creates a 3D rotation matrix from `angle` (in radians) around the z axis.
     #[inline]
     pub fn from_rotation_z(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec3A::new(cosa, sina, 0.0),
             Vec3A::new(-sina, cosa, 0.0),
@@ -267,7 +263,7 @@ impl Mat3A {
     /// [`Self::transform_point2()`] and [`Self::transform_vector2()`].
     #[inline]
     pub fn from_angle(angle: f32) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::from_cols(
             Vec3A::new(cos, sin, 0.0),
             Vec3A::new(-sin, cos, 0.0),
@@ -282,7 +278,7 @@ impl Mat3A {
     /// [`Self::transform_point2()`] and [`Self::transform_vector2()`].
     #[inline]
     pub fn from_scale_angle_translation(scale: Vec2, angle: f32, translation: Vec2) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::from_cols(
             Vec3A::new(cos * scale.x, sin * scale.x, 0.0),
             Vec3A::new(-sin * scale.y, cos * scale.y, 0.0),
diff --git a/src/f32/sse2/mat4.rs b/src/f32/sse2/mat4.rs
index aa69493..82e7635 100644
--- a/src/f32/sse2/mat4.rs
+++ b/src/f32/sse2/mat4.rs
@@ -1,6 +1,8 @@
 // Generated from mat.rs.tera template. Edit the template, not the generated file.
 
-use crate::{sse2::*, swizzles::*, DMat4, EulerRot, Mat3, Mat3A, Quat, Vec3, Vec3A, Vec4};
+use crate::{
+    f32::math, sse2::*, swizzles::*, DMat4, EulerRot, Mat3, Mat3A, Quat, Vec3, Vec3A, Vec4,
+};
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
@@ -11,10 +13,6 @@ use core::arch::x86::*;
 #[cfg(target_arch = "x86_64")]
 use core::arch::x86_64::*;
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 4x4 matrix from column vectors.
 #[inline(always)]
 pub const fn mat4(x_axis: Vec4, y_axis: Vec4, z_axis: Vec4, w_axis: Vec4) -> Mat4 {
@@ -241,7 +239,7 @@ impl Mat4 {
         glam_assert!(det != 0.0);
 
         let scale = Vec3::new(
-            self.x_axis.length() * det.signum(),
+            self.x_axis.length() * math::signum(det),
             self.y_axis.length(),
             self.z_axis.length(),
         );
@@ -332,7 +330,7 @@ impl Mat4 {
     pub fn from_axis_angle(axis: Vec3, angle: f32) -> Self {
         glam_assert!(axis.is_normalized());
 
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         let axis_sin = axis.mul(sin);
         let axis_sq = axis.mul(axis);
         let omc = 1.0 - cos;
@@ -380,7 +378,7 @@ impl Mat4 {
     /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
     #[inline]
     pub fn from_rotation_x(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec4::X,
             Vec4::new(0.0, cosa, sina, 0.0),
@@ -396,7 +394,7 @@ impl Mat4 {
     /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
     #[inline]
     pub fn from_rotation_y(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec4::new(cosa, 0.0, -sina, 0.0),
             Vec4::Y,
@@ -412,7 +410,7 @@ impl Mat4 {
     /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
     #[inline]
     pub fn from_rotation_z(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec4::new(cosa, sina, 0.0, 0.0),
             Vec4::new(-sina, cosa, 0.0, 0.0),
@@ -817,7 +815,7 @@ impl Mat4 {
         z_far: f32,
     ) -> Self {
         let inv_length = 1.0 / (z_near - z_far);
-        let f = 1.0 / (0.5 * fov_y_radians).tan();
+        let f = 1.0 / math::tan(0.5 * fov_y_radians);
         let a = f / aspect_ratio;
         let b = (z_near + z_far) * inv_length;
         let c = (2.0 * z_near * z_far) * inv_length;
@@ -838,7 +836,7 @@ impl Mat4 {
     #[inline]
     pub fn perspective_lh(fov_y_radians: f32, aspect_ratio: f32, z_near: f32, z_far: f32) -> Self {
         glam_assert!(z_near > 0.0 && z_far > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         let r = z_far / (z_far - z_near);
@@ -859,7 +857,7 @@ impl Mat4 {
     #[inline]
     pub fn perspective_rh(fov_y_radians: f32, aspect_ratio: f32, z_near: f32, z_far: f32) -> Self {
         glam_assert!(z_near > 0.0 && z_far > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         let r = z_far / (z_near - z_far);
@@ -879,7 +877,7 @@ impl Mat4 {
     #[inline]
     pub fn perspective_infinite_lh(fov_y_radians: f32, aspect_ratio: f32, z_near: f32) -> Self {
         glam_assert!(z_near > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         Self::from_cols(
@@ -902,7 +900,7 @@ impl Mat4 {
         z_near: f32,
     ) -> Self {
         glam_assert!(z_near > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         Self::from_cols(
@@ -918,7 +916,7 @@ impl Mat4 {
     #[inline]
     pub fn perspective_infinite_rh(fov_y_radians: f32, aspect_ratio: f32, z_near: f32) -> Self {
         glam_assert!(z_near > 0.0);
-        let f = 1.0 / (0.5 * fov_y_radians).tan();
+        let f = 1.0 / math::tan(0.5 * fov_y_radians);
         Self::from_cols(
             Vec4::new(f / aspect_ratio, 0.0, 0.0, 0.0),
             Vec4::new(0.0, f, 0.0, 0.0),
@@ -936,7 +934,7 @@ impl Mat4 {
         z_near: f32,
     ) -> Self {
         glam_assert!(z_near > 0.0);
-        let f = 1.0 / (0.5 * fov_y_radians).tan();
+        let f = 1.0 / math::tan(0.5 * fov_y_radians);
         Self::from_cols(
             Vec4::new(f / aspect_ratio, 0.0, 0.0, 0.0),
             Vec4::new(0.0, f, 0.0, 0.0),
diff --git a/src/f32/sse2/quat.rs b/src/f32/sse2/quat.rs
index 12fdf84..0ba0e91 100644
--- a/src/f32/sse2/quat.rs
+++ b/src/f32/sse2/quat.rs
@@ -2,14 +2,11 @@
 
 use crate::{
     euler::{EulerFromQuaternion, EulerRot, EulerToQuaternion},
+    f32::math,
     sse2::*,
-    DQuat, FloatEx, Mat3, Mat3A, Mat4, Vec2, Vec3, Vec3A, Vec4,
+    DQuat, Mat3, Mat3A, Mat4, Vec2, Vec3, Vec3A, Vec4,
 };
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 #[cfg(target_arch = "x86")]
 use core::arch::x86::*;
 #[cfg(target_arch = "x86_64")]
@@ -131,7 +128,7 @@ impl Quat {
     #[inline]
     pub fn from_axis_angle(axis: Vec3, angle: f32) -> Self {
         glam_assert!(axis.is_normalized());
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         let v = axis * s;
         Self::from_xyzw(v.x, v.y, v.z, c)
     }
@@ -152,21 +149,21 @@ impl Quat {
     /// Creates a quaternion from the `angle` (in radians) around the x axis.
     #[inline]
     pub fn from_rotation_x(angle: f32) -> Self {
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         Self::from_xyzw(s, 0.0, 0.0, c)
     }
 
     /// Creates a quaternion from the `angle` (in radians) around the y axis.
     #[inline]
     pub fn from_rotation_y(angle: f32) -> Self {
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         Self::from_xyzw(0.0, s, 0.0, c)
     }
 
     /// Creates a quaternion from the `angle` (in radians) around the z axis.
     #[inline]
     pub fn from_rotation_z(angle: f32) -> Self {
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         Self::from_xyzw(0.0, 0.0, s, c)
     }
 
@@ -190,7 +187,7 @@ impl Quat {
             if dif10 <= 0.0 {
                 // x^2 >= y^2
                 let four_xsq = omm22 - dif10;
-                let inv4x = 0.5 / four_xsq.sqrt();
+                let inv4x = 0.5 / math::sqrt(four_xsq);
                 Self::from_xyzw(
                     four_xsq * inv4x,
                     (m01 + m10) * inv4x,
@@ -200,7 +197,7 @@ impl Quat {
             } else {
                 // y^2 >= x^2
                 let four_ysq = omm22 + dif10;
-                let inv4y = 0.5 / four_ysq.sqrt();
+                let inv4y = 0.5 / math::sqrt(four_ysq);
                 Self::from_xyzw(
                     (m01 + m10) * inv4y,
                     four_ysq * inv4y,
@@ -215,7 +212,7 @@ impl Quat {
             if sum10 <= 0.0 {
                 // z^2 >= w^2
                 let four_zsq = opm22 - sum10;
-                let inv4z = 0.5 / four_zsq.sqrt();
+                let inv4z = 0.5 / math::sqrt(four_zsq);
                 Self::from_xyzw(
                     (m02 + m20) * inv4z,
                     (m12 + m21) * inv4z,
@@ -225,7 +222,7 @@ impl Quat {
             } else {
                 // w^2 >= z^2
                 let four_wsq = opm22 + sum10;
-                let inv4w = 0.5 / four_wsq.sqrt();
+                let inv4w = 0.5 / math::sqrt(four_wsq);
                 Self::from_xyzw(
                     (m12 - m21) * inv4w,
                     (m20 - m02) * inv4w,
@@ -345,7 +342,7 @@ impl Quat {
             let z = from.x * to.y - to.x * from.y;
             let w = 1.0 + dot;
             // calculate length with x=0 and y=0 to normalize
-            let len_rcp = 1.0 / (z * z + w * w).sqrt();
+            let len_rcp = 1.0 / math::sqrt(z * z + w * w);
             Self::from_xyzw(0.0, 0.0, z * len_rcp, w * len_rcp)
         }
     }
@@ -357,7 +354,7 @@ impl Quat {
         let v = Vec3::new(self.x, self.y, self.z);
         let length = v.length();
         if length >= EPSILON {
-            let angle = 2.0 * f32::atan2(length, self.w);
+            let angle = 2.0 * math::atan2(length, self.w);
             let axis = v / length;
             (axis, angle)
         } else {
@@ -497,7 +494,7 @@ impl Quat {
         // If the quat.w is close to -1.0, the angle will be near 2*PI which is close to
         // a negative 0 rotation. By forcing quat.w to be positive, we'll end up with
         // the shortest path.
-        let positive_w_angle = self.w.abs().acos_approx() * 2.0;
+        let positive_w_angle = math::acos_approx(math::abs(self.w)) * 2.0;
         positive_w_angle < threshold_angle
     }
 
@@ -512,7 +509,7 @@ impl Quat {
     #[inline]
     pub fn angle_between(self, rhs: Self) -> f32 {
         glam_assert!(self.is_normalized() && rhs.is_normalized());
-        self.dot(rhs).abs().acos_approx() * 2.0
+        math::acos_approx(math::abs(self.dot(rhs))) * 2.0
     }
 
     /// Returns true if the absolute difference of all elements between `self` and `rhs`
@@ -593,7 +590,7 @@ impl Quat {
             // assumes lerp returns a normalized quaternion
             self.lerp(end, s)
         } else {
-            let theta = dot.acos_approx();
+            let theta = math::acos_approx(dot);
 
             let x = 1.0 - s;
             let y = s;
diff --git a/src/f32/sse2/vec3a.rs b/src/f32/sse2/vec3a.rs
index 3d6d919..7c6267e 100644
--- a/src/f32/sse2/vec3a.rs
+++ b/src/f32/sse2/vec3a.rs
@@ -1,6 +1,6 @@
 // Generated from vec.rs.tera template. Edit the template, not the generated file.
 
-use crate::{sse2::*, BVec3A, Vec2, Vec3, Vec4};
+use crate::{f32::math, sse2::*, BVec3A, Vec2, Vec3, Vec4};
 
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
@@ -12,10 +12,6 @@ use core::arch::x86::*;
 #[cfg(target_arch = "x86_64")]
 use core::arch::x86_64::*;
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 union UnionCast {
     a: [f32; 4],
     v: Vec3A,
@@ -467,7 +463,7 @@ impl Vec3A {
     #[inline]
     pub fn is_normalized(self) -> bool {
         // TODO: do something with epsilon
-        (self.length_squared() - 1.0).abs() <= 1e-4
+        math::abs(self.length_squared() - 1.0) <= 1e-4
     }
 
     /// Returns the vector projection of `self` onto `rhs`.
@@ -565,13 +561,17 @@ impl Vec3A {
     /// `self`.
     #[inline]
     pub fn exp(self) -> Self {
-        Self::new(self.x.exp(), self.y.exp(), self.z.exp())
+        Self::new(math::exp(self.x), math::exp(self.y), math::exp(self.z))
     }
 
     /// Returns a vector containing each element of `self` raised to the power of `n`.
     #[inline]
     pub fn powf(self, n: f32) -> Self {
-        Self::new(self.x.powf(n), self.y.powf(n), self.z.powf(n))
+        Self::new(
+            math::powf(self.x, n),
+            math::powf(self.y, n),
+            math::powf(self.z, n),
+        )
     }
 
     /// Returns a vector containing the reciprocal `1.0/n` of each element of `self`.
@@ -615,9 +615,9 @@ impl Vec3A {
         glam_assert!(min <= max);
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -627,7 +627,7 @@ impl Vec3A {
     pub fn clamp_length_max(self, max: f32) -> Self {
         let length_sq = self.length_squared();
         if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -637,7 +637,7 @@ impl Vec3A {
     pub fn clamp_length_min(self, min: f32) -> Self {
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -658,9 +658,9 @@ impl Vec3A {
         }
         #[cfg(not(target_feature = "fma"))]
         Self::new(
-            self.x.mul_add(a.x, b.x),
-            self.y.mul_add(a.y, b.y),
-            self.z.mul_add(a.z, b.z),
+            math::mul_add(self.x, a.x, b.x),
+            math::mul_add(self.y, a.y, b.y),
+            math::mul_add(self.z, a.z, b.z),
         )
     }
 
@@ -669,10 +669,10 @@ impl Vec3A {
     /// The input vectors do not need to be unit length however they must be non-zero.
     #[inline]
     pub fn angle_between(self, rhs: Self) -> f32 {
-        use crate::FloatEx;
-        self.dot(rhs)
-            .div(self.length_squared().mul(rhs.length_squared()).sqrt())
-            .acos_approx()
+        math::acos_approx(
+            self.dot(rhs)
+                .div(math::sqrt(self.length_squared().mul(rhs.length_squared()))),
+        )
     }
 
     /// Returns some vector that is orthogonal to the given one.
@@ -684,7 +684,7 @@ impl Vec3A {
     #[inline]
     pub fn any_orthogonal_vector(&self) -> Self {
         // This can probably be optimized
-        if self.x.abs() > self.y.abs() {
+        if math::abs(self.x) > math::abs(self.y) {
             Self::new(-self.z, 0.0, self.x) // self.cross(Self::Y)
         } else {
             Self::new(0.0, self.z, -self.y) // self.cross(Self::X)
@@ -701,10 +701,7 @@ impl Vec3A {
     pub fn any_orthonormal_vector(&self) -> Self {
         glam_assert!(self.is_normalized());
         // From https://graphics.pixar.com/library/OrthonormalB/paper.pdf
-        #[cfg(feature = "std")]
-        let sign = (1.0_f32).copysign(self.z);
-        #[cfg(not(feature = "std"))]
-        let sign = self.z.signum();
+        let sign = math::signum(self.z);
         let a = -1.0 / (sign + self.z);
         let b = self.x * self.y * a;
         Self::new(b, sign + self.y * self.y * a, -self.y)
@@ -720,10 +717,7 @@ impl Vec3A {
     pub fn any_orthonormal_pair(&self) -> (Self, Self) {
         glam_assert!(self.is_normalized());
         // From https://graphics.pixar.com/library/OrthonormalB/paper.pdf
-        #[cfg(feature = "std")]
-        let sign = (1.0_f32).copysign(self.z);
-        #[cfg(not(feature = "std"))]
-        let sign = self.z.signum();
+        let sign = math::signum(self.z);
         let a = -1.0 / (sign + self.z);
         let b = self.x * self.y * a;
         (
diff --git a/src/f32/sse2/vec4.rs b/src/f32/sse2/vec4.rs
index d492c7a..d896dbd 100644
--- a/src/f32/sse2/vec4.rs
+++ b/src/f32/sse2/vec4.rs
@@ -1,6 +1,6 @@
 // Generated from vec.rs.tera template. Edit the template, not the generated file.
 
-use crate::{sse2::*, BVec4A, Vec2, Vec3, Vec3A};
+use crate::{f32::math, sse2::*, BVec4A, Vec2, Vec3, Vec3A};
 
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
@@ -12,10 +12,6 @@ use core::arch::x86::*;
 #[cfg(target_arch = "x86_64")]
 use core::arch::x86_64::*;
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 union UnionCast {
     a: [f32; 4],
     v: Vec4,
@@ -441,7 +437,7 @@ impl Vec4 {
     #[inline]
     pub fn is_normalized(self) -> bool {
         // TODO: do something with epsilon
-        (self.length_squared() - 1.0).abs() <= 1e-4
+        math::abs(self.length_squared() - 1.0) <= 1e-4
     }
 
     /// Returns the vector projection of `self` onto `rhs`.
@@ -539,17 +535,22 @@ impl Vec4 {
     /// `self`.
     #[inline]
     pub fn exp(self) -> Self {
-        Self::new(self.x.exp(), self.y.exp(), self.z.exp(), self.w.exp())
+        Self::new(
+            math::exp(self.x),
+            math::exp(self.y),
+            math::exp(self.z),
+            math::exp(self.w),
+        )
     }
 
     /// Returns a vector containing each element of `self` raised to the power of `n`.
     #[inline]
     pub fn powf(self, n: f32) -> Self {
         Self::new(
-            self.x.powf(n),
-            self.y.powf(n),
-            self.z.powf(n),
-            self.w.powf(n),
+            math::powf(self.x, n),
+            math::powf(self.y, n),
+            math::powf(self.z, n),
+            math::powf(self.w, n),
         )
     }
 
@@ -594,9 +595,9 @@ impl Vec4 {
         glam_assert!(min <= max);
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -606,7 +607,7 @@ impl Vec4 {
     pub fn clamp_length_max(self, max: f32) -> Self {
         let length_sq = self.length_squared();
         if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -616,7 +617,7 @@ impl Vec4 {
     pub fn clamp_length_min(self, min: f32) -> Self {
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -637,10 +638,10 @@ impl Vec4 {
         }
         #[cfg(not(target_feature = "fma"))]
         Self::new(
-            self.x.mul_add(a.x, b.x),
-            self.y.mul_add(a.y, b.y),
-            self.z.mul_add(a.z, b.z),
-            self.w.mul_add(a.w, b.w),
+            math::mul_add(self.x, a.x, b.x),
+            math::mul_add(self.y, a.y, b.y),
+            math::mul_add(self.z, a.z, b.z),
+            math::mul_add(self.w, a.w, b.w),
         )
     }
 
diff --git a/src/f32/vec2.rs b/src/f32/vec2.rs
index 5cc7d6d..ef2d77e 100644
--- a/src/f32/vec2.rs
+++ b/src/f32/vec2.rs
@@ -1,16 +1,12 @@
 // Generated from vec.rs.tera template. Edit the template, not the generated file.
 
-use crate::{BVec2, Vec3};
+use crate::{f32::math, BVec2, Vec3};
 
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
 use core::{f32, ops::*};
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 2-dimensional vector.
 #[inline(always)]
 pub const fn vec2(x: f32, y: f32) -> Vec2 {
@@ -246,8 +242,8 @@ impl Vec2 {
     #[inline]
     pub fn abs(self) -> Self {
         Self {
-            x: self.x.abs(),
-            y: self.y.abs(),
+            x: math::abs(self.x),
+            y: math::abs(self.y),
         }
     }
 
@@ -259,8 +255,8 @@ impl Vec2 {
     #[inline]
     pub fn signum(self) -> Self {
         Self {
-            x: self.x.signum(),
-            y: self.y.signum(),
+            x: math::signum(self.x),
+            y: math::signum(self.y),
         }
     }
 
@@ -268,8 +264,8 @@ impl Vec2 {
     #[inline]
     pub fn copysign(self, rhs: Self) -> Self {
         Self {
-            x: self.x.copysign(rhs.x),
-            y: self.y.copysign(rhs.y),
+            x: math::copysign(self.x, rhs.x),
+            y: math::copysign(self.y, rhs.y),
         }
     }
 
@@ -307,7 +303,7 @@ impl Vec2 {
     #[doc(alias = "magnitude")]
     #[inline]
     pub fn length(self) -> f32 {
-        self.dot(self).sqrt()
+        math::sqrt(self.dot(self))
     }
 
     /// Computes the squared length of `self`.
@@ -397,7 +393,7 @@ impl Vec2 {
     #[inline]
     pub fn is_normalized(self) -> bool {
         // TODO: do something with epsilon
-        (self.length_squared() - 1.0).abs() <= 1e-4
+        math::abs(self.length_squared() - 1.0) <= 1e-4
     }
 
     /// Returns the vector projection of `self` onto `rhs`.
@@ -466,8 +462,8 @@ impl Vec2 {
     #[inline]
     pub fn round(self) -> Self {
         Self {
-            x: self.x.round(),
-            y: self.y.round(),
+            x: math::round(self.x),
+            y: math::round(self.y),
         }
     }
 
@@ -476,8 +472,8 @@ impl Vec2 {
     #[inline]
     pub fn floor(self) -> Self {
         Self {
-            x: self.x.floor(),
-            y: self.y.floor(),
+            x: math::floor(self.x),
+            y: math::floor(self.y),
         }
     }
 
@@ -486,8 +482,8 @@ impl Vec2 {
     #[inline]
     pub fn ceil(self) -> Self {
         Self {
-            x: self.x.ceil(),
-            y: self.y.ceil(),
+            x: math::ceil(self.x),
+            y: math::ceil(self.y),
         }
     }
 
@@ -504,21 +500,21 @@ impl Vec2 {
     /// `self`.
     #[inline]
     pub fn exp(self) -> Self {
-        Self::new(self.x.exp(), self.y.exp())
+        Self::new(math::exp(self.x), math::exp(self.y))
     }
 
     /// Returns a vector containing each element of `self` raised to the power of `n`.
     #[inline]
     pub fn powf(self, n: f32) -> Self {
-        Self::new(self.x.powf(n), self.y.powf(n))
+        Self::new(math::powf(self.x, n), math::powf(self.y, n))
     }
 
     /// Returns a vector containing the reciprocal `1.0/n` of each element of `self`.
     #[inline]
     pub fn recip(self) -> Self {
         Self {
-            x: self.x.recip(),
-            y: self.y.recip(),
+            x: 1.0 / self.x,
+            y: 1.0 / self.y,
         }
     }
 
@@ -557,9 +553,9 @@ impl Vec2 {
         glam_assert!(min <= max);
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -569,7 +565,7 @@ impl Vec2 {
     pub fn clamp_length_max(self, max: f32) -> Self {
         let length_sq = self.length_squared();
         if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -579,7 +575,7 @@ impl Vec2 {
     pub fn clamp_length_min(self, min: f32) -> Self {
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -594,7 +590,10 @@ impl Vec2 {
     /// mind.
     #[inline]
     pub fn mul_add(self, a: Self, b: Self) -> Self {
-        Self::new(self.x.mul_add(a.x, b.x), self.y.mul_add(a.y, b.y))
+        Self::new(
+            math::mul_add(self.x, a.x, b.x),
+            math::mul_add(self.y, a.y, b.y),
+        )
     }
 
     /// Creates a 2D vector containing `[angle.cos(), angle.sin()]`. This can be used in
@@ -602,7 +601,7 @@ impl Vec2 {
     /// create the vector [-1, 0] and rotate `Vec2::Y` around it returning `-Vec2::Y`.
     #[inline]
     pub fn from_angle(angle: f32) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self { x: cos, y: sin }
     }
 
@@ -611,11 +610,11 @@ impl Vec2 {
     /// The input vectors do not need to be unit length however they must be non-zero.
     #[inline]
     pub fn angle_between(self, rhs: Self) -> f32 {
-        use crate::FloatEx;
-        let angle =
-            (self.dot(rhs) / (self.length_squared() * rhs.length_squared()).sqrt()).acos_approx();
+        let angle = math::acos_approx(
+            self.dot(rhs) / math::sqrt(self.length_squared() * rhs.length_squared()),
+        );
 
-        angle * self.perp_dot(rhs).signum()
+        angle * math::signum(self.perp_dot(rhs))
     }
 
     /// Returns a vector that is equal to `self` rotated by 90 degrees.
diff --git a/src/f32/vec3.rs b/src/f32/vec3.rs
index ae8e818..243a707 100644
--- a/src/f32/vec3.rs
+++ b/src/f32/vec3.rs
@@ -1,16 +1,12 @@
 // Generated from vec.rs.tera template. Edit the template, not the generated file.
 
-use crate::{BVec3, Vec2, Vec4};
+use crate::{f32::math, BVec3, Vec2, Vec4};
 
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
 use core::{f32, ops::*};
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 3-dimensional vector.
 #[inline(always)]
 pub const fn vec3(x: f32, y: f32, z: f32) -> Vec3 {
@@ -286,9 +282,9 @@ impl Vec3 {
     #[inline]
     pub fn abs(self) -> Self {
         Self {
-            x: self.x.abs(),
-            y: self.y.abs(),
-            z: self.z.abs(),
+            x: math::abs(self.x),
+            y: math::abs(self.y),
+            z: math::abs(self.z),
         }
     }
 
@@ -300,9 +296,9 @@ impl Vec3 {
     #[inline]
     pub fn signum(self) -> Self {
         Self {
-            x: self.x.signum(),
-            y: self.y.signum(),
-            z: self.z.signum(),
+            x: math::signum(self.x),
+            y: math::signum(self.y),
+            z: math::signum(self.z),
         }
     }
 
@@ -310,9 +306,9 @@ impl Vec3 {
     #[inline]
     pub fn copysign(self, rhs: Self) -> Self {
         Self {
-            x: self.x.copysign(rhs.x),
-            y: self.y.copysign(rhs.y),
-            z: self.z.copysign(rhs.z),
+            x: math::copysign(self.x, rhs.x),
+            y: math::copysign(self.y, rhs.y),
+            z: math::copysign(self.z, rhs.z),
         }
     }
 
@@ -352,7 +348,7 @@ impl Vec3 {
     #[doc(alias = "magnitude")]
     #[inline]
     pub fn length(self) -> f32 {
-        self.dot(self).sqrt()
+        math::sqrt(self.dot(self))
     }
 
     /// Computes the squared length of `self`.
@@ -442,7 +438,7 @@ impl Vec3 {
     #[inline]
     pub fn is_normalized(self) -> bool {
         // TODO: do something with epsilon
-        (self.length_squared() - 1.0).abs() <= 1e-4
+        math::abs(self.length_squared() - 1.0) <= 1e-4
     }
 
     /// Returns the vector projection of `self` onto `rhs`.
@@ -511,9 +507,9 @@ impl Vec3 {
     #[inline]
     pub fn round(self) -> Self {
         Self {
-            x: self.x.round(),
-            y: self.y.round(),
-            z: self.z.round(),
+            x: math::round(self.x),
+            y: math::round(self.y),
+            z: math::round(self.z),
         }
     }
 
@@ -522,9 +518,9 @@ impl Vec3 {
     #[inline]
     pub fn floor(self) -> Self {
         Self {
-            x: self.x.floor(),
-            y: self.y.floor(),
-            z: self.z.floor(),
+            x: math::floor(self.x),
+            y: math::floor(self.y),
+            z: math::floor(self.z),
         }
     }
 
@@ -533,9 +529,9 @@ impl Vec3 {
     #[inline]
     pub fn ceil(self) -> Self {
         Self {
-            x: self.x.ceil(),
-            y: self.y.ceil(),
-            z: self.z.ceil(),
+            x: math::ceil(self.x),
+            y: math::ceil(self.y),
+            z: math::ceil(self.z),
         }
     }
 
@@ -552,22 +548,26 @@ impl Vec3 {
     /// `self`.
     #[inline]
     pub fn exp(self) -> Self {
-        Self::new(self.x.exp(), self.y.exp(), self.z.exp())
+        Self::new(math::exp(self.x), math::exp(self.y), math::exp(self.z))
     }
 
     /// Returns a vector containing each element of `self` raised to the power of `n`.
     #[inline]
     pub fn powf(self, n: f32) -> Self {
-        Self::new(self.x.powf(n), self.y.powf(n), self.z.powf(n))
+        Self::new(
+            math::powf(self.x, n),
+            math::powf(self.y, n),
+            math::powf(self.z, n),
+        )
     }
 
     /// Returns a vector containing the reciprocal `1.0/n` of each element of `self`.
     #[inline]
     pub fn recip(self) -> Self {
         Self {
-            x: self.x.recip(),
-            y: self.y.recip(),
-            z: self.z.recip(),
+            x: 1.0 / self.x,
+            y: 1.0 / self.y,
+            z: 1.0 / self.z,
         }
     }
 
@@ -606,9 +606,9 @@ impl Vec3 {
         glam_assert!(min <= max);
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -618,7 +618,7 @@ impl Vec3 {
     pub fn clamp_length_max(self, max: f32) -> Self {
         let length_sq = self.length_squared();
         if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -628,7 +628,7 @@ impl Vec3 {
     pub fn clamp_length_min(self, min: f32) -> Self {
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -644,9 +644,9 @@ impl Vec3 {
     #[inline]
     pub fn mul_add(self, a: Self, b: Self) -> Self {
         Self::new(
-            self.x.mul_add(a.x, b.x),
-            self.y.mul_add(a.y, b.y),
-            self.z.mul_add(a.z, b.z),
+            math::mul_add(self.x, a.x, b.x),
+            math::mul_add(self.y, a.y, b.y),
+            math::mul_add(self.z, a.z, b.z),
         )
     }
 
@@ -655,10 +655,10 @@ impl Vec3 {
     /// The input vectors do not need to be unit length however they must be non-zero.
     #[inline]
     pub fn angle_between(self, rhs: Self) -> f32 {
-        use crate::FloatEx;
-        self.dot(rhs)
-            .div(self.length_squared().mul(rhs.length_squared()).sqrt())
-            .acos_approx()
+        math::acos_approx(
+            self.dot(rhs)
+                .div(math::sqrt(self.length_squared().mul(rhs.length_squared()))),
+        )
     }
 
     /// Returns some vector that is orthogonal to the given one.
@@ -670,7 +670,7 @@ impl Vec3 {
     #[inline]
     pub fn any_orthogonal_vector(&self) -> Self {
         // This can probably be optimized
-        if self.x.abs() > self.y.abs() {
+        if math::abs(self.x) > math::abs(self.y) {
             Self::new(-self.z, 0.0, self.x) // self.cross(Self::Y)
         } else {
             Self::new(0.0, self.z, -self.y) // self.cross(Self::X)
@@ -687,10 +687,7 @@ impl Vec3 {
     pub fn any_orthonormal_vector(&self) -> Self {
         glam_assert!(self.is_normalized());
         // From https://graphics.pixar.com/library/OrthonormalB/paper.pdf
-        #[cfg(feature = "std")]
-        let sign = (1.0_f32).copysign(self.z);
-        #[cfg(not(feature = "std"))]
-        let sign = self.z.signum();
+        let sign = math::signum(self.z);
         let a = -1.0 / (sign + self.z);
         let b = self.x * self.y * a;
         Self::new(b, sign + self.y * self.y * a, -self.y)
@@ -706,10 +703,7 @@ impl Vec3 {
     pub fn any_orthonormal_pair(&self) -> (Self, Self) {
         glam_assert!(self.is_normalized());
         // From https://graphics.pixar.com/library/OrthonormalB/paper.pdf
-        #[cfg(feature = "std")]
-        let sign = (1.0_f32).copysign(self.z);
-        #[cfg(not(feature = "std"))]
-        let sign = self.z.signum();
+        let sign = math::signum(self.z);
         let a = -1.0 / (sign + self.z);
         let b = self.x * self.y * a;
         (
diff --git a/src/f32/wasm32/mat2.rs b/src/f32/wasm32/mat2.rs
index 7e0f184..3a4829d 100644
--- a/src/f32/wasm32/mat2.rs
+++ b/src/f32/wasm32/mat2.rs
@@ -1,6 +1,6 @@
 // Generated from mat.rs.tera template. Edit the template, not the generated file.
 
-use crate::{swizzles::*, DMat2, Mat3, Mat3A, Vec2};
+use crate::{f32::math, swizzles::*, DMat2, Mat3, Mat3A, Vec2};
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
@@ -8,10 +8,6 @@ use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign};
 
 use core::arch::wasm32::*;
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 2x2 matrix from column vectors.
 #[inline(always)]
 pub const fn mat2(x_axis: Vec2, y_axis: Vec2) -> Mat2 {
@@ -90,14 +86,14 @@ impl Mat2 {
     /// `angle` (in radians).
     #[inline]
     pub fn from_scale_angle(scale: Vec2, angle: f32) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::new(cos * scale.x, sin * scale.x, -sin * scale.y, cos * scale.y)
     }
 
     /// Creates a 2x2 matrix containing a rotation of `angle` (in radians).
     #[inline]
     pub fn from_angle(angle: f32) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::new(cos, sin, -sin, cos)
     }
 
diff --git a/src/f32/wasm32/mat3a.rs b/src/f32/wasm32/mat3a.rs
index 573ea1f..98b1a01 100644
--- a/src/f32/wasm32/mat3a.rs
+++ b/src/f32/wasm32/mat3a.rs
@@ -1,6 +1,6 @@
 // Generated from mat.rs.tera template. Edit the template, not the generated file.
 
-use crate::{swizzles::*, DMat3, EulerRot, Mat2, Mat3, Mat4, Quat, Vec2, Vec3, Vec3A};
+use crate::{f32::math, swizzles::*, DMat3, EulerRot, Mat2, Mat3, Mat4, Quat, Vec2, Vec3, Vec3A};
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
@@ -8,10 +8,6 @@ use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign};
 
 use core::arch::wasm32::*;
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 3x3 matrix from column vectors.
 #[inline(always)]
 pub const fn mat3a(x_axis: Vec3A, y_axis: Vec3A, z_axis: Vec3A) -> Mat3A {
@@ -188,7 +184,7 @@ impl Mat3A {
     pub fn from_axis_angle(axis: Vec3, angle: f32) -> Self {
         glam_assert!(axis.is_normalized());
 
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         let (xsin, ysin, zsin) = axis.mul(sin).into();
         let (x, y, z) = axis.into();
         let (x2, y2, z2) = axis.mul(axis).into();
@@ -214,7 +210,7 @@ impl Mat3A {
     /// Creates a 3D rotation matrix from `angle` (in radians) around the x axis.
     #[inline]
     pub fn from_rotation_x(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec3A::X,
             Vec3A::new(0.0, cosa, sina),
@@ -225,7 +221,7 @@ impl Mat3A {
     /// Creates a 3D rotation matrix from `angle` (in radians) around the y axis.
     #[inline]
     pub fn from_rotation_y(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec3A::new(cosa, 0.0, -sina),
             Vec3A::Y,
@@ -236,7 +232,7 @@ impl Mat3A {
     /// Creates a 3D rotation matrix from `angle` (in radians) around the z axis.
     #[inline]
     pub fn from_rotation_z(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec3A::new(cosa, sina, 0.0),
             Vec3A::new(-sina, cosa, 0.0),
@@ -264,7 +260,7 @@ impl Mat3A {
     /// [`Self::transform_point2()`] and [`Self::transform_vector2()`].
     #[inline]
     pub fn from_angle(angle: f32) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::from_cols(
             Vec3A::new(cos, sin, 0.0),
             Vec3A::new(-sin, cos, 0.0),
@@ -279,7 +275,7 @@ impl Mat3A {
     /// [`Self::transform_point2()`] and [`Self::transform_vector2()`].
     #[inline]
     pub fn from_scale_angle_translation(scale: Vec2, angle: f32, translation: Vec2) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::from_cols(
             Vec3A::new(cos * scale.x, sin * scale.x, 0.0),
             Vec3A::new(-sin * scale.y, cos * scale.y, 0.0),
diff --git a/src/f32/wasm32/mat4.rs b/src/f32/wasm32/mat4.rs
index bdaeb29..36c6ac6 100644
--- a/src/f32/wasm32/mat4.rs
+++ b/src/f32/wasm32/mat4.rs
@@ -1,6 +1,8 @@
 // Generated from mat.rs.tera template. Edit the template, not the generated file.
 
-use crate::{swizzles::*, wasm32::*, DMat4, EulerRot, Mat3, Mat3A, Quat, Vec3, Vec3A, Vec4};
+use crate::{
+    f32::math, swizzles::*, wasm32::*, DMat4, EulerRot, Mat3, Mat3A, Quat, Vec3, Vec3A, Vec4,
+};
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
@@ -8,10 +10,6 @@ use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign};
 
 use core::arch::wasm32::*;
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 4x4 matrix from column vectors.
 #[inline(always)]
 pub const fn mat4(x_axis: Vec4, y_axis: Vec4, z_axis: Vec4, w_axis: Vec4) -> Mat4 {
@@ -238,7 +236,7 @@ impl Mat4 {
         glam_assert!(det != 0.0);
 
         let scale = Vec3::new(
-            self.x_axis.length() * det.signum(),
+            self.x_axis.length() * math::signum(det),
             self.y_axis.length(),
             self.z_axis.length(),
         );
@@ -329,7 +327,7 @@ impl Mat4 {
     pub fn from_axis_angle(axis: Vec3, angle: f32) -> Self {
         glam_assert!(axis.is_normalized());
 
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         let axis_sin = axis.mul(sin);
         let axis_sq = axis.mul(axis);
         let omc = 1.0 - cos;
@@ -377,7 +375,7 @@ impl Mat4 {
     /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
     #[inline]
     pub fn from_rotation_x(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec4::X,
             Vec4::new(0.0, cosa, sina, 0.0),
@@ -393,7 +391,7 @@ impl Mat4 {
     /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
     #[inline]
     pub fn from_rotation_y(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec4::new(cosa, 0.0, -sina, 0.0),
             Vec4::Y,
@@ -409,7 +407,7 @@ impl Mat4 {
     /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
     #[inline]
     pub fn from_rotation_z(angle: f32) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             Vec4::new(cosa, sina, 0.0, 0.0),
             Vec4::new(-sina, cosa, 0.0, 0.0),
@@ -808,7 +806,7 @@ impl Mat4 {
         z_far: f32,
     ) -> Self {
         let inv_length = 1.0 / (z_near - z_far);
-        let f = 1.0 / (0.5 * fov_y_radians).tan();
+        let f = 1.0 / math::tan(0.5 * fov_y_radians);
         let a = f / aspect_ratio;
         let b = (z_near + z_far) * inv_length;
         let c = (2.0 * z_near * z_far) * inv_length;
@@ -829,7 +827,7 @@ impl Mat4 {
     #[inline]
     pub fn perspective_lh(fov_y_radians: f32, aspect_ratio: f32, z_near: f32, z_far: f32) -> Self {
         glam_assert!(z_near > 0.0 && z_far > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         let r = z_far / (z_far - z_near);
@@ -850,7 +848,7 @@ impl Mat4 {
     #[inline]
     pub fn perspective_rh(fov_y_radians: f32, aspect_ratio: f32, z_near: f32, z_far: f32) -> Self {
         glam_assert!(z_near > 0.0 && z_far > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         let r = z_far / (z_near - z_far);
@@ -870,7 +868,7 @@ impl Mat4 {
     #[inline]
     pub fn perspective_infinite_lh(fov_y_radians: f32, aspect_ratio: f32, z_near: f32) -> Self {
         glam_assert!(z_near > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         Self::from_cols(
@@ -893,7 +891,7 @@ impl Mat4 {
         z_near: f32,
     ) -> Self {
         glam_assert!(z_near > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         Self::from_cols(
@@ -909,7 +907,7 @@ impl Mat4 {
     #[inline]
     pub fn perspective_infinite_rh(fov_y_radians: f32, aspect_ratio: f32, z_near: f32) -> Self {
         glam_assert!(z_near > 0.0);
-        let f = 1.0 / (0.5 * fov_y_radians).tan();
+        let f = 1.0 / math::tan(0.5 * fov_y_radians);
         Self::from_cols(
             Vec4::new(f / aspect_ratio, 0.0, 0.0, 0.0),
             Vec4::new(0.0, f, 0.0, 0.0),
@@ -927,7 +925,7 @@ impl Mat4 {
         z_near: f32,
     ) -> Self {
         glam_assert!(z_near > 0.0);
-        let f = 1.0 / (0.5 * fov_y_radians).tan();
+        let f = 1.0 / math::tan(0.5 * fov_y_radians);
         Self::from_cols(
             Vec4::new(f / aspect_ratio, 0.0, 0.0, 0.0),
             Vec4::new(0.0, f, 0.0, 0.0),
diff --git a/src/f32/wasm32/quat.rs b/src/f32/wasm32/quat.rs
index e5eed2e..04a0bca 100644
--- a/src/f32/wasm32/quat.rs
+++ b/src/f32/wasm32/quat.rs
@@ -2,14 +2,11 @@
 
 use crate::{
     euler::{EulerFromQuaternion, EulerRot, EulerToQuaternion},
+    f32::math,
     wasm32::*,
-    DQuat, FloatEx, Mat3, Mat3A, Mat4, Vec2, Vec3, Vec3A, Vec4,
+    DQuat, Mat3, Mat3A, Mat4, Vec2, Vec3, Vec3A, Vec4,
 };
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 use core::arch::wasm32::*;
 
 #[cfg(not(target_arch = "spirv"))]
@@ -124,7 +121,7 @@ impl Quat {
     #[inline]
     pub fn from_axis_angle(axis: Vec3, angle: f32) -> Self {
         glam_assert!(axis.is_normalized());
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         let v = axis * s;
         Self::from_xyzw(v.x, v.y, v.z, c)
     }
@@ -145,21 +142,21 @@ impl Quat {
     /// Creates a quaternion from the `angle` (in radians) around the x axis.
     #[inline]
     pub fn from_rotation_x(angle: f32) -> Self {
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         Self::from_xyzw(s, 0.0, 0.0, c)
     }
 
     /// Creates a quaternion from the `angle` (in radians) around the y axis.
     #[inline]
     pub fn from_rotation_y(angle: f32) -> Self {
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         Self::from_xyzw(0.0, s, 0.0, c)
     }
 
     /// Creates a quaternion from the `angle` (in radians) around the z axis.
     #[inline]
     pub fn from_rotation_z(angle: f32) -> Self {
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         Self::from_xyzw(0.0, 0.0, s, c)
     }
 
@@ -183,7 +180,7 @@ impl Quat {
             if dif10 <= 0.0 {
                 // x^2 >= y^2
                 let four_xsq = omm22 - dif10;
-                let inv4x = 0.5 / four_xsq.sqrt();
+                let inv4x = 0.5 / math::sqrt(four_xsq);
                 Self::from_xyzw(
                     four_xsq * inv4x,
                     (m01 + m10) * inv4x,
@@ -193,7 +190,7 @@ impl Quat {
             } else {
                 // y^2 >= x^2
                 let four_ysq = omm22 + dif10;
-                let inv4y = 0.5 / four_ysq.sqrt();
+                let inv4y = 0.5 / math::sqrt(four_ysq);
                 Self::from_xyzw(
                     (m01 + m10) * inv4y,
                     four_ysq * inv4y,
@@ -208,7 +205,7 @@ impl Quat {
             if sum10 <= 0.0 {
                 // z^2 >= w^2
                 let four_zsq = opm22 - sum10;
-                let inv4z = 0.5 / four_zsq.sqrt();
+                let inv4z = 0.5 / math::sqrt(four_zsq);
                 Self::from_xyzw(
                     (m02 + m20) * inv4z,
                     (m12 + m21) * inv4z,
@@ -218,7 +215,7 @@ impl Quat {
             } else {
                 // w^2 >= z^2
                 let four_wsq = opm22 + sum10;
-                let inv4w = 0.5 / four_wsq.sqrt();
+                let inv4w = 0.5 / math::sqrt(four_wsq);
                 Self::from_xyzw(
                     (m12 - m21) * inv4w,
                     (m20 - m02) * inv4w,
@@ -338,7 +335,7 @@ impl Quat {
             let z = from.x * to.y - to.x * from.y;
             let w = 1.0 + dot;
             // calculate length with x=0 and y=0 to normalize
-            let len_rcp = 1.0 / (z * z + w * w).sqrt();
+            let len_rcp = 1.0 / math::sqrt(z * z + w * w);
             Self::from_xyzw(0.0, 0.0, z * len_rcp, w * len_rcp)
         }
     }
@@ -350,7 +347,7 @@ impl Quat {
         let v = Vec3::new(self.x, self.y, self.z);
         let length = v.length();
         if length >= EPSILON {
-            let angle = 2.0 * f32::atan2(length, self.w);
+            let angle = 2.0 * math::atan2(length, self.w);
             let axis = v / length;
             (axis, angle)
         } else {
@@ -490,7 +487,7 @@ impl Quat {
         // If the quat.w is close to -1.0, the angle will be near 2*PI which is close to
         // a negative 0 rotation. By forcing quat.w to be positive, we'll end up with
         // the shortest path.
-        let positive_w_angle = self.w.abs().acos_approx() * 2.0;
+        let positive_w_angle = math::acos_approx(math::abs(self.w)) * 2.0;
         positive_w_angle < threshold_angle
     }
 
@@ -505,7 +502,7 @@ impl Quat {
     #[inline]
     pub fn angle_between(self, rhs: Self) -> f32 {
         glam_assert!(self.is_normalized() && rhs.is_normalized());
-        self.dot(rhs).abs().acos_approx() * 2.0
+        math::acos_approx(math::abs(self.dot(rhs))) * 2.0
     }
 
     /// Returns true if the absolute difference of all elements between `self` and `rhs`
@@ -584,7 +581,7 @@ impl Quat {
             // assumes lerp returns a normalized quaternion
             self.lerp(end, s)
         } else {
-            let theta = dot.acos_approx();
+            let theta = math::acos_approx(dot);
 
             // TODO: v128_sin is broken
             // let x = 1.0 - s;
@@ -593,9 +590,9 @@ impl Quat {
             // let w = 0.0;
             // let tmp = f32x4_mul(f32x4_splat(theta), f32x4(x, y, z, w));
             // let tmp = v128_sin(tmp);
-            let x = (theta * (1.0 - s)).sin();
-            let y = (theta * s).sin();
-            let z = theta.sin();
+            let x = math::sin(theta * (1.0 - s));
+            let y = math::sin(theta * s);
+            let z = math::sin(theta);
             let w = 0.0;
             let tmp = f32x4(x, y, z, w);
 
diff --git a/src/f32/wasm32/vec3a.rs b/src/f32/wasm32/vec3a.rs
index 4b4d585..227ce27 100644
--- a/src/f32/wasm32/vec3a.rs
+++ b/src/f32/wasm32/vec3a.rs
@@ -1,6 +1,6 @@
 // Generated from vec.rs.tera template. Edit the template, not the generated file.
 
-use crate::{wasm32::*, BVec3A, Vec2, Vec3, Vec4};
+use crate::{f32::math, wasm32::*, BVec3A, Vec2, Vec3, Vec4};
 
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
@@ -9,10 +9,6 @@ use core::{f32, ops::*};
 
 use core::arch::wasm32::*;
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 union UnionCast {
     a: [f32; 4],
     v: Vec3A,
@@ -442,7 +438,7 @@ impl Vec3A {
     #[inline]
     pub fn is_normalized(self) -> bool {
         // TODO: do something with epsilon
-        (self.length_squared() - 1.0).abs() <= 1e-4
+        math::abs(self.length_squared() - 1.0) <= 1e-4
     }
 
     /// Returns the vector projection of `self` onto `rhs`.
@@ -540,13 +536,17 @@ impl Vec3A {
     /// `self`.
     #[inline]
     pub fn exp(self) -> Self {
-        Self::new(self.x.exp(), self.y.exp(), self.z.exp())
+        Self::new(math::exp(self.x), math::exp(self.y), math::exp(self.z))
     }
 
     /// Returns a vector containing each element of `self` raised to the power of `n`.
     #[inline]
     pub fn powf(self, n: f32) -> Self {
-        Self::new(self.x.powf(n), self.y.powf(n), self.z.powf(n))
+        Self::new(
+            math::powf(self.x, n),
+            math::powf(self.y, n),
+            math::powf(self.z, n),
+        )
     }
 
     /// Returns a vector containing the reciprocal `1.0/n` of each element of `self`.
@@ -590,9 +590,9 @@ impl Vec3A {
         glam_assert!(min <= max);
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -602,7 +602,7 @@ impl Vec3A {
     pub fn clamp_length_max(self, max: f32) -> Self {
         let length_sq = self.length_squared();
         if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -612,7 +612,7 @@ impl Vec3A {
     pub fn clamp_length_min(self, min: f32) -> Self {
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -628,9 +628,9 @@ impl Vec3A {
     #[inline]
     pub fn mul_add(self, a: Self, b: Self) -> Self {
         Self::new(
-            self.x.mul_add(a.x, b.x),
-            self.y.mul_add(a.y, b.y),
-            self.z.mul_add(a.z, b.z),
+            math::mul_add(self.x, a.x, b.x),
+            math::mul_add(self.y, a.y, b.y),
+            math::mul_add(self.z, a.z, b.z),
         )
     }
 
@@ -639,10 +639,10 @@ impl Vec3A {
     /// The input vectors do not need to be unit length however they must be non-zero.
     #[inline]
     pub fn angle_between(self, rhs: Self) -> f32 {
-        use crate::FloatEx;
-        self.dot(rhs)
-            .div(self.length_squared().mul(rhs.length_squared()).sqrt())
-            .acos_approx()
+        math::acos_approx(
+            self.dot(rhs)
+                .div(math::sqrt(self.length_squared().mul(rhs.length_squared()))),
+        )
     }
 
     /// Returns some vector that is orthogonal to the given one.
@@ -654,7 +654,7 @@ impl Vec3A {
     #[inline]
     pub fn any_orthogonal_vector(&self) -> Self {
         // This can probably be optimized
-        if self.x.abs() > self.y.abs() {
+        if math::abs(self.x) > math::abs(self.y) {
             Self::new(-self.z, 0.0, self.x) // self.cross(Self::Y)
         } else {
             Self::new(0.0, self.z, -self.y) // self.cross(Self::X)
@@ -671,10 +671,7 @@ impl Vec3A {
     pub fn any_orthonormal_vector(&self) -> Self {
         glam_assert!(self.is_normalized());
         // From https://graphics.pixar.com/library/OrthonormalB/paper.pdf
-        #[cfg(feature = "std")]
-        let sign = (1.0_f32).copysign(self.z);
-        #[cfg(not(feature = "std"))]
-        let sign = self.z.signum();
+        let sign = math::signum(self.z);
         let a = -1.0 / (sign + self.z);
         let b = self.x * self.y * a;
         Self::new(b, sign + self.y * self.y * a, -self.y)
@@ -690,10 +687,7 @@ impl Vec3A {
     pub fn any_orthonormal_pair(&self) -> (Self, Self) {
         glam_assert!(self.is_normalized());
         // From https://graphics.pixar.com/library/OrthonormalB/paper.pdf
-        #[cfg(feature = "std")]
-        let sign = (1.0_f32).copysign(self.z);
-        #[cfg(not(feature = "std"))]
-        let sign = self.z.signum();
+        let sign = math::signum(self.z);
         let a = -1.0 / (sign + self.z);
         let b = self.x * self.y * a;
         (
diff --git a/src/f32/wasm32/vec4.rs b/src/f32/wasm32/vec4.rs
index 927a4f8..1b61be6 100644
--- a/src/f32/wasm32/vec4.rs
+++ b/src/f32/wasm32/vec4.rs
@@ -1,6 +1,6 @@
 // Generated from vec.rs.tera template. Edit the template, not the generated file.
 
-use crate::{wasm32::*, BVec4A, Vec2, Vec3, Vec3A};
+use crate::{f32::math, wasm32::*, BVec4A, Vec2, Vec3, Vec3A};
 
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
@@ -9,10 +9,6 @@ use core::{f32, ops::*};
 
 use core::arch::wasm32::*;
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 union UnionCast {
     a: [f32; 4],
     v: Vec4,
@@ -423,7 +419,7 @@ impl Vec4 {
     #[inline]
     pub fn is_normalized(self) -> bool {
         // TODO: do something with epsilon
-        (self.length_squared() - 1.0).abs() <= 1e-4
+        math::abs(self.length_squared() - 1.0) <= 1e-4
     }
 
     /// Returns the vector projection of `self` onto `rhs`.
@@ -521,17 +517,22 @@ impl Vec4 {
     /// `self`.
     #[inline]
     pub fn exp(self) -> Self {
-        Self::new(self.x.exp(), self.y.exp(), self.z.exp(), self.w.exp())
+        Self::new(
+            math::exp(self.x),
+            math::exp(self.y),
+            math::exp(self.z),
+            math::exp(self.w),
+        )
     }
 
     /// Returns a vector containing each element of `self` raised to the power of `n`.
     #[inline]
     pub fn powf(self, n: f32) -> Self {
         Self::new(
-            self.x.powf(n),
-            self.y.powf(n),
-            self.z.powf(n),
-            self.w.powf(n),
+            math::powf(self.x, n),
+            math::powf(self.y, n),
+            math::powf(self.z, n),
+            math::powf(self.w, n),
         )
     }
 
@@ -576,9 +577,9 @@ impl Vec4 {
         glam_assert!(min <= max);
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -588,7 +589,7 @@ impl Vec4 {
     pub fn clamp_length_max(self, max: f32) -> Self {
         let length_sq = self.length_squared();
         if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -598,7 +599,7 @@ impl Vec4 {
     pub fn clamp_length_min(self, min: f32) -> Self {
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -614,10 +615,10 @@ impl Vec4 {
     #[inline]
     pub fn mul_add(self, a: Self, b: Self) -> Self {
         Self::new(
-            self.x.mul_add(a.x, b.x),
-            self.y.mul_add(a.y, b.y),
-            self.z.mul_add(a.z, b.z),
-            self.w.mul_add(a.w, b.w),
+            math::mul_add(self.x, a.x, b.x),
+            math::mul_add(self.y, a.y, b.y),
+            math::mul_add(self.z, a.z, b.z),
+            math::mul_add(self.w, a.w, b.w),
         )
     }
 
diff --git a/src/f64.rs b/src/f64.rs
index 307cd1a..3d05ede 100644
--- a/src/f64.rs
+++ b/src/f64.rs
@@ -7,6 +7,7 @@ mod dquat;
 mod dvec2;
 mod dvec3;
 mod dvec4;
+pub(crate) mod math;
 
 pub use daffine2::DAffine2;
 pub use daffine3::DAffine3;
diff --git a/src/f64/daffine3.rs b/src/f64/daffine3.rs
index a5f4535..c0a56a0 100644
--- a/src/f64/daffine3.rs
+++ b/src/f64/daffine3.rs
@@ -264,15 +264,12 @@ impl DAffine3 {
     /// vector contains any zero elements when `glam_assert` is enabled.
     #[inline]
     pub fn to_scale_rotation_translation(&self) -> (DVec3, DQuat, DVec3) {
-        #[cfg(feature = "libm")]
-        #[allow(unused_imports)]
-        use num_traits::Float;
-
+        use crate::f64::math;
         let det = self.matrix3.determinant();
         glam_assert!(det != 0.0);
 
         let scale = DVec3::new(
-            self.matrix3.x_axis.length() * det.signum(),
+            self.matrix3.x_axis.length() * math::signum(det),
             self.matrix3.y_axis.length(),
             self.matrix3.z_axis.length(),
         );
diff --git a/src/f64/dmat2.rs b/src/f64/dmat2.rs
index bd52c8c..daa7d91 100644
--- a/src/f64/dmat2.rs
+++ b/src/f64/dmat2.rs
@@ -1,15 +1,11 @@
 // Generated from mat.rs.tera template. Edit the template, not the generated file.
 
-use crate::{swizzles::*, DMat3, DVec2, Mat2};
+use crate::{f64::math, swizzles::*, DMat3, DVec2, Mat2};
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
 use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign};
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 2x2 matrix from column vectors.
 #[inline(always)]
 pub const fn dmat2(x_axis: DVec2, y_axis: DVec2) -> DMat2 {
@@ -91,14 +87,14 @@ impl DMat2 {
     /// `angle` (in radians).
     #[inline]
     pub fn from_scale_angle(scale: DVec2, angle: f64) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::new(cos * scale.x, sin * scale.x, -sin * scale.y, cos * scale.y)
     }
 
     /// Creates a 2x2 matrix containing a rotation of `angle` (in radians).
     #[inline]
     pub fn from_angle(angle: f64) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::new(cos, sin, -sin, cos)
     }
 
diff --git a/src/f64/dmat3.rs b/src/f64/dmat3.rs
index 600223c..926fe03 100644
--- a/src/f64/dmat3.rs
+++ b/src/f64/dmat3.rs
@@ -1,15 +1,11 @@
 // Generated from mat.rs.tera template. Edit the template, not the generated file.
 
-use crate::{swizzles::*, DMat2, DMat4, DQuat, DVec2, DVec3, EulerRot, Mat3};
+use crate::{f64::math, swizzles::*, DMat2, DMat4, DQuat, DVec2, DVec3, EulerRot, Mat3};
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
 use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign};
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 3x3 matrix from column vectors.
 #[inline(always)]
 pub const fn dmat3(x_axis: DVec3, y_axis: DVec3, z_axis: DVec3) -> DMat3 {
@@ -189,7 +185,7 @@ impl DMat3 {
     pub fn from_axis_angle(axis: DVec3, angle: f64) -> Self {
         glam_assert!(axis.is_normalized());
 
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         let (xsin, ysin, zsin) = axis.mul(sin).into();
         let (x, y, z) = axis.into();
         let (x2, y2, z2) = axis.mul(axis).into();
@@ -215,7 +211,7 @@ impl DMat3 {
     /// Creates a 3D rotation matrix from `angle` (in radians) around the x axis.
     #[inline]
     pub fn from_rotation_x(angle: f64) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             DVec3::X,
             DVec3::new(0.0, cosa, sina),
@@ -226,7 +222,7 @@ impl DMat3 {
     /// Creates a 3D rotation matrix from `angle` (in radians) around the y axis.
     #[inline]
     pub fn from_rotation_y(angle: f64) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             DVec3::new(cosa, 0.0, -sina),
             DVec3::Y,
@@ -237,7 +233,7 @@ impl DMat3 {
     /// Creates a 3D rotation matrix from `angle` (in radians) around the z axis.
     #[inline]
     pub fn from_rotation_z(angle: f64) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             DVec3::new(cosa, sina, 0.0),
             DVec3::new(-sina, cosa, 0.0),
@@ -265,7 +261,7 @@ impl DMat3 {
     /// [`Self::transform_point2()`] and [`Self::transform_vector2()`].
     #[inline]
     pub fn from_angle(angle: f64) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::from_cols(
             DVec3::new(cos, sin, 0.0),
             DVec3::new(-sin, cos, 0.0),
@@ -280,7 +276,7 @@ impl DMat3 {
     /// [`Self::transform_point2()`] and [`Self::transform_vector2()`].
     #[inline]
     pub fn from_scale_angle_translation(scale: DVec2, angle: f64, translation: DVec2) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self::from_cols(
             DVec3::new(cos * scale.x, sin * scale.x, 0.0),
             DVec3::new(-sin * scale.y, cos * scale.y, 0.0),
diff --git a/src/f64/dmat4.rs b/src/f64/dmat4.rs
index cc4283c..3e2e4de 100644
--- a/src/f64/dmat4.rs
+++ b/src/f64/dmat4.rs
@@ -1,15 +1,11 @@
 // Generated from mat.rs.tera template. Edit the template, not the generated file.
 
-use crate::{swizzles::*, DMat3, DQuat, DVec3, DVec4, EulerRot, Mat4};
+use crate::{f64::math, swizzles::*, DMat3, DQuat, DVec3, DVec4, EulerRot, Mat4};
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
 use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign};
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 4x4 matrix from column vectors.
 #[inline(always)]
 pub const fn dmat4(x_axis: DVec4, y_axis: DVec4, z_axis: DVec4, w_axis: DVec4) -> DMat4 {
@@ -249,7 +245,7 @@ impl DMat4 {
         glam_assert!(det != 0.0);
 
         let scale = DVec3::new(
-            self.x_axis.length() * det.signum(),
+            self.x_axis.length() * math::signum(det),
             self.y_axis.length(),
             self.z_axis.length(),
         );
@@ -325,7 +321,7 @@ impl DMat4 {
     pub fn from_axis_angle(axis: DVec3, angle: f64) -> Self {
         glam_assert!(axis.is_normalized());
 
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         let axis_sin = axis.mul(sin);
         let axis_sq = axis.mul(axis);
         let omc = 1.0 - cos;
@@ -373,7 +369,7 @@ impl DMat4 {
     /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
     #[inline]
     pub fn from_rotation_x(angle: f64) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             DVec4::X,
             DVec4::new(0.0, cosa, sina, 0.0),
@@ -389,7 +385,7 @@ impl DMat4 {
     /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
     #[inline]
     pub fn from_rotation_y(angle: f64) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             DVec4::new(cosa, 0.0, -sina, 0.0),
             DVec4::Y,
@@ -405,7 +401,7 @@ impl DMat4 {
     /// [`Self::transform_point3()`] and [`Self::transform_vector3()`].
     #[inline]
     pub fn from_rotation_z(angle: f64) -> Self {
-        let (sina, cosa) = angle.sin_cos();
+        let (sina, cosa) = math::sin_cos(angle);
         Self::from_cols(
             DVec4::new(cosa, sina, 0.0, 0.0),
             DVec4::new(-sina, cosa, 0.0, 0.0),
@@ -714,7 +710,7 @@ impl DMat4 {
         z_far: f64,
     ) -> Self {
         let inv_length = 1.0 / (z_near - z_far);
-        let f = 1.0 / (0.5 * fov_y_radians).tan();
+        let f = 1.0 / math::tan(0.5 * fov_y_radians);
         let a = f / aspect_ratio;
         let b = (z_near + z_far) * inv_length;
         let c = (2.0 * z_near * z_far) * inv_length;
@@ -735,7 +731,7 @@ impl DMat4 {
     #[inline]
     pub fn perspective_lh(fov_y_radians: f64, aspect_ratio: f64, z_near: f64, z_far: f64) -> Self {
         glam_assert!(z_near > 0.0 && z_far > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         let r = z_far / (z_far - z_near);
@@ -756,7 +752,7 @@ impl DMat4 {
     #[inline]
     pub fn perspective_rh(fov_y_radians: f64, aspect_ratio: f64, z_near: f64, z_far: f64) -> Self {
         glam_assert!(z_near > 0.0 && z_far > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         let r = z_far / (z_near - z_far);
@@ -776,7 +772,7 @@ impl DMat4 {
     #[inline]
     pub fn perspective_infinite_lh(fov_y_radians: f64, aspect_ratio: f64, z_near: f64) -> Self {
         glam_assert!(z_near > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         Self::from_cols(
@@ -799,7 +795,7 @@ impl DMat4 {
         z_near: f64,
     ) -> Self {
         glam_assert!(z_near > 0.0);
-        let (sin_fov, cos_fov) = (0.5 * fov_y_radians).sin_cos();
+        let (sin_fov, cos_fov) = math::sin_cos(0.5 * fov_y_radians);
         let h = cos_fov / sin_fov;
         let w = h / aspect_ratio;
         Self::from_cols(
@@ -815,7 +811,7 @@ impl DMat4 {
     #[inline]
     pub fn perspective_infinite_rh(fov_y_radians: f64, aspect_ratio: f64, z_near: f64) -> Self {
         glam_assert!(z_near > 0.0);
-        let f = 1.0 / (0.5 * fov_y_radians).tan();
+        let f = 1.0 / math::tan(0.5 * fov_y_radians);
         Self::from_cols(
             DVec4::new(f / aspect_ratio, 0.0, 0.0, 0.0),
             DVec4::new(0.0, f, 0.0, 0.0),
@@ -833,7 +829,7 @@ impl DMat4 {
         z_near: f64,
     ) -> Self {
         glam_assert!(z_near > 0.0);
-        let f = 1.0 / (0.5 * fov_y_radians).tan();
+        let f = 1.0 / math::tan(0.5 * fov_y_radians);
         Self::from_cols(
             DVec4::new(f / aspect_ratio, 0.0, 0.0, 0.0),
             DVec4::new(0.0, f, 0.0, 0.0),
diff --git a/src/f64/dquat.rs b/src/f64/dquat.rs
index 66a175c..c8dd7e3 100644
--- a/src/f64/dquat.rs
+++ b/src/f64/dquat.rs
@@ -2,13 +2,10 @@
 
 use crate::{
     euler::{EulerFromQuaternion, EulerRot, EulerToQuaternion},
-    DMat3, DMat4, DVec2, DVec3, DVec4, FloatEx, Quat,
+    f64::math,
+    DMat3, DMat4, DVec2, DVec3, DVec4, Quat,
 };
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
@@ -128,7 +125,7 @@ impl DQuat {
     #[inline]
     pub fn from_axis_angle(axis: DVec3, angle: f64) -> Self {
         glam_assert!(axis.is_normalized());
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         let v = axis * s;
         Self::from_xyzw(v.x, v.y, v.z, c)
     }
@@ -149,21 +146,21 @@ impl DQuat {
     /// Creates a quaternion from the `angle` (in radians) around the x axis.
     #[inline]
     pub fn from_rotation_x(angle: f64) -> Self {
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         Self::from_xyzw(s, 0.0, 0.0, c)
     }
 
     /// Creates a quaternion from the `angle` (in radians) around the y axis.
     #[inline]
     pub fn from_rotation_y(angle: f64) -> Self {
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         Self::from_xyzw(0.0, s, 0.0, c)
     }
 
     /// Creates a quaternion from the `angle` (in radians) around the z axis.
     #[inline]
     pub fn from_rotation_z(angle: f64) -> Self {
-        let (s, c) = (angle * 0.5).sin_cos();
+        let (s, c) = math::sin_cos(angle * 0.5);
         Self::from_xyzw(0.0, 0.0, s, c)
     }
 
@@ -187,7 +184,7 @@ impl DQuat {
             if dif10 <= 0.0 {
                 // x^2 >= y^2
                 let four_xsq = omm22 - dif10;
-                let inv4x = 0.5 / four_xsq.sqrt();
+                let inv4x = 0.5 / math::sqrt(four_xsq);
                 Self::from_xyzw(
                     four_xsq * inv4x,
                     (m01 + m10) * inv4x,
@@ -197,7 +194,7 @@ impl DQuat {
             } else {
                 // y^2 >= x^2
                 let four_ysq = omm22 + dif10;
-                let inv4y = 0.5 / four_ysq.sqrt();
+                let inv4y = 0.5 / math::sqrt(four_ysq);
                 Self::from_xyzw(
                     (m01 + m10) * inv4y,
                     four_ysq * inv4y,
@@ -212,7 +209,7 @@ impl DQuat {
             if sum10 <= 0.0 {
                 // z^2 >= w^2
                 let four_zsq = opm22 - sum10;
-                let inv4z = 0.5 / four_zsq.sqrt();
+                let inv4z = 0.5 / math::sqrt(four_zsq);
                 Self::from_xyzw(
                     (m02 + m20) * inv4z,
                     (m12 + m21) * inv4z,
@@ -222,7 +219,7 @@ impl DQuat {
             } else {
                 // w^2 >= z^2
                 let four_wsq = opm22 + sum10;
-                let inv4w = 0.5 / four_wsq.sqrt();
+                let inv4w = 0.5 / math::sqrt(four_wsq);
                 Self::from_xyzw(
                     (m12 - m21) * inv4w,
                     (m20 - m02) * inv4w,
@@ -336,7 +333,7 @@ impl DQuat {
             let z = from.x * to.y - to.x * from.y;
             let w = 1.0 + dot;
             // calculate length with x=0 and y=0 to normalize
-            let len_rcp = 1.0 / (z * z + w * w).sqrt();
+            let len_rcp = 1.0 / math::sqrt(z * z + w * w);
             Self::from_xyzw(0.0, 0.0, z * len_rcp, w * len_rcp)
         }
     }
@@ -348,7 +345,7 @@ impl DQuat {
         let v = DVec3::new(self.x, self.y, self.z);
         let length = v.length();
         if length >= EPSILON {
-            let angle = 2.0 * f64::atan2(length, self.w);
+            let angle = 2.0 * math::atan2(length, self.w);
             let axis = v / length;
             (axis, angle)
         } else {
@@ -492,7 +489,7 @@ impl DQuat {
         // If the quat.w is close to -1.0, the angle will be near 2*PI which is close to
         // a negative 0 rotation. By forcing quat.w to be positive, we'll end up with
         // the shortest path.
-        let positive_w_angle = self.w.abs().acos_approx() * 2.0;
+        let positive_w_angle = math::acos_approx(math::abs(self.w)) * 2.0;
         positive_w_angle < threshold_angle
     }
 
@@ -507,7 +504,7 @@ impl DQuat {
     #[inline]
     pub fn angle_between(self, rhs: Self) -> f64 {
         glam_assert!(self.is_normalized() && rhs.is_normalized());
-        self.dot(rhs).abs().acos_approx() * 2.0
+        math::acos_approx(math::abs(self.dot(rhs))) * 2.0
     }
 
     /// Returns true if the absolute difference of all elements between `self` and `rhs`
@@ -579,13 +576,13 @@ impl DQuat {
             // assumes lerp returns a normalized quaternion
             self.lerp(end, s)
         } else {
-            let theta = dot.acos_approx();
+            let theta = math::acos_approx(dot);
 
-            let scale1 = (theta * (1.0 - s)).sin();
-            let scale2 = (theta * s).sin();
-            let theta_sin = theta.sin();
+            let scale1 = math::sin(theta * (1.0 - s));
+            let scale2 = math::sin(theta * s);
+            let theta_sin = math::sin(theta);
 
-            self.mul(scale1).add(end.mul(scale2)).mul(theta_sin.recip())
+            self.mul(scale1).add(end.mul(scale2)).mul(1.0 / theta_sin)
         }
     }
 
diff --git a/src/f64/dvec2.rs b/src/f64/dvec2.rs
index 01b0f22..47f6701 100644
--- a/src/f64/dvec2.rs
+++ b/src/f64/dvec2.rs
@@ -1,16 +1,12 @@
 // Generated from vec.rs.tera template. Edit the template, not the generated file.
 
-use crate::{BVec2, DVec3};
+use crate::{f64::math, BVec2, DVec3};
 
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
 use core::{f32, ops::*};
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 2-dimensional vector.
 #[inline(always)]
 pub const fn dvec2(x: f64, y: f64) -> DVec2 {
@@ -246,8 +242,8 @@ impl DVec2 {
     #[inline]
     pub fn abs(self) -> Self {
         Self {
-            x: self.x.abs(),
-            y: self.y.abs(),
+            x: math::abs(self.x),
+            y: math::abs(self.y),
         }
     }
 
@@ -259,8 +255,8 @@ impl DVec2 {
     #[inline]
     pub fn signum(self) -> Self {
         Self {
-            x: self.x.signum(),
-            y: self.y.signum(),
+            x: math::signum(self.x),
+            y: math::signum(self.y),
         }
     }
 
@@ -268,8 +264,8 @@ impl DVec2 {
     #[inline]
     pub fn copysign(self, rhs: Self) -> Self {
         Self {
-            x: self.x.copysign(rhs.x),
-            y: self.y.copysign(rhs.y),
+            x: math::copysign(self.x, rhs.x),
+            y: math::copysign(self.y, rhs.y),
         }
     }
 
@@ -307,7 +303,7 @@ impl DVec2 {
     #[doc(alias = "magnitude")]
     #[inline]
     pub fn length(self) -> f64 {
-        self.dot(self).sqrt()
+        math::sqrt(self.dot(self))
     }
 
     /// Computes the squared length of `self`.
@@ -397,7 +393,7 @@ impl DVec2 {
     #[inline]
     pub fn is_normalized(self) -> bool {
         // TODO: do something with epsilon
-        (self.length_squared() - 1.0).abs() <= 1e-4
+        math::abs(self.length_squared() - 1.0) <= 1e-4
     }
 
     /// Returns the vector projection of `self` onto `rhs`.
@@ -466,8 +462,8 @@ impl DVec2 {
     #[inline]
     pub fn round(self) -> Self {
         Self {
-            x: self.x.round(),
-            y: self.y.round(),
+            x: math::round(self.x),
+            y: math::round(self.y),
         }
     }
 
@@ -476,8 +472,8 @@ impl DVec2 {
     #[inline]
     pub fn floor(self) -> Self {
         Self {
-            x: self.x.floor(),
-            y: self.y.floor(),
+            x: math::floor(self.x),
+            y: math::floor(self.y),
         }
     }
 
@@ -486,8 +482,8 @@ impl DVec2 {
     #[inline]
     pub fn ceil(self) -> Self {
         Self {
-            x: self.x.ceil(),
-            y: self.y.ceil(),
+            x: math::ceil(self.x),
+            y: math::ceil(self.y),
         }
     }
 
@@ -504,21 +500,21 @@ impl DVec2 {
     /// `self`.
     #[inline]
     pub fn exp(self) -> Self {
-        Self::new(self.x.exp(), self.y.exp())
+        Self::new(math::exp(self.x), math::exp(self.y))
     }
 
     /// Returns a vector containing each element of `self` raised to the power of `n`.
     #[inline]
     pub fn powf(self, n: f64) -> Self {
-        Self::new(self.x.powf(n), self.y.powf(n))
+        Self::new(math::powf(self.x, n), math::powf(self.y, n))
     }
 
     /// Returns a vector containing the reciprocal `1.0/n` of each element of `self`.
     #[inline]
     pub fn recip(self) -> Self {
         Self {
-            x: self.x.recip(),
-            y: self.y.recip(),
+            x: 1.0 / self.x,
+            y: 1.0 / self.y,
         }
     }
 
@@ -557,9 +553,9 @@ impl DVec2 {
         glam_assert!(min <= max);
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -569,7 +565,7 @@ impl DVec2 {
     pub fn clamp_length_max(self, max: f64) -> Self {
         let length_sq = self.length_squared();
         if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -579,7 +575,7 @@ impl DVec2 {
     pub fn clamp_length_min(self, min: f64) -> Self {
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -594,7 +590,10 @@ impl DVec2 {
     /// mind.
     #[inline]
     pub fn mul_add(self, a: Self, b: Self) -> Self {
-        Self::new(self.x.mul_add(a.x, b.x), self.y.mul_add(a.y, b.y))
+        Self::new(
+            math::mul_add(self.x, a.x, b.x),
+            math::mul_add(self.y, a.y, b.y),
+        )
     }
 
     /// Creates a 2D vector containing `[angle.cos(), angle.sin()]`. This can be used in
@@ -602,7 +601,7 @@ impl DVec2 {
     /// create the vector [-1, 0] and rotate `Vec2::Y` around it returning `-Vec2::Y`.
     #[inline]
     pub fn from_angle(angle: f64) -> Self {
-        let (sin, cos) = angle.sin_cos();
+        let (sin, cos) = math::sin_cos(angle);
         Self { x: cos, y: sin }
     }
 
@@ -611,11 +610,11 @@ impl DVec2 {
     /// The input vectors do not need to be unit length however they must be non-zero.
     #[inline]
     pub fn angle_between(self, rhs: Self) -> f64 {
-        use crate::FloatEx;
-        let angle =
-            (self.dot(rhs) / (self.length_squared() * rhs.length_squared()).sqrt()).acos_approx();
+        let angle = math::acos_approx(
+            self.dot(rhs) / math::sqrt(self.length_squared() * rhs.length_squared()),
+        );
 
-        angle * self.perp_dot(rhs).signum()
+        angle * math::signum(self.perp_dot(rhs))
     }
 
     /// Returns a vector that is equal to `self` rotated by 90 degrees.
diff --git a/src/f64/dvec3.rs b/src/f64/dvec3.rs
index 7493fd2..4d6dccc 100644
--- a/src/f64/dvec3.rs
+++ b/src/f64/dvec3.rs
@@ -1,16 +1,12 @@
 // Generated from vec.rs.tera template. Edit the template, not the generated file.
 
-use crate::{BVec3, DVec2, DVec4};
+use crate::{f64::math, BVec3, DVec2, DVec4};
 
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
 use core::{f32, ops::*};
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 3-dimensional vector.
 #[inline(always)]
 pub const fn dvec3(x: f64, y: f64, z: f64) -> DVec3 {
@@ -286,9 +282,9 @@ impl DVec3 {
     #[inline]
     pub fn abs(self) -> Self {
         Self {
-            x: self.x.abs(),
-            y: self.y.abs(),
-            z: self.z.abs(),
+            x: math::abs(self.x),
+            y: math::abs(self.y),
+            z: math::abs(self.z),
         }
     }
 
@@ -300,9 +296,9 @@ impl DVec3 {
     #[inline]
     pub fn signum(self) -> Self {
         Self {
-            x: self.x.signum(),
-            y: self.y.signum(),
-            z: self.z.signum(),
+            x: math::signum(self.x),
+            y: math::signum(self.y),
+            z: math::signum(self.z),
         }
     }
 
@@ -310,9 +306,9 @@ impl DVec3 {
     #[inline]
     pub fn copysign(self, rhs: Self) -> Self {
         Self {
-            x: self.x.copysign(rhs.x),
-            y: self.y.copysign(rhs.y),
-            z: self.z.copysign(rhs.z),
+            x: math::copysign(self.x, rhs.x),
+            y: math::copysign(self.y, rhs.y),
+            z: math::copysign(self.z, rhs.z),
         }
     }
 
@@ -352,7 +348,7 @@ impl DVec3 {
     #[doc(alias = "magnitude")]
     #[inline]
     pub fn length(self) -> f64 {
-        self.dot(self).sqrt()
+        math::sqrt(self.dot(self))
     }
 
     /// Computes the squared length of `self`.
@@ -442,7 +438,7 @@ impl DVec3 {
     #[inline]
     pub fn is_normalized(self) -> bool {
         // TODO: do something with epsilon
-        (self.length_squared() - 1.0).abs() <= 1e-4
+        math::abs(self.length_squared() - 1.0) <= 1e-4
     }
 
     /// Returns the vector projection of `self` onto `rhs`.
@@ -511,9 +507,9 @@ impl DVec3 {
     #[inline]
     pub fn round(self) -> Self {
         Self {
-            x: self.x.round(),
-            y: self.y.round(),
-            z: self.z.round(),
+            x: math::round(self.x),
+            y: math::round(self.y),
+            z: math::round(self.z),
         }
     }
 
@@ -522,9 +518,9 @@ impl DVec3 {
     #[inline]
     pub fn floor(self) -> Self {
         Self {
-            x: self.x.floor(),
-            y: self.y.floor(),
-            z: self.z.floor(),
+            x: math::floor(self.x),
+            y: math::floor(self.y),
+            z: math::floor(self.z),
         }
     }
 
@@ -533,9 +529,9 @@ impl DVec3 {
     #[inline]
     pub fn ceil(self) -> Self {
         Self {
-            x: self.x.ceil(),
-            y: self.y.ceil(),
-            z: self.z.ceil(),
+            x: math::ceil(self.x),
+            y: math::ceil(self.y),
+            z: math::ceil(self.z),
         }
     }
 
@@ -552,22 +548,26 @@ impl DVec3 {
     /// `self`.
     #[inline]
     pub fn exp(self) -> Self {
-        Self::new(self.x.exp(), self.y.exp(), self.z.exp())
+        Self::new(math::exp(self.x), math::exp(self.y), math::exp(self.z))
     }
 
     /// Returns a vector containing each element of `self` raised to the power of `n`.
     #[inline]
     pub fn powf(self, n: f64) -> Self {
-        Self::new(self.x.powf(n), self.y.powf(n), self.z.powf(n))
+        Self::new(
+            math::powf(self.x, n),
+            math::powf(self.y, n),
+            math::powf(self.z, n),
+        )
     }
 
     /// Returns a vector containing the reciprocal `1.0/n` of each element of `self`.
     #[inline]
     pub fn recip(self) -> Self {
         Self {
-            x: self.x.recip(),
-            y: self.y.recip(),
-            z: self.z.recip(),
+            x: 1.0 / self.x,
+            y: 1.0 / self.y,
+            z: 1.0 / self.z,
         }
     }
 
@@ -606,9 +606,9 @@ impl DVec3 {
         glam_assert!(min <= max);
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -618,7 +618,7 @@ impl DVec3 {
     pub fn clamp_length_max(self, max: f64) -> Self {
         let length_sq = self.length_squared();
         if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -628,7 +628,7 @@ impl DVec3 {
     pub fn clamp_length_min(self, min: f64) -> Self {
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -644,9 +644,9 @@ impl DVec3 {
     #[inline]
     pub fn mul_add(self, a: Self, b: Self) -> Self {
         Self::new(
-            self.x.mul_add(a.x, b.x),
-            self.y.mul_add(a.y, b.y),
-            self.z.mul_add(a.z, b.z),
+            math::mul_add(self.x, a.x, b.x),
+            math::mul_add(self.y, a.y, b.y),
+            math::mul_add(self.z, a.z, b.z),
         )
     }
 
@@ -655,10 +655,10 @@ impl DVec3 {
     /// The input vectors do not need to be unit length however they must be non-zero.
     #[inline]
     pub fn angle_between(self, rhs: Self) -> f64 {
-        use crate::FloatEx;
-        self.dot(rhs)
-            .div(self.length_squared().mul(rhs.length_squared()).sqrt())
-            .acos_approx()
+        math::acos_approx(
+            self.dot(rhs)
+                .div(math::sqrt(self.length_squared().mul(rhs.length_squared()))),
+        )
     }
 
     /// Returns some vector that is orthogonal to the given one.
@@ -670,7 +670,7 @@ impl DVec3 {
     #[inline]
     pub fn any_orthogonal_vector(&self) -> Self {
         // This can probably be optimized
-        if self.x.abs() > self.y.abs() {
+        if math::abs(self.x) > math::abs(self.y) {
             Self::new(-self.z, 0.0, self.x) // self.cross(Self::Y)
         } else {
             Self::new(0.0, self.z, -self.y) // self.cross(Self::X)
@@ -687,10 +687,7 @@ impl DVec3 {
     pub fn any_orthonormal_vector(&self) -> Self {
         glam_assert!(self.is_normalized());
         // From https://graphics.pixar.com/library/OrthonormalB/paper.pdf
-        #[cfg(feature = "std")]
-        let sign = (1.0_f64).copysign(self.z);
-        #[cfg(not(feature = "std"))]
-        let sign = self.z.signum();
+        let sign = math::signum(self.z);
         let a = -1.0 / (sign + self.z);
         let b = self.x * self.y * a;
         Self::new(b, sign + self.y * self.y * a, -self.y)
@@ -706,10 +703,7 @@ impl DVec3 {
     pub fn any_orthonormal_pair(&self) -> (Self, Self) {
         glam_assert!(self.is_normalized());
         // From https://graphics.pixar.com/library/OrthonormalB/paper.pdf
-        #[cfg(feature = "std")]
-        let sign = (1.0_f64).copysign(self.z);
-        #[cfg(not(feature = "std"))]
-        let sign = self.z.signum();
+        let sign = math::signum(self.z);
         let a = -1.0 / (sign + self.z);
         let b = self.x * self.y * a;
         (
diff --git a/src/f64/dvec4.rs b/src/f64/dvec4.rs
index 4683a97..fc9423c 100644
--- a/src/f64/dvec4.rs
+++ b/src/f64/dvec4.rs
@@ -1,16 +1,12 @@
 // Generated from vec.rs.tera template. Edit the template, not the generated file.
 
-use crate::{BVec4, DVec2, DVec3};
+use crate::{f64::math, BVec4, DVec2, DVec3};
 
 #[cfg(not(target_arch = "spirv"))]
 use core::fmt;
 use core::iter::{Product, Sum};
 use core::{f32, ops::*};
 
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
 /// Creates a 4-dimensional vector.
 #[inline(always)]
 pub const fn dvec4(x: f64, y: f64, z: f64, w: f64) -> DVec4 {
@@ -309,10 +305,10 @@ impl DVec4 {
     #[inline]
     pub fn abs(self) -> Self {
         Self {
-            x: self.x.abs(),
-            y: self.y.abs(),
-            z: self.z.abs(),
-            w: self.w.abs(),
+            x: math::abs(self.x),
+            y: math::abs(self.y),
+            z: math::abs(self.z),
+            w: math::abs(self.w),
         }
     }
 
@@ -324,10 +320,10 @@ impl DVec4 {
     #[inline]
     pub fn signum(self) -> Self {
         Self {
-            x: self.x.signum(),
-            y: self.y.signum(),
-            z: self.z.signum(),
-            w: self.w.signum(),
+            x: math::signum(self.x),
+            y: math::signum(self.y),
+            z: math::signum(self.z),
+            w: math::signum(self.w),
         }
     }
 
@@ -335,10 +331,10 @@ impl DVec4 {
     #[inline]
     pub fn copysign(self, rhs: Self) -> Self {
         Self {
-            x: self.x.copysign(rhs.x),
-            y: self.y.copysign(rhs.y),
-            z: self.z.copysign(rhs.z),
-            w: self.w.copysign(rhs.w),
+            x: math::copysign(self.x, rhs.x),
+            y: math::copysign(self.y, rhs.y),
+            z: math::copysign(self.z, rhs.z),
+            w: math::copysign(self.w, rhs.w),
         }
     }
 
@@ -384,7 +380,7 @@ impl DVec4 {
     #[doc(alias = "magnitude")]
     #[inline]
     pub fn length(self) -> f64 {
-        self.dot(self).sqrt()
+        math::sqrt(self.dot(self))
     }
 
     /// Computes the squared length of `self`.
@@ -474,7 +470,7 @@ impl DVec4 {
     #[inline]
     pub fn is_normalized(self) -> bool {
         // TODO: do something with epsilon
-        (self.length_squared() - 1.0).abs() <= 1e-4
+        math::abs(self.length_squared() - 1.0) <= 1e-4
     }
 
     /// Returns the vector projection of `self` onto `rhs`.
@@ -543,10 +539,10 @@ impl DVec4 {
     #[inline]
     pub fn round(self) -> Self {
         Self {
-            x: self.x.round(),
-            y: self.y.round(),
-            z: self.z.round(),
-            w: self.w.round(),
+            x: math::round(self.x),
+            y: math::round(self.y),
+            z: math::round(self.z),
+            w: math::round(self.w),
         }
     }
 
@@ -555,10 +551,10 @@ impl DVec4 {
     #[inline]
     pub fn floor(self) -> Self {
         Self {
-            x: self.x.floor(),
-            y: self.y.floor(),
-            z: self.z.floor(),
-            w: self.w.floor(),
+            x: math::floor(self.x),
+            y: math::floor(self.y),
+            z: math::floor(self.z),
+            w: math::floor(self.w),
         }
     }
 
@@ -567,10 +563,10 @@ impl DVec4 {
     #[inline]
     pub fn ceil(self) -> Self {
         Self {
-            x: self.x.ceil(),
-            y: self.y.ceil(),
-            z: self.z.ceil(),
-            w: self.w.ceil(),
+            x: math::ceil(self.x),
+            y: math::ceil(self.y),
+            z: math::ceil(self.z),
+            w: math::ceil(self.w),
         }
     }
 
@@ -587,17 +583,22 @@ impl DVec4 {
     /// `self`.
     #[inline]
     pub fn exp(self) -> Self {
-        Self::new(self.x.exp(), self.y.exp(), self.z.exp(), self.w.exp())
+        Self::new(
+            math::exp(self.x),
+            math::exp(self.y),
+            math::exp(self.z),
+            math::exp(self.w),
+        )
     }
 
     /// Returns a vector containing each element of `self` raised to the power of `n`.
     #[inline]
     pub fn powf(self, n: f64) -> Self {
         Self::new(
-            self.x.powf(n),
-            self.y.powf(n),
-            self.z.powf(n),
-            self.w.powf(n),
+            math::powf(self.x, n),
+            math::powf(self.y, n),
+            math::powf(self.z, n),
+            math::powf(self.w, n),
         )
     }
 
@@ -605,10 +606,10 @@ impl DVec4 {
     #[inline]
     pub fn recip(self) -> Self {
         Self {
-            x: self.x.recip(),
-            y: self.y.recip(),
-            z: self.z.recip(),
-            w: self.w.recip(),
+            x: 1.0 / self.x,
+            y: 1.0 / self.y,
+            z: 1.0 / self.z,
+            w: 1.0 / self.w,
         }
     }
 
@@ -647,9 +648,9 @@ impl DVec4 {
         glam_assert!(min <= max);
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -659,7 +660,7 @@ impl DVec4 {
     pub fn clamp_length_max(self, max: f64) -> Self {
         let length_sq = self.length_squared();
         if length_sq > max * max {
-            self * (length_sq.sqrt().recip() * max)
+            max * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -669,7 +670,7 @@ impl DVec4 {
     pub fn clamp_length_min(self, min: f64) -> Self {
         let length_sq = self.length_squared();
         if length_sq < min * min {
-            self * (length_sq.sqrt().recip() * min)
+            min * (self / math::sqrt(length_sq))
         } else {
             self
         }
@@ -685,10 +686,10 @@ impl DVec4 {
     #[inline]
     pub fn mul_add(self, a: Self, b: Self) -> Self {
         Self::new(
-            self.x.mul_add(a.x, b.x),
-            self.y.mul_add(a.y, b.y),
-            self.z.mul_add(a.z, b.z),
-            self.w.mul_add(a.w, b.w),
+            math::mul_add(self.x, a.x, b.x),
+            math::mul_add(self.y, a.y, b.y),
+            math::mul_add(self.z, a.z, b.z),
+            math::mul_add(self.w, a.w, b.w),
         )
     }
 
diff --git a/src/f64/math.rs b/src/f64/math.rs
new file mode 100644
index 0000000..309be33
--- /dev/null
+++ b/src/f64/math.rs
@@ -0,0 +1,181 @@
+#[cfg(feature = "libm")]
+mod libm_math {
+    #[inline(always)]
+    pub(crate) fn abs(f: f64) -> f64 {
+        if f.is_sign_positive() {
+            return f;
+        }
+        if f.is_sign_negative() {
+            return -f;
+        }
+        f64::NAN
+    }
+
+    #[inline(always)]
+    pub(crate) fn acos_approx(f: f64) -> f64 {
+        libm::acos(f.clamp(-1.0, 1.0))
+    }
+
+    #[inline(always)]
+    pub(crate) fn asin_clamped(f: f64) -> f64 {
+        libm::asin(f.clamp(-1.0, 1.0))
+    }
+
+    #[inline(always)]
+    pub(crate) fn atan2(f: f64, other: f64) -> f64 {
+        libm::atan2(f, other)
+    }
+
+    #[inline(always)]
+    pub(crate) fn sin(f: f64) -> f64 {
+        libm::sin(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn sin_cos(f: f64) -> (f64, f64) {
+        libm::sincos(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn tan(f: f64) -> f64 {
+        libm::tan(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn sqrt(f: f64) -> f64 {
+        libm::sqrt(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn copysign(f: f64, sign: f64) -> f64 {
+        libm::copysign(f, sign)
+    }
+
+    #[inline(always)]
+    pub(crate) fn signum(f: f64) -> f64 {
+        if f.is_nan() {
+            f64::NAN
+        } else {
+            copysign(1.0, f)
+        }
+    }
+
+    #[inline(always)]
+    pub(crate) fn round(f: f64) -> f64 {
+        libm::round(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn ceil(f: f64) -> f64 {
+        libm::ceil(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn floor(f: f64) -> f64 {
+        libm::floor(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn exp(f: f64) -> f64 {
+        libm::exp(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn powf(f: f64, n: f64) -> f64 {
+        libm::pow(f, n)
+    }
+
+    #[inline(always)]
+    pub(crate) fn mul_add(a: f64, b: f64, c: f64) -> f64 {
+        libm::fma(a, b, c)
+    }
+}
+
+#[cfg(not(feature = "libm"))]
+mod std_math {
+    #[inline(always)]
+    pub(crate) fn abs(f: f64) -> f64 {
+        f64::abs(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn acos_approx(f: f64) -> f64 {
+        f64::acos(f64::clamp(f, -1.0, 1.0))
+    }
+
+    #[inline(always)]
+    pub(crate) fn asin_clamped(f: f64) -> f64 {
+        f64::asin(f64::clamp(f, -1.0, 1.0))
+    }
+
+    #[inline(always)]
+    pub(crate) fn atan2(f: f64, other: f64) -> f64 {
+        f64::atan2(f, other)
+    }
+
+    #[inline(always)]
+    pub(crate) fn sin(f: f64) -> f64 {
+        f64::sin(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn sin_cos(f: f64) -> (f64, f64) {
+        f64::sin_cos(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn tan(f: f64) -> f64 {
+        f64::tan(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn sqrt(f: f64) -> f64 {
+        f64::sqrt(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn copysign(f: f64, sign: f64) -> f64 {
+        f64::copysign(f, sign)
+    }
+
+    #[inline(always)]
+    pub(crate) fn signum(f: f64) -> f64 {
+        f64::signum(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn round(f: f64) -> f64 {
+        f64::round(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn ceil(f: f64) -> f64 {
+        f64::ceil(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn floor(f: f64) -> f64 {
+        f64::floor(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn exp(f: f64) -> f64 {
+        f64::exp(f)
+    }
+
+    #[inline(always)]
+    pub(crate) fn powf(f: f64, n: f64) -> f64 {
+        f64::powf(f, n)
+    }
+
+    #[inline(always)]
+    pub(crate) fn mul_add(a: f64, b: f64, c: f64) -> f64 {
+        f64::mul_add(a, b, c)
+    }
+}
+
+#[cfg(feature = "libm")]
+pub(crate) use libm_math::*;
+
+#[cfg(not(feature = "libm"))]
+pub(crate) use std_math::*;
diff --git a/src/float_ex.rs b/src/float_ex.rs
deleted file mode 100644
index 70f382a..0000000
--- a/src/float_ex.rs
+++ /dev/null
@@ -1,51 +0,0 @@
-#[cfg(feature = "libm")]
-#[allow(unused_imports)]
-use num_traits::Float;
-
-pub(crate) trait FloatEx {
-    /// Returns a very close approximation of `self.clamp(-1.0, 1.0).acos()`.
-    fn acos_approx(self) -> Self;
-}
-
-impl FloatEx for f32 {
-    #[inline(always)]
-    fn acos_approx(self) -> Self {
-        // Based on https://github.com/microsoft/DirectXMath `XMScalarAcos`
-        // Clamp input to [-1,1].
-        let nonnegative = self >= 0.0;
-        let x = self.abs();
-        let mut omx = 1.0 - x;
-        if omx < 0.0 {
-            omx = 0.0;
-        }
-        let root = omx.sqrt();
-
-        // 7-degree minimax approximation
-        #[allow(clippy::approx_constant)]
-        let mut result = ((((((-0.001_262_491_1 * x + 0.006_670_09) * x - 0.017_088_126) * x
-            + 0.030_891_88)
-            * x
-            - 0.050_174_303)
-            * x
-            + 0.088_978_99)
-            * x
-            - 0.214_598_8)
-            * x
-            + 1.570_796_3;
-        result *= root;
-
-        // acos(x) = pi - acos(-x) when x < 0
-        if nonnegative {
-            result
-        } else {
-            core::f32::consts::PI - result
-        }
-    }
-}
-
-impl FloatEx for f64 {
-    #[inline(always)]
-    fn acos_approx(self) -> Self {
-        f64::acos(self.clamp(-1.0, 1.0))
-    }
-}
diff --git a/src/lib.rs b/src/lib.rs
index f9b1185..e556de9 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -218,10 +218,8 @@ and benchmarks.
 * `std` - the default feature, has no dependencies.
 * `approx` - traits and macros for approximate float comparisons
 * `bytemuck` - for casting into slices of bytes
-* `libm` - required to compile with `no_std`
+* `libm` - uses `libm` math functions instead of `std`, required to compile with `no_std`
 * `mint` - for interoperating with other 3D math libraries
-* `num-traits` - required to compile `no_std`, will be included when enabling
-  the `libm` feature
 * `rand` - implementations of `Distribution` trait for all `glam` types.
 * `rkyv` - implementations of `Archive`, `Serialize` and `Deserialize` for all
   `glam` types. Note that serialization is not interoperable with and without the
@@ -272,7 +270,6 @@ mod align16;
 mod deref;
 mod euler;
 mod features;
-mod float_ex;
 
 #[cfg(target_arch = "spirv")]
 mod spirv;
@@ -298,8 +295,6 @@ mod coresimd;
 ))]
 use align16::Align16;
 
-use float_ex::FloatEx;
-
 /** `bool` vector mask types. */
 pub mod bool;
 pub use self::bool::*;


Test patch:
diff --git a/build_and_test_features.sh b/build_and_test_features.sh
index 47d82a5..98c39a6 100755
--- a/build_and_test_features.sh
+++ b/build_and_test_features.sh
@@ -10,6 +10,8 @@ FEATURE_SETS=(
   "std scalar-math approx bytemuck mint rand serde debug-glam-assert"
   "std cuda"
   "std scalar-math cuda"
+  "std libm"
+  "std scalar-math libm"
   # no_std
   "libm"
   "libm scalar-math approx bytemuck mint rand serde debug-glam-assert"
@@ -26,4 +28,4 @@ done
 
 RUSTFLAGS='-C target-feature=+fma' cargo check
 
-cargo check -p glam-no_std
\ No newline at end of file
+cargo check -p glam-no_std
diff --git a/tests/mat3.rs b/tests/mat3.rs
index baf1e42..6119793 100644
--- a/tests/mat3.rs
+++ b/tests/mat3.rs
@@ -193,7 +193,7 @@ macro_rules! impl_mat3_tests {
 
             let yx0 = y0 * x0;
             let yx1 = $mat3::from_euler(EulerRot::YXZ, yaw, pitch, zero);
-            assert_approx_eq!(yx0, yx1);
+            assert_approx_eq!(yx0, yx1, 1e-6);
 
             let yxz0 = y0 * x0 * z0;
             let yxz1 = $mat3::from_euler(EulerRot::YXZ, yaw, pitch, roll);
diff --git a/tests/mat4.rs b/tests/mat4.rs
index 164bdd6..87703f9 100644
--- a/tests/mat4.rs
+++ b/tests/mat4.rs
@@ -242,7 +242,7 @@ macro_rules! impl_mat4_tests {
 
             let yx0 = y0 * x0;
             let yx1 = $mat4::from_euler(EulerRot::YXZ, yaw, pitch, zero);
-            assert_approx_eq!(yx0, yx1);
+            assert_approx_eq!(yx0, yx1, 1e-6);
 
             let yxz0 = y0 * x0 * z0;
             let yxz1 = $mat4::from_euler(EulerRot::YXZ, yaw, pitch, roll);