diff --git a/lib/gis/lrand48.c b/lib/gis/lrand48.c index c6259411c02..1491dfe63a7 100644 --- a/lib/gis/lrand48.c +++ b/lib/gis/lrand48.c @@ -3,10 +3,25 @@ * * \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 @@ -14,6 +29,15 @@ #include #include +#if defined(__STDC_VERSION__) && __STDC_VERSION__ >= 201112L && \ + !defined(__STDC_NO_ATOMICS__) +#include +#include +#define LRAND48_ATOMIC 1 +#else +#define LRAND48_ATOMIC 0 +#endif + #include #include @@ -30,6 +54,21 @@ 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; @@ -37,6 +76,8 @@ 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 } /* diff --git a/lib/gis/testsuite/test_lrand48.py b/lib/gis/testsuite/test_lrand48.py new file mode 100644 index 00000000000..f63f30eef9d --- /dev/null +++ b/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()