Repository: bitshifter/glam-rs (rust)

Issue:
Expose FMA via "mul_add" for f32/f64 Vec types.
This should allow users to ensure that a FMA instruction is used (with the SSE2 instruction if possible), which should both improve the accuracy and potentially speed of the resultant code.

Pull request:
Expose mul_add for vector types.
Fixes #269.

This PR does the following:

 * Exposes "mul_add" on all vector types.
 * Updates documentation for the function
 * Forces the scalar implementations to use `f32::mul_add` and `f64::mul_add` for the types where it's relevant.
 * Updates the relevant test names as "mul_add" is no longer private.

Patch:
diff --git a/Cargo.toml b/Cargo.toml
index 7694d96..09131f1 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -36,6 +36,13 @@ libm = ["num-traits/libm"]
 # align types to match CUDA requirements
 cuda = []
 
+# Enables platform specific optimizations that might speed-up certain operations.
+# This will cause APIs to output different results depending on the platform used
+# and will likely break cross-platform determinism.
+# This should NOT be enabled by intermediate libraries, deferring the decision to
+# the end binary build instead.
+fast-math = []
+
 [dependencies]
 approx = { version = "0.5", optional = true, default-features = false }
 bytemuck = { version = "1.5", optional = true, default-features = false }
diff --git a/README.md b/README.md
index 80dd79b..8c2893d 100644
--- a/README.md
+++ b/README.md
@@ -100,6 +100,11 @@ glam = { version = "0.20.2", default-features = false }
   `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
+* [`fast-math`] - By default, glam attempts to provide bit-for-bit identical 
+  results on all platforms. Using this feature will enable platform specific 
+  optimizations that may not be identical to other platforms. **Intermediate 
+  libraries should not use this feature and defer the decision to the final 
+  binary build**.
 
 [`approx`]: https://docs.rs/approx
 [`bytemuck`]: https://docs.rs/bytemuck
diff --git a/src/core/scalar/vector.rs b/src/core/scalar/vector.rs
index 218f03c..4712502 100644
--- a/src/core/scalar/vector.rs
+++ b/src/core/scalar/vector.rs
@@ -210,14 +210,6 @@ impl<T: NumEx> Vector<T> for XY<T> {
         }
     }
 
-    #[inline]
-    fn mul_add(self, b: Self, c: Self) -> Self {
-        Self {
-            x: self.x * b.x + c.x,
-            y: self.y * b.y + c.y,
-        }
-    }
-
     #[inline]
     fn sub(self, other: Self) -> Self {
         Self {
@@ -389,15 +381,6 @@ impl<T: NumEx> Vector<T> for XYZ<T> {
         }
     }
 
-    #[inline]
-    fn mul_add(self, b: Self, c: Self) -> Self {
-        Self {
-            x: self.x * b.x + c.x,
-            y: self.y * b.y + c.y,
-            z: self.z * b.z + c.z,
-        }
-    }
-
     #[inline]
     fn sub(self, other: Self) -> Self {
         Self {
@@ -591,16 +574,6 @@ impl<T: NumEx> Vector<T> for XYZW<T> {
         }
     }
 
-    #[inline]
-    fn mul_add(self, b: Self, c: Self) -> Self {
-        Self {
-            x: self.x * b.x + c.x,
-            y: self.y * b.y + c.y,
-            z: self.z * b.z + c.z,
-            w: self.w * b.w + c.w,
-        }
-    }
-
     #[inline]
     fn sub(self, other: Self) -> Self {
         Self {
@@ -1022,11 +995,6 @@ impl Vector<f32> for XYZF32A16 {
         XYZ::mul(self.into(), other.into()).into()
     }
 
-    #[inline]
-    fn mul_add(self, b: Self, c: Self) -> Self {
-        XYZ::mul_add(self.into(), b.into(), c.into()).into()
-    }
-
     #[inline]
     fn rem(self, other: Self) -> Self {
         XYZ::rem(self.into(), other.into()).into()
diff --git a/src/core/sse2/float.rs b/src/core/sse2/float.rs
index 36a3fa4..9c3ac8e 100644
--- a/src/core/sse2/float.rs
+++ b/src/core/sse2/float.rs
@@ -88,12 +88,14 @@ pub(crate) unsafe fn m128_ceil(v: __m128) -> __m128 {
 
 #[inline(always)]
 pub(crate) unsafe fn m128_mul_add(a: __m128, b: __m128, c: __m128) -> __m128 {
-    #[cfg(target_feature = "fma")]
+    // Only enable fused multiply-adds here if "fast-math" is enabled and the
+    // platform supports it. Otherwise this may break cross-platform determinism.
+    #[cfg(all(feature = "fast-math", target_feature = "fma"))]
     {
         _mm_fmadd_ps(a, b, c)
     }
 
-    #[cfg(not(target_feature = "fma"))]
+    #[cfg(any(not(feature = "fast-math"), not(target_feature = "fma")))]
     {
         _mm_add_ps(_mm_mul_ps(a, b), c)
     }
diff --git a/src/core/sse2/vector.rs b/src/core/sse2/vector.rs
index 9cee05c..d5380ce 100644
--- a/src/core/sse2/vector.rs
+++ b/src/core/sse2/vector.rs
@@ -237,11 +237,6 @@ impl Vector<f32> for __m128 {
         unsafe { _mm_mul_ps(self, other) }
     }
 
-    #[inline(always)]
-    fn mul_add(self, a: Self, b: Self) -> Self {
-        unsafe { m128_mul_add(self, a, b) }
-    }
-
     #[inline(always)]
     fn sub(self, other: Self) -> Self {
         unsafe { _mm_sub_ps(self, other) }
@@ -656,6 +651,12 @@ impl FloatVector3<f32> for __m128 {
         unsafe { _mm_cmpunord_ps(self, self) }
     }
 
+    #[inline(always)]
+    #[cfg(target_feature = "fma")]
+    fn mul_add(self, b: Self, c: Self) -> Self {
+        unsafe { _mm_fmadd_ps(self, b, c) }
+    }
+
     #[inline]
     fn floor(self) -> Self {
         unsafe { m128_floor(self) }
@@ -749,6 +750,12 @@ impl FloatVector4<f32> for __m128 {
         unsafe { _mm_cmpunord_ps(self, self) }
     }
 
+    #[inline(always)]
+    #[cfg(target_feature = "fma")]
+    fn mul_add(self, b: Self, c: Self) -> Self {
+        unsafe { _mm_fmadd_ps(self, b, c) }
+    }
+
     #[inline]
     fn floor(self) -> Self {
         unsafe { m128_floor(self) }
diff --git a/src/core/traits/matrix.rs b/src/core/traits/matrix.rs
index b251c8c..f1b00eb 100644
--- a/src/core/traits/matrix.rs
+++ b/src/core/traits/matrix.rs
@@ -297,8 +297,8 @@ pub trait Matrix3x3<T: NumEx, V3: Vector3<T>>: Matrix<T> {
     #[inline]
     fn mul_vector(&self, other: V3) -> V3 {
         let mut res = self.x_axis().mul(other.splat_x());
-        res = self.y_axis().mul_add(other.splat_y(), res);
-        res = self.z_axis().mul_add(other.splat_z(), res);
+        res = res.add(self.y_axis().mul(other.splat_y()));
+        res = res.add(self.z_axis().mul(other.splat_z()));
         res
     }
 
@@ -641,9 +641,9 @@ pub trait Matrix4x4<T: NumEx, V4: Vector4<T>>: Matrix<T> {
     #[inline]
     fn mul_vector(&self, other: V4) -> V4 {
         let mut res = self.x_axis().mul(other.splat_x());
-        res = self.y_axis().mul_add(other.splat_y(), res);
-        res = self.z_axis().mul_add(other.splat_z(), res);
-        res = self.w_axis().mul_add(other.splat_w(), res);
+        res = res.add(self.y_axis().mul(other.splat_y()));
+        res = res.add(self.z_axis().mul(other.splat_z()));
+        res = res.add(self.w_axis().mul(other.splat_w()));
         res
     }
 
diff --git a/src/core/traits/scalar.rs b/src/core/traits/scalar.rs
index 1f1f954..bf2b135 100644
--- a/src/core/traits/scalar.rs
+++ b/src/core/traits/scalar.rs
@@ -29,6 +29,7 @@ pub trait Float: Num + Copy + core::ops::Neg<Output = Self> {
     fn floor(self) -> Self;
     fn is_finite(self) -> bool;
     fn is_nan(self) -> bool;
+    fn mul_add(self, b: Self, c: Self) -> Self;
     fn powf(self, n: Self) -> Self;
     fn recip(self) -> Self;
     fn round(self) -> Self;
@@ -98,6 +99,10 @@ macro_rules! impl_float_trait {
                 $t::is_nan(self)
             }
             #[inline(always)]
+            fn mul_add(self, b: Self, c: Self) -> Self {
+                $t::mul_add(self, b, c)
+            }
+            #[inline(always)]
             fn powf(self, n: Self) -> Self {
                 $t::powf(self, n)
             }
diff --git a/src/core/traits/vector.rs b/src/core/traits/vector.rs
index ac67bbc..16ff5fe 100644
--- a/src/core/traits/vector.rs
+++ b/src/core/traits/vector.rs
@@ -92,7 +92,6 @@ pub trait Vector<T>: Sized + Copy + Clone {
     fn add(self, other: Self) -> Self;
     fn div(self, other: Self) -> Self;
     fn mul(self, other: Self) -> Self;
-    fn mul_add(self, a: Self, b: Self) -> Self;
     fn rem(self, rhs: Self) -> Self;
     fn sub(self, other: Self) -> Self;
 
@@ -531,6 +530,14 @@ pub trait FloatVector2<T: FloatEx>: SignedVector2<T> {
         self.x().is_nan() || self.y().is_nan()
     }
 
+    #[inline]
+    fn mul_add(self, b: Self, c: Self) -> Self {
+        Self::new(
+            self.x().mul_add(b.x(), c.x()),
+            self.y().mul_add(b.y(), c.y()),
+        )
+    }
+
     #[inline]
     fn is_nan_mask(self) -> Self::Mask
     where
@@ -630,6 +637,15 @@ pub trait FloatVector3<T: FloatEx>: SignedVector3<T> {
         self.x().is_nan() || self.y().is_nan() || self.z().is_nan()
     }
 
+    #[inline]
+    fn mul_add(self, b: Self, c: Self) -> Self {
+        Self::new(
+            self.x().mul_add(b.x(), c.x()),
+            self.y().mul_add(b.y(), c.y()),
+            self.z().mul_add(b.z(), c.z()),
+        )
+    }
+
     #[inline]
     fn is_nan_mask(self) -> Self::Mask
     where
@@ -753,6 +769,16 @@ pub trait FloatVector4<T: FloatEx>: SignedVector4<T> {
         self.x().is_nan() || self.y().is_nan() || self.z().is_nan() || self.w().is_nan()
     }
 
+    #[inline]
+    fn mul_add(self, b: Self, c: Self) -> Self {
+        Self::new(
+            self.x().mul_add(b.x(), c.x()),
+            self.y().mul_add(b.y(), c.y()),
+            self.z().mul_add(b.z(), c.z()),
+            self.w().mul_add(b.w(), c.w()),
+        )
+    }
+
     #[inline]
     fn is_nan_mask(self) -> Self::Mask
     where
diff --git a/src/core/wasm32/matrix.rs b/src/core/wasm32/matrix.rs
index eca5602..1797a1d 100644
--- a/src/core/wasm32/matrix.rs
+++ b/src/core/wasm32/matrix.rs
@@ -482,8 +482,8 @@ impl FloatMatrix4x4<f32, v128> for Columns4<v128> {
     #[inline]
     fn transform_float4_as_point3(&self, other: v128) -> v128 {
         let mut res = self.x_axis.mul(Vector4::splat_x(other));
-        res = self.y_axis.mul_add(Vector4::splat_y(other), res);
-        res = self.z_axis.mul_add(Vector4::splat_z(other), res);
+        res = res.add(self.y_axis.mul(Vector4::splat_y(other)));
+        res = res.add(self.z_axis.mul(Vector4::splat_z(other)));
         res = self.w_axis.add(res);
         res
     }
@@ -491,16 +491,16 @@ impl FloatMatrix4x4<f32, v128> for Columns4<v128> {
     #[inline]
     fn transform_float4_as_vector3(&self, other: v128) -> v128 {
         let mut res = self.x_axis.mul(Vector4::splat_x(other));
-        res = self.y_axis.mul_add(Vector4::splat_y(other), res);
-        res = self.z_axis.mul_add(Vector4::splat_z(other), res);
+        res = res.add(self.y_axis.mul(Vector4::splat_y(other)));
+        res = res.add(self.z_axis.mul(Vector4::splat_z(other)));
         res
     }
 
     #[inline]
     fn project_float4_as_point3(&self, other: v128) -> v128 {
         let mut res = self.x_axis.mul(Vector4::splat_x(other));
-        res = self.y_axis.mul_add(Vector4::splat_y(other), res);
-        res = self.z_axis.mul_add(Vector4::splat_z(other), res);
+        res = res.add(self.y_axis.mul(Vector4::splat_y(other)));
+        res = res.add(self.z_axis.mul(Vector4::splat_z(other)));
         res = self.w_axis.add(res);
         res = res.mul(res.splat_w().recip());
         res
diff --git a/src/core/wasm32/vector.rs b/src/core/wasm32/vector.rs
index b13ace2..9ecaa77 100644
--- a/src/core/wasm32/vector.rs
+++ b/src/core/wasm32/vector.rs
@@ -223,11 +223,6 @@ impl Vector<f32> for v128 {
         f32x4_mul(self, other)
     }
 
-    #[inline(always)]
-    fn mul_add(self, b: Self, c: Self) -> Self {
-        f32x4_add(f32x4_mul(self, b), c)
-    }
-
     #[inline(always)]
     fn sub(self, other: Self) -> Self {
         f32x4_sub(self, other)
diff --git a/src/quat.rs b/src/quat.rs
index 0fb4abe..135c942 100644
--- a/src/quat.rs
+++ b/src/quat.rs
@@ -282,6 +282,12 @@ macro_rules! impl_quat_methods {
             euler.convert_quat(self)
         }
 
+        /// `[x, y, z, w]`
+        #[inline(always)]
+        pub fn to_array(&self) -> [$t; 4] {
+            [self.x, self.y, self.z, self.w]
+        }
+
         /// Returns the vector part of the quaternion.
         #[inline(always)]
         pub fn xyz(self) -> $vec3 {
diff --git a/src/vec.rs b/src/vec.rs
index 6064d16..3c5a226 100644
--- a/src/vec.rs
+++ b/src/vec.rs
@@ -148,13 +148,6 @@ macro_rules! impl_vecn_common_methods {
         pub fn write_to_slice(self, slice: &mut [$t]) {
             $vectrait::write_to_slice_unaligned(self.0, slice)
         }
-
-        /// Per element multiplication/addition of the three inputs: b + (self * a)
-        #[inline(always)]
-        #[allow(dead_code)]
-        pub(crate) fn mul_add(self, a: Self, b: Self) -> Self {
-            Self(self.0.mul_add(a.0, b.0))
-        }
     };
 }
 
@@ -418,7 +411,8 @@ macro_rules! impl_vecn_float_methods {
         /// Performs a linear interpolation between `self` and `other` based on the value `s`.
         ///
         /// When `s` is `0.0`, the result will be equal to `self`.  When `s` is `1.0`, the result
-        /// will be equal to `other`.
+        /// will be equal to `other`. When `s` is outside of range [0,1], the result is linearly
+        /// extrapolated.
         #[doc(alias = "mix")]
         #[inline]
         pub fn lerp(self, other: Self, s: $t) -> Self {
@@ -476,6 +470,18 @@ macro_rules! impl_vecn_float_methods {
                 self
             }
         }
+
+        /// Fused multiply-add. Computes `(self * a) + b` element-wise with only one rounding
+        /// error, yielding a more accurate result than an unfused multiply-add.
+        ///
+        /// Using `mul_add` *may* be more performant than an unfused multiply-add if the target
+        /// architecture has a dedicated fma CPU instruction. However, this is not always true,
+        /// and will be heavily dependant on designing algorithms with specific target hardware in
+        /// mind.
+        #[inline(always)]
+        pub fn mul_add(self, a: Self, b: Self) -> Self {
+            Self($flttrait::mul_add(self.0, a.0, b.0))
+        }
     };
 }
 
diff --git a/src/vec3.rs b/src/vec3.rs
index a86c481..941aea6 100644
--- a/src/vec3.rs
+++ b/src/vec3.rs
@@ -417,7 +417,7 @@ impl_vecn_scalar_bit_op_traits!(UVec3, u32, XYZU32);
 impl_vecn_bit_op_traits!(UVec3, XYZU32);
 
 #[test]
-fn test_vec3_private() {
+fn test_vec3_mul_add() {
     assert_eq!(
         vec3a(1.0, 1.0, 1.0).mul_add(vec3a(0.5, 2.0, -4.0), vec3a(-1.0, -1.0, -1.0)),
         vec3a(-0.5, 1.0, -5.0)
diff --git a/src/vec4.rs b/src/vec4.rs
index 4ce6f19..2a0b50c 100644
--- a/src/vec4.rs
+++ b/src/vec4.rs
@@ -386,18 +386,19 @@ impl_vecn_scalar_bit_op_traits!(UVec4, u32, XYZWU32);
 
 impl_vecn_bit_op_traits!(UVec4, XYZWU32);
 
-#[test]
-fn test_vec4_private() {
-    assert_eq!(
-        vec4(1.0, 1.0, 1.0, 1.0).mul_add(vec4(0.5, 2.0, -4.0, 0.0), vec4(-1.0, -1.0, -1.0, -1.0)),
-        vec4(-0.5, 1.0, -5.0, -1.0)
-    );
-}
-
 #[cfg(test)]
 mod tests {
     use super::{vec4, Vec3};
 
+    #[test]
+    fn test_vec4_mul_add() {
+        assert_eq!(
+            vec4(1.0, 1.0, 1.0, 1.0)
+                .mul_add(vec4(0.5, 2.0, -4.0, 0.0), vec4(-1.0, -1.0, -1.0, -1.0)),
+            vec4(-0.5, 1.0, -5.0, -1.0)
+        );
+    }
+
     #[test]
     fn from_vec3() {
         assert_eq!(


Test patch:
diff --git a/tests/quat.rs b/tests/quat.rs
index 7ec942c..831bf0f 100644
--- a/tests/quat.rs
+++ b/tests/quat.rs
@@ -461,6 +461,10 @@ macro_rules! impl_quat_tests {
             should_glam_assert!({ $quat::from_rotation_arc_2d($vec2::ZERO, $vec2::X) });
             should_glam_assert!({ $quat::from_rotation_arc_2d($vec2::X, $vec2::ZERO) });
         });
+
+        glam_test!(test_to_array, {
+            assert!($new(1.0, 2.0, 3.0, 4.0).to_array() == [1.0, 2.0, 3.0, 4.0]);
+        });
     };
 }