lib: Make random number generation thread safe

lib/gis/lrand48.c

@@ -3,17 +3,41 @@
  *
  * \brief GIS Library - Pseudo-random number generation
  *
- * SPDX-FileCopyrightText: 2014 GRASS Development Team
+ * The generator is the standard drand48 linear congruential generator
+ * X' = (A * X + B) mod 2^48 with A = 0x5DEECE66D and B = 0xB.
+ *
+ * When C11 atomic operations are available, the generator state is
+ * advanced with an atomic compare-and-swap and the generating functions
+ * are thread-safe: the sequence of generated values for a given seed is
+ * the same as in a single-threaded run. Which thread receives which
+ * value depends on scheduling, so results are fully reproducible only
+ * with single-threaded execution. Without C11 atomics (notably MSVC,
+ * which defines __STDC_NO_ATOMICS__), the generator falls back to plain
+ * state updates, so multi-threaded usage is safe only when compiled
+ * with C11 atomics.
+ *
+ * The seeding functions are not thread-safe; see G_srand48().
+ *
+ * SPDX-FileCopyrightText: 2014-2026 GRASS Development Team
  * SPDX-License-Identifier: GPL-2.0-or-later
  *
- * \author Glynn Clements
+ * \authors Glynn Clements, Maris Nartiss (thread safety)
  */
 
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
 #include <errno.h>
 
+#if defined(__STDC_VERSION__) && __STDC_VERSION__ >= 201112L && \
+    !defined(__STDC_NO_ATOMICS__)
+#include <stdatomic.h>
+#include <stdint.h>
+#define LRAND48_ATOMIC 1
+#else
+#define LRAND48_ATOMIC 0
+#endif
+
 #include <grass/gis.h>
 #include <grass/glocale.h>
 
@@ -30,13 +54,30 @@ typedef unsigned short uint16;
 typedef unsigned int uint32;
 typedef signed int int32;
 
+#if LRAND48_ATOMIC
+
+/* The whole 48-bit state is kept in one atomic integer so that it can be
+ * advanced in one compare-and-swap: a successful swap is exactly one
+ * generator step, giving the same sequence of states as a single-threaded
+ * run. The multiplication may wrap around at 2^64; that does not change
+ * the result modulo 2^48. */
+static atomic_uint_least64_t state;
+
+#define LCG_A  UINT64_C(0x5DEECE66D)
+#define LCG_B  UINT64_C(0xB)
+#define MASK48 UINT64_C(0xFFFFFFFFFFFF)
+
+#else
+
 static uint16 x0, x1, x2;
 static const uint32 a0 = 0xE66D;
 static const uint32 a1 = 0xDEEC;
 static const uint32 a2 = 0x5;
 
 static const uint32 b0 = 0xB;
 
+#endif /* LRAND48_ATOMIC */
+
 static int seeded;
 
 #define LO(x) ((x) & 0xFFFFU)
@@ -45,15 +86,23 @@ static int seeded;
 /*!
  * \brief Seed the pseudo-random number generator
  *
+ * This function is not thread-safe. In a multi-threaded program, call
+ * `G_srand48()` once *before* starting the worker threads; it must not
+ * run concurrently with another thread seeding or generating values.
+ *
  * \param[in] seedval 32-bit integer used to seed the PRNG
  */
 void G_srand48(long seedval)
 {
     uint32 x = (uint32) * (unsigned long *)&seedval;
 
+#if LRAND48_ATOMIC
+    atomic_store(&state, ((uint_least64_t)x << 16) | 0x330E);
+#else
     x2 = (uint16)HI(x);
     x1 = (uint16)LO(x);
     x0 = (uint16)0x330E;
+#endif
     seeded = 1;
 }
 
@@ -63,6 +112,10 @@ void G_srand48(long seedval)
  * A weak hash of the current time and PID is generated and used to
  * seed the PRNG
  *
+ * This function is not thread-safe. In a multi-threaded program, call
+ * `G_srand48_auto()` once *before* starting the worker threads; it must
+ * not run concurrently with another thread seeding or generating values.
+ *
  * \return generated seed value passed to G_srand48()
  */
 long G_srand48_auto(void)
@@ -100,6 +153,29 @@ long G_srand48_auto(void)
     return (long)seed;
 }
 
