Use a platform-independent and thread-safe RNG

[wegank]

Motivation

msolve used rand()/srand() for all its random choices: the primes of the multi-modular algorithms, random linear forms, the multipliers of the probabilistic linear algebra, and the start vectors and verification values of FGLM. rand() is implementation-defined. Its sequence differs between C libraries, and RAND_MAX is only 32767 on Windows. So the same --random-seed led to different computations on different operating systems. In addition, rand() was called concurrently from OpenMP threads in the probabilistic linear algebra and in the loop over primes, so multi-threaded runs were not reproducible even on a single platform.

Changes

• src/neogb/mt.c, src/neogb/mt.h: the Mersenne Twister MT19937 from GSL 1.9 (rng/mt.c), the last GSL release licensed under GPL-2.0-or-later. GSL's license and copyright notices are kept, and the (dated) modifications are listed at the top of mt.c:

  • the gsl_rng framework and the 1998/1999 seeding variants were removed;
  • mt_set/mt_get are declared in mt.h, so they are no longer static;
  • the code uses uint32_t instead of unsigned long.

  The generated sequence is identical to GSL's gsl_rng_mt19937.

• src/neogb/rng.c, src/neogb/rng.h: msolve's interface to the generator.

  • Each thread owns its own thread-local generator.
  • msolve_srand()/msolve_rand() replace srand()/rand(). msolve_rand() returns values in [0, 2^31 - 1], the same range as glibc's rand().
  • msolve_rng_t provides local generators.

• Parallel loops: random numbers no longer depend on the thread scheduling.

  • Probabilistic linear algebra (the 9 loops in la_ff_{8,16,32}.c): each block draws its multipliers from a local generator seeded with a base seed plus the block index. The base seed is drawn before the parallel loop.
  • Loop over primes (secondary_modular_steps): each prime reseeds the generator of the executing thread with the base seed plus its index. The generator of the calling thread is saved before the loop and restored after it.

• --random-seed (or time(0) by default) seeds the main thread's generator.

• New test neogb_rng (test/neogb/rng/mt19937.c).

Behaviour

• For a given seed, the random numbers are the same on all platforms.
• For a given seed and number of threads, the random choices are reproducible: they do not depend on the thread scheduling.
• For a given seed, the random choices differ from previous msolve versions, since the generator is different. Runs with different numbers of threads may also make different random choices. The final results (reduced Gröbner bases, rational parametrizations) are canonical and were unchanged in all tests.
• libneogb installs two new headers (rng.h, mt.h) and exports the new msolve_rng_*/msolve_rand* functions as well as mt_set/mt_get.

Testing

• make check passes (69 tests, including the new one).
• neogb_rng checks:
  • GSL's reference value (seed 4357: the 1000th output is 1186927261);
  • the value required by the C++ standard for std::mt19937 (seed 5489: the 10000th output is 4123659995);
  • that seed 0 is replaced by 4357, as in GSL.
• The output is identical to std::mt19937 for 7 seeds × 10,000 draws.
• mt.c and rng.c build without warnings (-Wall -Wextra -pedantic) under gnu99, c11, gnu17 and gnu23.
• Reproducibility with several threads:
  • With nonradical-radicalshape-31.ms -P 1 -l 42, master gave different outputs between runs with 4 and 8 threads. With this branch, every seed gives a single output across 1, 4 and 8 threads (40 seeds, 6 runs each).
  • Outputs were also identical across 12 runs with 1 to 8 threads for eco10-31, eco11-31, kat7-qq, henrion5-qq and cyclic5-qq with -l 42/-l 44.
• Tested on macOS (arm64, LLVM clang). The Windows behaviour (32-bit unsigned long) was emulated and gives identical outputs; it has not been tested natively.

Not addressed

• Probabilistic linear algebra: a block stops when a random combination reduces to zero against the pivots known at that moment, including those found by other threads. The multipliers are now reproducible, but in the rare case where this criterion fails, the failure may still depend on timing.
• Multi-modular loop: how often rational reconstruction is attempted adapts to measured wall-clock time (the nbdoit heuristic in msolve.c). So the number of primes used may vary between runs, even with one thread. Outputs are not affected.

Use of generative AI

This PR was prepared with the assistance of generative AI (Claude Code): parts of the code, the tests and this description were written with its help. The commit carries a corresponding Co-Authored-By trailer.

🤖 Generated with Claude Code

[wegank]

This should be rebased once #370 is merged.