fix: improve sqrt implementations for f32 and f64 without std (#192)

Handles special cases (NaN, zero, negative, infinity) and uses
Newton-Raphson refinement for better accuracy.

- Updated RMS doc tests to account for >6 decimal precision.
- Inlined std sqrt functions for f32 and f64 for performance.

Author: roderickvd

dasp_rms/src/lib.rs

@@ -117,9 +117,17 @@ where
     /// fn main() {
     ///     let window = ring_buffer::Fixed::from([[0.0]; 4]);
     ///     let mut rms = Rms::new(window);
+    ///
     ///     assert_eq!(rms.next([1.0]), [0.5]);
-    ///     assert_eq!(rms.next([-1.0]), [0.7071067811865476]);
-    ///     assert_eq!(rms.next([1.0]), [0.8660254037844386]);
+    ///
+    ///     let result = rms.next([-1.0])[0];
+    ///     assert!(f64::abs(result - 0.7071067811865476) < 0.000001,
+    ///             "Expected ~0.7071067811865476, got {}", result);
+    ///
+    ///     let result = rms.next([1.0])[0];
+    ///     assert!(f64::abs(result - 0.8660254037844386) < 0.000001,
+    ///             "Expected ~0.8660254037844386, got {}", result);
+    ///
     ///     assert_eq!(rms.next([-1.0]), [1.0]);
     /// }
     /// ```

dasp_sample/src/ops.rs

@@ -1,34 +1,68 @@
 pub mod f32 {
-    #[allow(unused_imports)]
-    use core;
-
+    /// Newton-Raphson square root implementation for f32.
+    /// Uses bit manipulation for initial guess, then 3 iterations for ~6-7 decimal places.
+    /// Accuracy: ~6-7 decimal places
     #[cfg(not(feature = "std"))]
     pub fn sqrt(x: f32) -> f32 {
-        if x >= 0.0 {
-            f32::from_bits((x.to_bits() + 0x3f80_0000) >> 1)
-        } else {
-            f32::NAN
+        if x < 0.0 {
+            return f32::NAN;
+        }
+        if x == 0.0 {
+            return x; // preserves +0.0 and -0.0
+        }
+
+        // Initial guess from bit manipulation: halve exponent, shift mantissa
+        let bits = x.to_bits();
+        let exp = (bits >> 23) & 0xff;
+        let mant = bits & 0x7fffff;
+
+        let unbiased = exp as i32 - 127;
+        let sqrt_exp = (unbiased / 2 + 127) as u32;
+        let guess_bits = (sqrt_exp << 23) | (mant >> 1);
+        let mut guess = f32::from_bits(guess_bits);
+
+        for _ in 0..3 {
+            guess = 0.5 * (guess + x / guess);
         }
+        guess
     }
     #[cfg(feature = "std")]
+    #[inline]
     pub fn sqrt(x: f32) -> f32 {
         x.sqrt()
     }
 }
 
 pub mod f64 {
-    #[allow(unused_imports)]
-    use core;
-
+    /// Newton-Raphson square root implementation for f64.
+    /// Uses bit manipulation for initial guess, then 4 iterations for ~14-15 decimal places.
+    /// Accuracy: ~14-15 decimal places
     #[cfg(not(feature = "std"))]
     pub fn sqrt(x: f64) -> f64 {
-        if x >= 0.0 {
-            f64::from_bits((x.to_bits() + 0x3f80_0000) >> 1)
-        } else {
-            f64::NAN
+        if x < 0.0 {
+            return f64::NAN;
+        }
+        if x == 0.0 {
+            return x; // preserves +0.0 and -0.0
+        }
+
+        // Initial guess from bit manipulation: halve exponent, shift mantissa
+        let bits = x.to_bits();
+        let exp = (bits >> 52) & 0x7ff;
+        let mant = bits & 0x000f_ffff_ffff_ffff;
+
+        let unbiased = exp as i32 - 1023;
+        let sqrt_exp = (unbiased / 2 + 1023) as u64;
+        let guess_bits = (sqrt_exp << 52) | (mant >> 1);
+        let mut guess = f64::from_bits(guess_bits);
+
+        for _ in 0..4 {
+            guess = 0.5 * (guess + x / guess);
         }
+        guess
     }
     #[cfg(feature = "std")]
+    #[inline]
     pub fn sqrt(x: f64) -> f64 {
         x.sqrt()
     }

dasp_signal/src/rms.rs

@@ -33,10 +33,18 @@ pub trait SignalRms: Signal {
     ///     let signal = signal::from_iter(frames.iter().cloned());
     ///     let ring_buffer = ring_buffer::Fixed::from([[0.0]; 2]);
     ///     let mut rms_signal = signal.rms(ring_buffer);
-    ///     assert_eq!(
-    ///         [rms_signal.next(), rms_signal.next(), rms_signal.next()],
-    ///         [[0.6363961030678927], [0.8514693182963201], [0.7071067811865476]]
-    ///     );
+    ///
+    ///     let result = rms_signal.next()[0];
+    ///     assert!(f64::abs(result - 0.6363961030678927) < 0.000001,
+    ///             "Expected ~0.6363961030678927, got {}", result);
+    ///
+    ///     let result = rms_signal.next()[0];
+    ///     assert!(f64::abs(result - 0.8514693182963201) < 0.000001,
+    ///             "Expected ~0.8514693182963201, got {}", result);
+    ///
+    ///     let result = rms_signal.next()[0];
+    ///     assert!(f64::abs(result - 0.7071067811865476) < 0.000001,
+    ///             "Expected ~0.7071067811865476, got {}", result);
     /// }
     /// ```
     ///