+#if LRAND48_ATOMIC
+
+/* Advance the generator by one step and return the new state. Callers
+ * derive their result from the returned value, not from shared state, so
+ * concurrent calls each get a distinct step of the sequence. */
+static uint_least64_t G__next(void)
+{
+    uint_least64_t cur = atomic_load_explicit(&state, memory_order_relaxed);
+    uint_least64_t next;
+
+    if (!seeded)
+        G_fatal_error(_("Pseudo-random number generator not seeded"));
+
+    do {
+        next = (LCG_A * cur + LCG_B) & MASK48;
+    } while (!atomic_compare_exchange_weak_explicit(
+        &state, &cur, next, memory_order_relaxed, memory_order_relaxed));
+
+    return next;
+}
+
+#else
+
 static void G__next(void)
 {
     uint32 a0x0 = a0 * x0;
@@ -123,41 +199,66 @@ static void G__next(void)
     x2 = (uint16)LO(y2);
 }
 
+#endif /* LRAND48_ATOMIC */
+
 /*!
  * \brief Generate an integer in the range [0, 2^31)
  *
+ * This function is thread-safe only when compiled with C11 atomics
+ * (see the comment at the top of the file).
+ *
  * \return the generated value
  */
 long G_lrand48(void)
 {
+#if LRAND48_ATOMIC
+    return (long)(G__next() >> 17);
+#else
     uint32 r;
 
     G__next();
     r = ((uint32)x2 << 15) | ((uint32)x1 >> 1);
     return (long)r;
+#endif
 }
 
 /*!
  * \brief Generate an integer in the range [-2^31, 2^31)
  *
+ * This function is thread-safe only when compiled with C11 atomics
+ * (see the comment at the top of the file).
+ *
  * \return the generated value
  */
 long G_mrand48(void)
 {
+#if LRAND48_ATOMIC
+    uint32 r = (uint32)(G__next() >> 16);
+
+    return (long)(int32)r;
+#else
     uint32 r;
 
     G__next();
     r = ((uint32)x2 << 16) | ((uint32)x1);
     return (long)(int32)r;
+#endif
 }
 
 /*!
  * \brief Generate a floating-point value in the range [0,1)
  *
+ * This function is thread-safe only when compiled with C11 atomics
+ * (see the comment at the top of the file).
+ *
  * \return the generated value
  */
 double G_drand48(void)
 {
+#if LRAND48_ATOMIC
+    /* The state is below 2^53, so the conversion to double is exact. */
+    return (double)G__next() / 281474976710656.0; /* 2^48 */
+#else
     double r = 0.0;
 
     G__next();
@@ -168,6 +269,7 @@ double G_drand48(void)
     r += x0;
     r /= 281474976710656.0; /* 2^48 */
     return r;
+#endif
 }
 
 /*

lib/gis/testsuite/test_lrand48.py

@@ -0,0 +1,75 @@
+"""Test of gis library lrand48 PRNG thread-safety
+
+@author Maris Nartiss
+
+SPDX-FileCopyrightText: 2025 GRASS Development Team
+SPDX-License-Identifier: GPL-2.0-or-later
+"""
+
+import threading
+
+from grass.gunittest.case import TestCase
+from grass.gunittest.main import test
+from grass.lib.gis import G_lrand48, G_srand48
+
+
+class Lrand48ThreadSafetyTestCase(TestCase):
+    """Test case for lrand48 thread-safety and reproducibility."""
+
+    def test_thread_safety_and_reproducibility(self):
+        """Verify that multi-threaded execution produces the same set of
+        random numbers as single-threaded execution."""
+
+        seed = 1337
+        num_values = 10000
+        num_threads = 4
+        values_per_thread = num_values // num_threads
+
+        list_single = []
+        G_srand48(seed)
+        for _ in range(num_values):
+            list_single.append(G_lrand48())
+
+        list_multi_raw = []
+        lock = threading.Lock()
+
+        def worker():
+            """Calls G_lrand48 and appends the result to a shared list."""
+            local_results = []
+            for _ in range(values_per_thread):
+                local_results.append(G_lrand48())
+
+            # The lock protects only the Python result list, not the C
+            # generator under test (list.extend is atomic under the GIL,
+            # but that is a CPython implementation detail).
+            with lock:
+                list_multi_raw.extend(local_results)
+
+        # Reset the seed to ensure the sequence starts from the beginning.
+        G_srand48(seed)
+
+        threads = []
+        for _ in range(num_threads):
+            thread = threading.Thread(target=worker)
+            threads.append(thread)
+            thread.start()
+
+        for thread in threads:
+            thread.join()
+
+        # The generator serializes state updates, so the threads together
+        # must consume exactly the single-threaded sequence; only the
+        # distribution of values between threads may differ. Sorting both
+        # lists removes the scheduling-dependent order before comparison.
+        # A set-based comparison would not work, because a correct
+        # sequence can contain legitimate duplicates.
+        self.assertListEqual(
+            sorted(list_single),
+            sorted(list_multi_raw),
+            "The set of generated numbers differs between "
+            "single-threaded and multi-threaded runs.",
+        )
+
+
+if __name__ == "__main__":
+    test()