Move Cython random number generation to C++ mechanism

brian2/codegen/generators/cython_generator.py

@@ -555,30 +555,23 @@ def determine_keywords(self):
     name="_exprel",
     availability_check=C99Check("exprel"),
 )
-_BUFFER_SIZE = 20000
+
+# ==============================================================================
+# Random Number Generation
+# ==============================================================================
+# We use a pre-compiled global RandomGenerator via cythonrng module.
+# This ensures:
+# 1. All code objects share the same RNG (like C++ standalone)
+# 2. Proper state save/restore is possible
+# 3. Identical sequences between Cython and C++ standalone backends
 
 rand_code = """
-cdef double _rand(int _idx):
-    cdef double **buffer_pointer = <double**>_namespace_rand_buffer
-    cdef double *buffer = buffer_pointer[0]
-    cdef _numpy.ndarray _new_rand
-
-    if(_namespace_rand_buffer_index[0] == 0):
-        if buffer != NULL:
-            free(buffer)
-        _new_rand = _numpy.random.rand(_BUFFER_SIZE)
-        buffer = <double *>_numpy.PyArray_DATA(_new_rand)
-        PyArray_CLEARFLAGS(<_numpy.PyArrayObject*>_new_rand, _numpy.NPY_ARRAY_OWNDATA)
-        buffer_pointer[0] = buffer
-
-    cdef double val = buffer[_namespace_rand_buffer_index[0]]
-    _namespace_rand_buffer_index[0] += 1
-    if _namespace_rand_buffer_index[0] == _BUFFER_SIZE:
-        _namespace_rand_buffer_index[0] = 0
-    return val
-""".replace("_BUFFER_SIZE", str(_BUFFER_SIZE))
-
-randn_code = rand_code.replace("rand", "randn").replace("randnom", "random")
+# Import the global random number generator functions
+from brian2.random.cythonrng cimport _rand, _randn
+"""
+
+# randn is included in the same import
+randn_code = rand_code
 
 poisson_code = """
 cdef double _loggam(double x):
@@ -665,20 +658,14 @@ def determine_keywords(self):
     CythonCodeGenerator,
     code=rand_code,
     name="_rand",
