fix: improve sqrt implementations for f32 and f64 without std

Handle special cases (NaN, zero, negative, infinity) and use
Newton-Raphson refinement for better accuracy. Update RMS doc tests
for floating-point precision.

Author: roderickvd

dasp_rms/src/lib.rs

@@ -117,10 +117,18 @@ 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]);
-    ///     assert_eq!(rms.next([-1.0]), [1.0]);
+    ///
+    ///     assert_eq!(rms.next([1.0f32]), [0.5f32]);
+    ///
+    ///     let result = rms.next([-1.0f32])[0];
+    ///     assert!((result - 0.7071067811865476).abs() < 0.0001,
+    ///             "Expected ~0.7071067811865476, got {}", result);
+    ///
+    ///     let result = rms.next([1.0f32])[0];
+    ///     assert!((result - 0.8660254037844386).abs() < 0.0001,
+    ///             "Expected ~0.8660254037844386, got {}", result);
+    ///
+    ///     assert_eq!(rms.next([-1.0f32]), [1.0f32]);
     /// }
     /// ```
     #[inline]

dasp_sample/src/ops.rs

@@ -4,11 +4,29 @@ pub mod f32 {
 
     #[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.is_nan() || x < 0.0 {
+            return f32::NAN;
         }
+        if x == 0.0 {
+            return x; // preserves +0.0 and -0.0
+        }
+        if x.is_infinite() {
+            return f32::INFINITY;
+        }
+
+        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")]
     pub fn sqrt(x: f32) -> f32 {
@@ -22,11 +40,29 @@ pub mod f64 {
 
     #[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.is_nan() || x < 0.0 {
+            return f64::NAN;
+        }
+        if x == 0.0 {
+            return x; // preserves +0.0 and -0.0
+        }
+        if x.is_infinite() {
+            return f64::INFINITY;
+        }
+
+        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")]
     pub fn sqrt(x: f64) -> f64 {