Move sin() into StdFloat trait.

Author: andy-thomason

crates/core_simd/src/vector/float.rs

@@ -45,20 +45,6 @@ macro_rules! impl_float_vector {
                 Self::splat(1.0) / self
             }
 
-            /// Calculate the sine of the angle
-            #[inline]
-            pub fn sin(self) -> Self {
-                let x = Self::splat(1.0 / (3.14159265358 * 2.0)) * self;
-                let x = x - x.floor() - 0.5;
-                Self::splat(12.268859941019306)
-                    .mul_add(x * x, Self::splat(-41.216241051002875))
-                    .mul_add(x * x, Self::splat(76.58672703334098))
-                    .mul_add(x * x, Self::splat(-81.59746095374902))
-                    .mul_add(x * x, Self::splat(41.34151143437585))
-                    .mul_add(x * x, Self::splat(-6.283184525811273))
-                    * x
-            }
-
             /// Converts each lane from radians to degrees.
             #[inline]
             #[must_use = "method returns a new vector and does not mutate the original value"]

crates/core_simd/tests/ops_macros.rs

@@ -397,14 +397,6 @@ macro_rules! impl_float_tests {
                     )
                 }
 
-                fn sin<const LANES: usize>() {
-                    test_helpers::test_unary_elementwise(
-                        &Vector::<LANES>::sin,
-                        &Scalar::sin,
-                        &|_| true,
-                    )
-                }
-
                 fn to_degrees<const LANES: usize>() {
                     test_helpers::test_unary_elementwise(
                         &Vector::<LANES>::to_degrees,

crates/std_float/src/lib.rs

@@ -115,6 +115,8 @@ pub trait StdFloat: Sealed + Sized {
     /// Returns the floating point's fractional value, with its integer part removed.
     #[must_use = "method returns a new vector and does not mutate the original value"]
     fn fract(self) -> Self;
+
+    fn sin(self) -> Self;
 }
 
 impl<const N: usize> Sealed for Simd<f32, N> where LaneCount<N>: SupportedLaneCount {}
@@ -131,6 +133,23 @@ where
     fn fract(self) -> Self {
         self - self.trunc()
     }
+
+    /// Calculate the sine of the angle
+    #[inline]
+    fn sin(self) -> Self {
+        #[allow(non_snake_case)]
+        let RECIP_2PI = Self::splat(0.15915494);
+
+        let scaled = self * RECIP_2PI;
+        let x = scaled - scaled.round();
+        Self::splat(-12.26885994095919635608)
+            .mul_add(x * x, Self::splat(41.21624105096574396575))
+            .mul_add(x * x, Self::splat(-76.58672703333290836700))
+            .mul_add(x * x, Self::splat(81.59746095374827019356))
+            .mul_add(x * x, Self::splat(-41.34151143437582891705))
+            .mul_add(x * x, Self::splat(6.28318452581127506328))
+            * x
+    }
 }
 
 impl<const N: usize> StdFloat for Simd<f64, N>
@@ -143,6 +162,11 @@ where
     fn fract(self) -> Self {
         self - self.trunc()
     }
+
+    #[inline]
+    fn sin(self) -> Self {
+        self
+    }
 }
 
 #[cfg(test)]
@@ -161,5 +185,77 @@ mod tests {
         let _xfma = x.mul_add(x, x);
         let _xsqrt = x.sqrt();
         let _ = x2.abs() * x2;
+        let _ = x.sin();
+    }
+
+    macro_rules! test_range {
+        (
+                min: $min: expr,
+                max: $max: expr,
+                limit: $limit: expr,
+                scalar_fn: $scalar_fn: expr,
+                vector_fn: $vector_fn: expr,
+                scalar_type: $scalar_type: ty,
+                vector_type: $vector_type: ty,
+            ) => {{
+            const NUM_ITER: usize = 0x10000;
+            let limit = <$vector_type>::splat($limit);
+            let b = (($max) - ($min)) * (1.0 / NUM_ITER as $scalar_type);
+            let a = $min;
+            let sf = $scalar_fn;
+            let vf = $vector_fn;
+            for i in (0..NUM_ITER / 4) {
+                let fi = (i * 4) as $scalar_type;
+                let x = <$vector_type>::from_array([
+                    (fi + 0.0) * b + a,
+                    (fi + 1.0) * b + a,
+                    (fi + 2.0) * b + a,
+                    (fi + 3.0) * b + a,
+                ]);
+                let yref = <$vector_type>::from_array([sf(x[0]), sf(x[1]), sf(x[2]), sf(x[3])]);
+                assert!(((vf(x) - yref).abs().lanes_le(limit)).all());
+            }
+        }};
+    }
+
+    #[test]
+    fn sin_f32() {
+        use core::f32::consts::PI;
+        let ulp = (2.0_f32).powi(-23);
+
+        // In the range +/- pi/4 the input has 1 ulp of error.
+        test_range!(
+            min: -PI/4.0,
+            max: PI/4.0,
+            limit: ulp * 1.0,
+            scalar_fn: |x : f32| x.sin(),
+            vector_fn: |x : f32x4| x.sin(),
+            scalar_type: f32,
+            vector_type: f32x4,
+        );
+
+        // In the range +/- pi/2 the input and output has 2 ulp of error.
+        test_range!(
+            min: -PI/2.0,
+            max: PI/2.0,
+            limit: ulp * 2.0,
+            scalar_fn: |x : f32| x.sin(),
+            vector_fn: |x : f32x4| x.sin(),
+            scalar_type: f32,
+            vector_type: f32x4,
+        );
+
+        // In the range +/- pi the input has 2 ulp of error and the output has 5.
+        // Note that the scalar sin also has this error but the implementation
+        // is different.
+        test_range!(
+            min: -PI,
+            max: PI,
+            limit: ulp * 5.0,
+            scalar_fn: |x : f32| x.sin(),
+            vector_fn: |x : f32x4| x.sin(),
+            scalar_type: f32,
+            vector_type: f32x4,
+        );
     }
 }