-    namespace={
-        "_rand_buffer": device.rand_buffer,
-        "_rand_buffer_index": device.rand_buffer_index,
-    },
+    namespace={},
 )
 
 DEFAULT_FUNCTIONS["randn"].implementations.add_implementation(
     CythonCodeGenerator,
     code=randn_code,
     name="_randn",
-    namespace={
-        "_randn_buffer": device.randn_buffer,
-        "_randn_buffer_index": device.randn_buffer_index,
-    },
+    namespace={},
 )
 DEFAULT_FUNCTIONS["poisson"].implementations.add_implementation(
     CythonCodeGenerator,

brian2/codegen/runtime/cython_rt/templates/common_group.pyx

@@ -69,6 +69,7 @@ cdef extern from "dynamic_array.h":
         size_t rows()
         size_t cols()
         size_t stride()
+
 {% endmacro %}
 
 {% macro before_run() %}

brian2/devices/cpp_standalone/brianlib/randomgenerator.h

@@ -0,0 +1,157 @@
+#ifndef _BRIAN_RANDOMGENERATOR_H
+#define _BRIAN_RANDOMGENERATOR_H
+
+#include <ostream>
+#include <random>
+#include <cmath>
+#include <iostream>
+#include <sstream>
+#include <string>
+
+/**
+ * @brief Random number generator class that provides reproducible
+ * random sequences across different Brian2 backends.
+ *
+ * Uses std::mt19937 (Mersenne Twister) with the Jean-Sebastien Roy
+ * algorithm for converting to uniform doubles, ensuring identical
+ * sequences when seeded with the same value.
+ *
+ * This class is used by both C++ standalone mode and Cython runtime mode
+ * to ensure cross-backend reproducibility.
+ */
+class RandomGenerator
+{
+private:
+    std::mt19937 gen;
+    double stored_gauss;
+    bool has_stored_gauss;
+
+public:
+    RandomGenerator() : has_stored_gauss(false)
+    {
+        seed();
+    }
+
+    /**
+     * @brief Seed with a random value from the system.
+     */
+    void seed()
+    {
+        std::random_device rd;
+        gen.seed(rd());
+        has_stored_gauss = false;
+    }
+
+    /**
+     * @brief Seed with a specific value for reproducibility.
+     * @param seed The seed value.
+     */
+    void seed(unsigned long seed_value)
+    {
+        gen.seed(seed_value);
+        has_stored_gauss = false;
+    }
+
+    /**
+     * @brief Generate a uniform random double in [0, 1)
+     *
+     * Uses the Jean-Sebastien Roy algorithm to extract 53 bits
+     * of randomness from two consecutive MT19937 outputs.
+     * This ensures reproducibility with older Brian2 versions
+     * and across different backends.
+     *
+     * The algorithm:
+     * - Takes two 32-bit outputs from MT19937
+     * - Extracts 27 bits from the first (shift right 5)
+     * - Extracts 26 bits from the second (shift right 6)
+     * - Combines them into a 53-bit integer (full double precision mantissa)
+     * - Divides by 2^53 to get a value in [0, 1)
+     */
+    double rand()
+    {
+        /* shifts : 67108864 = 0x4000000 = 2^26
+         *          9007199254740992 = 0x20000000000000 = 2^53 */
+        const long a = gen() >> 5; // Upper 27 bits
+        const long b = gen() >> 6; // Upper 26 bits
+        return (a * 67108864.0 + b) / 9007199254740992.0;
+    }
+
+    /**
+     * @brief Generate a standard normal random double (mean=0, std=1)
+     *
+     * Uses the Box-Muller transform with rejection sampling.
+     * Generates two values at once and caches one for the next call,
+     * making every other call essentially free.
+     */
+    double randn()
+    {
+        if (has_stored_gauss)
+        {
+            const double tmp = stored_gauss;
+            has_stored_gauss = false;
+            return tmp;
+        }
+        else
+        {
+            double f, x1, x2, r2;
+
+            do
+            {
+                x1 = 2.0 * rand() - 1.0;
+                x2 = 2.0 * rand() - 1.0;
+                r2 = x1 * x1 + x2 * x2;
+            } while (r2 >= 1.0 || r2 == 0.0);
+
+            /* Box-Muller transform */
+            f = sqrt(-2.0 * log(r2) / r2);
+            /* Keep for next call */
+            stored_gauss = f * x1;
+            has_stored_gauss = true;
+            return f * x2;
+        }
+    }
+
+    /**
+     * @brief Get the complete internal state as a string
+     *
+     * This captures the MT19937 state plus the Box-Muller cache,
+     * allowing exact restoration of the generator state.
+     */
+     std::string get_state() const {
+         std::ostringstream oss;
+         oss << gen << " " << stored_gauss << " " << has_stored_gauss;
+         return oss.str();
+     }
+
+     /**
+      * @brief Set the complete internal state from a string.
+      *
+      * Restores the exact state previously captured by get_state().
+      */
+      void set_state(const std::string& state){
+          std::istringstream iss(state);
+          iss >> gen >> stored_gauss >> has_stored_gauss;
+      }
+
+    // Allow exporting/setting the internal state of the random generator
+    friend std::ostream &operator<<(std::ostream &out, const RandomGenerator &rng);
+    friend std::istream &operator>>(std::istream &in, RandomGenerator &rng);
+};
+
+/**
+ * @brief Stream output operator for serializing generator state.
+ */
+inline std::ostream &operator<<(std::ostream &out, const RandomGenerator &rng)
+{
+    return out << rng.gen << " " << rng.stored_gauss << " " << rng.has_stored_gauss;
+}
+
+/**
+ * @brief Stream input operator for deserializing generator state.
+ */
+inline std::istream &operator>>(std::istream &in, RandomGenerator &rng)
+{
+    return in >> rng.gen >> rng.stored_gauss >> rng.has_stored_gauss;
+}
+
+#endif // _BRIAN_RANDOMGENERATOR_H

brian2/devices/cpp_standalone/templates/objects.cpp

@@ -19,6 +19,7 @@ set_variable_from_value(name, {{array_name}}, var_size, (char)atoi(s_value.c_str
 #include "brianlib/clocks.h"
 #include "brianlib/dynamic_array.h"
 #include "brianlib/stdint_compat.h"
+#include "brianlib/randomgenerator.h"
 #include "network.h"
 #include<chrono>
 #include<random>
@@ -37,16 +38,6 @@ std::string results_dir = "results/";  // can be overwritten by --results_dir co
 // For multhreading, we need one generator for each thread.
 std::vector< RandomGenerator > _random_generators;
 
-std::ostream& operator<<(std::ostream& out, const RandomGenerator& rng)
-{
-    return out << rng.gen;
-}
-
-std::istream& operator>>(std::istream& in, RandomGenerator& rng)
-{
-    return in >> rng.gen;
-}
-
 //////////////// networks /////////////////
 {% for net in networks | sort(attribute='name') %}
 Network {{net.name}};
@@ -413,6 +404,7 @@ void _dealloc_arrays()
 #include "brianlib/clocks.h"
 #include "brianlib/dynamic_array.h"
 #include "brianlib/stdint_compat.h"
+#include "brianlib/randomgenerator.h"
 #include "network.h"
 #include<chrono>
 #include<random>
@@ -423,64 +415,6 @@ namespace brian {
 
 extern std::string results_dir;
 
-class RandomGenerator {
-    private:
-        std::mt19937 gen;
-        double stored_gauss;
-        bool has_stored_gauss = false;
-    public:
-        RandomGenerator() {
-            seed();
-        }
-        void seed() {
-            std::random_device rd;
-            gen.seed(rd());
-            has_stored_gauss = false;
-        }
-        void seed(unsigned long seed) {
-            gen.seed(seed);
-            has_stored_gauss = false;
-        }
-        // Allow exporting/setting the internal state of the random generator
-        friend std::ostream& operator<<(std::ostream& out, const RandomGenerator& rng);
-        friend std::istream& operator>>(std::istream& in, RandomGenerator& rng);
-
-        double rand() {
-            /* shifts : 67108864 = 0x4000000, 9007199254740992 = 0x20000000000000 */
-            const long a = gen() >> 5;
-            const long b = gen() >> 6;
-            return (a * 67108864.0 + b) / 9007199254740992.0;
-        }
-
-        double randn() {
-            if (has_stored_gauss) {
-                const double tmp = stored_gauss;
-                has_stored_gauss = false;
-                return tmp;
-            }
-            else {
-                double f, x1, x2, r2;
-
-                do {
-                    x1 = 2.0*rand() - 1.0;
-                    x2 = 2.0*rand() - 1.0;
-                    r2 = x1*x1 + x2*x2;
-                }
-                while (r2 >= 1.0 || r2 == 0.0);
-
-                /* Box-Muller transform */
-                f = sqrt(-2.0*log(r2)/r2);
-                /* Keep for next call */
-                stored_gauss = f*x1;
-                has_stored_gauss = true;
-                return f*x2;
-            }
-        }
-};
-
-extern std::ostream& operator<<(std::ostream& out, const RandomGenerator& rng);
-extern std::istream& operator>>(std::istream& in, RandomGenerator& rng);
-
 // In OpenMP we need one state per thread
 extern std::vector< RandomGenerator > _random_generators;