initial commit

src/unirand.cpp

@@ -0,0 +1,106 @@
+#include "ReUnirand/unirand.hpp"
+#include <stdexcept>  // For std::runtime_error
+
+namespace ReUnirand {
+
+    // Constructor: initialises the internal array and variables.
+    MarsagliaUniRng::MarsagliaUniRng()
+        : uni_c(0.0f), uni_cd(0.0f), uni_cm(0.0f), uni_ui(0), uni_uj(0)
+    {
+        for (int i = 0; i < LEN_U; ++i) {
+            uni_u[i] = 0.0f;
+        }
+    }
+
+    // Generates a new random number.
+    float MarsagliaUniRng::uni() {
+        // Compute the difference between two array values.
+        float luni = uni_u[uni_ui] - uni_u[uni_uj];
+        if (luni < 0.0f) {
+            luni += 1.0f;
+        }
+        uni_u[uni_ui] = luni;
+
+        // Adjust the indices circularly.
+        if (uni_ui == 0)
+            uni_ui = LEN_U - 1;
+        else
+            --uni_ui;
+
+        if (uni_uj == 0)
+            uni_uj = LEN_U - 1;
+        else
+            --uni_uj;
+
+        // Update the correction value.
+        uni_c -= uni_cd;
+        if (uni_c < 0.0f) {
+            uni_c += uni_cm;
+        }
+
+        luni -= uni_c;
+        if (luni < 0.0f) {
+            luni += 1.0f;
+        }
+        return luni;
+    }
+
+    // Initialises the internal state using four seed values.
+    void MarsagliaUniRng::rstart(int i, int j, int k, int l) {
+        // Populate the uni_u array with computed values.
+        // Note: We fill indices 1 to 97, leaving index 0 as previously initialised.
+        for (int ii = 1; ii <= 97; ++ii) {
+            float s = 0.0f;
+            float t = 0.5f;
+            for (int jj = 1; jj <= 24; ++jj) {
+                // Calculate m as per Marsaglia's algorithm.
+                int m = (((i * j) % 179) * k) % 179;
+                i = j;
+                j = k;
+                k = m;
+                l = (53 * l + 1) % 169;
+                if (((l * m) % 64) >= 32) {
+                    s += t;
+                }
+                t *= 0.5f;
+            }
+            uni_u[ii] = s;
+        }
+        // Set fixed correction parameters.
+        uni_c = 362436.0f / 16777216.0f;
+        uni_cd = 7654321.0f / 16777216.0f;
+        uni_cm = 16777213.0f / 16777216.0f;
+        uni_ui = 97;
+        uni_uj = 33;
+    }
+
+    // Decomposes a single seed into four components and initialises the generator.
+    void MarsagliaUniRng::rinit(int ijkl) {
+        if (ijkl < 0 || ijkl > 900000000) {
+            throw std::runtime_error("rinit: seed out of range");
+        }
+
+        int ij = ijkl / 30082;
+        int kl = ijkl - (30082 * ij);
+        int i = ((ij / 177) % 177) + 2;
+        int j = (ij % 177) + 2;
+        int k = ((kl / 169) % 178) + 1;
+        int l = kl % 169;
+
+        // Validate seed components.
+        if (i <= 0 || i > 178)
+            throw std::runtime_error("rinit: i seed out of range");
+        if (j <= 0 || j > 178)
+            throw std::runtime_error("rinit: j seed out of range");
+        if (k <= 0 || k > 178)
+            throw std::runtime_error("rinit: k seed out of range");
+        if (l < 0 || l > 168)
+            throw std::runtime_error("rinit: l seed out of range");
+        if (i == 1 && j == 1 && k == 1)
+            throw std::runtime_error("rinit: 1 1 1 not allowed for first three seeds");
+
+        // Initialise the random number array.
+        rstart(i, j, k, l);
+    }
+
+} // namespace ReUnirand