From: Joshua Simmons <josh@nega.tv>
Date: Sun, 5 May 2024 14:07:04 +0000 (+0200)
Subject: narcissus-maths: Add additional rng functions
X-Git-Url: https://git.nega.tv//gitweb.cgi?a=commitdiff_plain;h=2976afd4220991901faabbbcda72a9f572579d34;p=josh%2Fnarcissus

narcissus-maths: Add additional rng functions

Add `next_[f32,f64]_s`, returning uniformly distibuted floats in the
range `-1.0..1.0`.
Add `next_uniform_unit_disc_f32`, generating a uniformly distributed
point on the unit disc.
Add `next_uniform_unit_circle_f32`, generating a uniformly distributed
point on the unit circle.
---

diff --git a/engine/narcissus-core/src/rand.rs b/engine/narcissus-core/src/rand.rs
index 1ad9147..ff52513 100644
--- a/engine/narcissus-core/src/rand.rs
+++ b/engine/narcissus-core/src/rand.rs
@@ -68,8 +68,8 @@ impl Pcg64 {
     #[inline]
     #[must_use]
     pub fn next_f32(&mut self) -> f32 {
-        let value = (self.next_u64() >> (64 - 24)) as i64 as f32;
-        value * 5.960_464_5e-8 // 0x1p-24f
+        let value = (self.next_u64() >> (64 - 24)) as f32;
+        value * 5.960_464_5e-8 // 0x1p-24
     }
 
     /// Generates a uniformly distributed random f64 in the range `0.0..1.0`
@@ -78,10 +78,71 @@ impl Pcg64 {
     #[inline]
     #[must_use]
     pub fn next_f64(&mut self) -> f64 {
-        let value = (self.next_u64() >> (64 - 53)) as i64 as f64;
+        let value = (self.next_u64() >> (64 - 53)) as f64;
         value * 1.110_223_024_625_156_5e-16 // 0x1p-53
     }
 
+    /// Generates a uniformly distributed random f32 in the range `-1.0..1.0`
+    ///
+    /// Always draws one 64 bit word from the PRNG.
+    #[inline]
+    #[must_use]
+    pub fn next_f32_s(&mut self) -> f32 {
+        let value = (self.next_u64() as i64 >> (64 - 25)) as f32;
+        value * 5.960_464_5e-8 // 0x1p-24
+    }
+
+    /// Generates a uniformly distributed random f64 in the range `-1.0..1.0`
+    ///
+    /// Always draws one 64 bit word from the PRNG.
+    #[inline]
+    #[must_use]
+    pub fn next_f64_s(&mut self) -> f64 {
+        let value = (self.next_u64() as i64 >> (64 - 54)) as f64;
+        value * 1.110_223_024_625_156_5e-16 // 0x1p-53
+    }
+
+    /// Generate a uniformly distributed point on the unit disc using rejection
+    /// sampling.
+    ///
+    /// Returns a tuple containing the dot product of the point with itself, as
+    /// well as the point.
+    ///
+    /// (p.p, [px, py])
+    ///
+    /// # Notes
+    ///
+    /// Uniform point on unit disc by Marc B. Reynolds:
+    /// https://marc-b-reynolds.github.io/distribution/2016/11/28/Uniform.html
+    pub fn next_uniform_unit_disc_f32(&mut self) -> (f32, [f32; 2]) {
+        let mut x;
+        let mut y;
+        let mut d;
+        loop {
+            x = self.next_f32_s();
+            y = self.next_f32_s();
+            d = x * x + y * y;
+            if d < 1.0 {
+                break;
+            }
+        }
+        (d, [x, y])
+    }
+
+    /// Generate a uniformly distributed point on the unit circle.
+    ///
+    /// # Notes
+    ///
+    /// Uniform point on unit circle by Marc B. Reynolds:
+    /// https://marc-b-reynolds.github.io/distribution/2016/11/28/Uniform.html
+    pub fn next_uniform_unit_circle_f32(&mut self) -> [f32; 2] {
+        const BIAS: f32 = 1.0 / 68719476736.0; // 0x1p-36
+        let (d, [x, y]) = self.next_uniform_unit_disc_f32();
+        let s = (d + BIAS * BIAS).sqrt().recip();
+        let x = x + BIAS;
+        [x * s, y * s]
+    }
+
     /// Randomly select an an element from `slice` with uniform probability.
     ///
     /// Always draws two 64 bit words from the PRNG.
@@ -247,4 +308,22 @@ mod tests {
             assert!(x >= 0.0 && x < 1.0);
         }
     }
+
+    #[test]
+    fn uniform_f32_s() {
+        let mut rng = Pcg64::new();
+        for _ in 0..100_000 {
+            let x = rng.next_f32_s();
+            assert!(x >= -1.0 && x < 1.0);
+        }
+    }
+
+    #[test]
+    fn uniform_f64_s() {
+        let mut rng = Pcg64::new();
+        for _ in 0..100_000 {
+            let x = rng.next_f64_s();
+            assert!(x >= -1.0 && x < 1.0);
+        }
+    }
 }