From a38ba26e93701ca8c4cfe6d07023d901b6b27443 Mon Sep 17 00:00:00 2001 From: Weijia Wang <9713184+wegank@users.noreply.github.com> Date: Thu, 1 Oct 2026 19:47:33 +0200 Subject: [PATCH] Use a platform-independent and thread-safe RNG Replace the calls to rand() in the library by a pseudo-random number generator that produces the same sequence on all platforms for a given seed: rand() is implementation-defined (e.g. RAND_MAX is only 32767 on Windows), so the random choices of msolve, such as primes and random linear forms, differed between operating systems. The generator is MT19937, taken from GSL 1.9 (src/neogb/mt.c), the last GSL release licensed under GPL version 2 or later. The imported file keeps GSL's license and copyright notices; the modifications (the GSL framework and the 1998/1999 seeding variants were removed, the code uses fixed-width types) are listed at the top of the file. The sequence coincides with GSL's gsl_rng_mt19937, which is checked by a new test. Each thread owns its own generator, the one of the main thread is seeded via --random-seed. Random numbers drawn inside OpenMP parallel loops (the random multipliers of the probabilistic linear algebra and the computations over several primes) now come from generators seeded from a base seed and the loop index, so that they do not depend on the thread scheduling: for a given seed and number of threads, the random choices are reproducible. The seed of diff_hankel-degenerate-16.sh, which was chosen so that glibc's rand() exercises the degenerate Hankel case, is replaced by one that does so with the new generator. Co-Authored-By: Claude Opus 5.5 --- .gitignore | 1 + Makefile.am | 2 + src/fglm/fglm_core.c | 17 +-- src/msolve/lifting-gb.c | 6 +- src/msolve/main.c | 4 +- src/msolve/msolve.c | 24 ++-- src/neogb/Makefile.am | 6 +- src/neogb/gb.c | 2 + src/neogb/la_ff_16.c | 2 +- src/neogb/la_ff_32.c | 5 +- src/neogb/la_ff_64.c | 11 +- src/neogb/la_ff_8.c | 2 +- src/neogb/la_ff_common.c | 23 ++-- src/neogb/libneogb.h | 1 + src/neogb/mt.c | 153 +++++++++++++++++++++++++ src/neogb/mt.h | 78 +++++++++++++ src/neogb/rng.c | 78 +++++++++++++ src/neogb/rng.h | 75 ++++++++++++ test/diff/diff_hankel-degenerate-16.sh | 4 +- test/neogb/rng/mt19937.c | 42 +++++++ 20 files changed, 497 insertions(+), 39 deletions(-) create mode 100644 src/neogb/mt.c create mode 100644 src/neogb/mt.h create mode 100644 src/neogb/rng.c create mode 100644 src/neogb/rng.h create mode 100644 test/neogb/rng/mt19937.c diff --git a/.gitignore b/.gitignore index bfb42bd3..893b6685 100644 --- a/.gitignore +++ b/.gitignore @@ -52,6 +52,7 @@ fglm_build_matrixn_radical-shape-31 line_endings_support neogb_io neogb_f4_52_bit +neogb_rng msolve-*.tar.gz msolve-*/ src/msolve/msolve-version.h diff --git a/Makefile.am b/Makefile.am index b025b44b..84890384 100644 --- a/Makefile.am +++ b/Makefile.am @@ -9,6 +9,7 @@ msolve_SOURCES = src/msolve/main.c check_PROGRAMS = neogb_io \ neogb_f4_52_bit \ + neogb_rng \ fglm_build_matrixn_radical-shape-31 \ fglm_build_matrixn_nonradical-shape-31 \ fglm_build_matrixn_nonradical-radicalshape-31 \ @@ -87,6 +88,7 @@ checkdiff = test/diff/diff_bug-2nd-prime-bad.sh \ # dist_check_DATA = test/input_files neogb_io_SOURCES = test/neogb/io/validate_input_data.c neogb_f4_52_bit_SOURCES = test/neogb/f4/f4_52_bit.c +neogb_rng_SOURCES = test/neogb/rng/mt19937.c fglm_build_matrixn_radical_shape_31_SOURCES = test/fglm/build_matrixn_radical-shape-31.c fglm_build_matrixn_nonradical_shape_31_SOURCES = test/fglm/build_matrixn_nonradical-shape-31.c fglm_build_matrixn_nonradical_radicalshape_31_SOURCES = test/fglm/build_matrixn_nonradical-radicalshape-31.c diff --git a/src/fglm/fglm_core.c b/src/fglm/fglm_core.c index f11c0fe5..7d58557f 100644 --- a/src/fglm/fglm_core.c +++ b/src/fglm/fglm_core.c @@ -29,6 +29,7 @@ #include /* for timing functions */ #include "../neogb/tools.h" +#include "../neogb/rng.h" #include "../msolve/streams.h" #ifdef _OPENMP @@ -618,8 +619,8 @@ static void generate_matrix_sequence(sp_matfglm_t *matxn, fglm_data_t *data, Rmat[i] = 0; } for(szmat_t i = 0; i < matxn->ncols; i++){ - Rmat[i] = (CF_t)rand() % prime; - Rmat[i] += (CF_t)rand() % prime; + Rmat[i] = (CF_t)msolve_rand() % prime; + Rmat[i] += (CF_t)msolve_rand() % prime; } /* allocates result matrix (matxn * Rmat) */ @@ -1150,7 +1151,7 @@ int compute_parametrizations_non_shape_position_case(param_t *param, if(linvars[nvars - 2 - nc] == 0 && squvars[nvars - 2 - nc] != 0){ - uint64_t lambda= 1 + ((uint64_t) rand() % (prime-1)); + uint64_t lambda= 1 + ((uint64_t) msolve_rand() % (prime-1)); /* needed for verification */ invert_table_polynomial (param, data, data_bms, block_size, @@ -1257,7 +1258,7 @@ static int compute_parametrizations_colon(param_t *param, if(linvars[nvars - 2 - nc] == 0 && squvars[nvars - 2 - nc] != 0){ - uint64_t lambda= 1 + ((uint64_t) rand() % (prime-1)); + uint64_t lambda= 1 + ((uint64_t) msolve_rand() % (prime-1)); /* needed for verification */ invert_table_polynomial (param, data, data_bms, block_size, @@ -1323,9 +1324,9 @@ static inline long initialize_fglm_data(sp_matfglm_t *matrix, if(matrix->dense_mat[i]==0) nb++; } - data->vecinit[0] = 1 + ((CF_t) rand() % (prime-1)); /* random, nonzero */ + data->vecinit[0] = 1 + ((CF_t) msolve_rand() % (prime-1)); /* random, nonzero */ for(szmat_t i = 1; i < matrix->ncols; i++){ - data->vecinit[i] = (CF_t)rand() % prime; + data->vecinit[i] = (CF_t)msolve_rand() % prime; } data->res[0] = data->vecinit[0]; for(szmat_t i = 1; i < block_size; i++){ @@ -1345,8 +1346,8 @@ static inline long initialize_fglm_colon_data(sp_matfglmcol_t *matrix, nb++; } for(szmat_t i = 0; i < matrix->ncols; i++){ - data->vecinit[i] = (CF_t)rand() % prime; - data->vecinit[i] += (CF_t)rand() % prime; + data->vecinit[i] = (CF_t)msolve_rand() % prime; + data->vecinit[i] += (CF_t)msolve_rand() % prime; /* data->vecinit[i] = (CF_t)(i+1) % prime; */ } data->res[0] = data->vecinit[0]; diff --git a/src/msolve/lifting-gb.c b/src/msolve/lifting-gb.c index f529cfb4..e2cc7549 100644 --- a/src/msolve/lifting-gb.c +++ b/src/msolve/lifting-gb.c @@ -1522,9 +1522,9 @@ gb_modpoly_t *core_groebner_qq( uint32_t primeinit = 0; uint32_t lprime = 1303905299; - prime = next_prime(rand() % (1303905301 - (1<<30) + 1) + (1<<30)); + prime = next_prime(msolve_rand() % (1303905301 - (1<<30) + 1) + (1<<30)); while(fc == 0 && is_lucky_prime_ui(prime, bs)){ - prime = next_prime(rand() % (1303905301 - (1<<30) + 1) + (1<<30)); + prime = next_prime(msolve_rand() % (1303905301 - (1<<30) + 1) + (1<<30)); } primeinit = prime; @@ -1640,7 +1640,7 @@ gb_modpoly_t *core_groebner_qq( learn = 0; - prime = next_prime(rand() % (1303905301 - (1<<30) + 1) + (1<<30)); + prime = next_prime(msolve_rand() % (1303905301 - (1<<30) + 1) + (1<<30)); if(info_level){ fprintf(VERBSTREAM, "New prime = %d\n", prime); } diff --git a/src/msolve/main.c b/src/msolve/main.c index c2fa1579..6ee64cb6 100644 --- a/src/msolve/main.c +++ b/src/msolve/main.c @@ -561,14 +561,14 @@ int main(int argc, char **argv){ &precision, &generate_pbm, &seed, &info_level, files); - /* srand initialization */ + /* pseudo-random number generator initialization */ uint32_t true_seed; if (seed < 0) { true_seed = time(0); } else { true_seed = (uint32_t) seed; } - srand(true_seed); + msolve_srand(true_seed); if (info_level) { fprintf (VERBSTREAM,"Initial seed for pseudo-random number generator "); fprintf (VERBSTREAM,"is %u\n",true_seed); diff --git a/src/msolve/msolve.c b/src/msolve/msolve.c index 983344a8..1360dab7 100644 --- a/src/msolve/msolve.c +++ b/src/msolve/msolve.c @@ -756,9 +756,9 @@ static int add_random_linear_form_to_input_system(data_gens_ff_t *gens, if (gens->field_char > 0) { int j = 0; for (i = len_old; i < len_new; i++, j++) { - int32_t rdval = ((int8_t)(rand()) % gens->field_char); + int32_t rdval = ((int8_t)(msolve_rand()) % gens->field_char); while (rdval == 0) { - rdval = ((int8_t)(rand()) % gens->field_char); + rdval = ((int8_t)(msolve_rand()) % gens->field_char); } gens->random_linear_form[j] = rdval; gens->cfs[i] = rdval; @@ -766,9 +766,9 @@ static int add_random_linear_form_to_input_system(data_gens_ff_t *gens, } else { int j = 0; for (i = len_old; i < len_new; i++, j++) { - int32_t rdval = rand(); + int32_t rdval = msolve_rand(); while (rdval == 0) { - rdval = rand(); + rdval = msolve_rand(); } gens->random_linear_form[j] = rdval; mpz_set_si(*(gens->mpz_cfs[2*i]), rdval); @@ -2230,9 +2230,12 @@ static void secondary_modular_steps(sp_matfglm_t **bmatrix, for(nvars_t i = 0; i < st->nprimes; i++){ bad_primes[i] = 0; } + const uint32_t rng_base = msolve_rand_u32(); + const msolve_rng_t rng_saved = *msolve_rng_current(); #pragma omp parallel for num_threads(nthrds) \ private(i) schedule(static) for (i = 0; i < st->nprimes; ++i){ + msolve_srand(msolve_rng_derive_seed(rng_base, i)); if (trace_det->mat_lifted < 2 || trace_det->lin_lifted < 2) { bs[i] = core_gba(bs_qq, st, &error, lp->p[i]); *stf4 = realtime()-rt; @@ -2320,6 +2323,7 @@ static void secondary_modular_steps(sp_matfglm_t **bmatrix, free_basis_and_only_local_hash_table_data(&(bs[i])); } } + *msolve_rng_current() = rng_saved; st->nthrds = nthrds; } @@ -2819,9 +2823,9 @@ int msolve_trace_qq(mpz_param_t *mpz_paramp, uint32_t lprime = 1303905299; /* choose next_prime of some random number between 1<<30 and 1303905301 */ /* (roughly, there are ~ 10**7 such primes) */ - prime = next_prime(rand() % (1303905301 - (1 << 30) + 1) + (1 << 30)); + prime = next_prime(msolve_rand() % (1303905301 - (1 << 30) + 1) + (1 << 30)); while (gens->field_char == 0 && (is_lucky_prime_ui(prime, bs_qq) || is_member(prime, (*binit_primes)))) { - prime = next_prime(rand() % (1303905301 - (1 << 30) + 1) + (1 << 30)); + prime = next_prime(msolve_rand() % (1303905301 - (1 << 30) + 1) + (1 << 30)); } primeinit = prime; lp->p[0] = primeinit; @@ -5696,9 +5700,9 @@ int core_msolve( trace_t *trace = initialize_trace(bs_qq, st); uint32_t prime = next_prime(1<<30); - prime = next_prime(rand() % (1303905301 - (1<<30) + 1) + (1<<30)); + prime = next_prime(msolve_rand() % (1303905301 - (1<<30) + 1) + (1<<30)); while(is_lucky_prime_ui(prime, bs_qq)){ - prime = next_prime(rand() % (1303905301 - (1<<30) + 1) + (1<<30)); + prime = next_prime(msolve_rand() % (1303905301 - (1<<30) + 1) + (1<<30)); } uint32_t primeinit = prime; @@ -5880,9 +5884,9 @@ int core_msolve( st->tr = trace; uint32_t prime = next_prime(1<<30); - prime = next_prime(rand() % (1303905301 - (1<<30) + 1) + (1<<30)); + prime = next_prime(msolve_rand() % (1303905301 - (1<<30) + 1) + (1<<30)); while(is_lucky_prime_ui(prime, bs_qq)){ - prime = next_prime(rand() % (1303905301 - (1<<30) + 1) + (1<<30)); + prime = next_prime(msolve_rand() % (1303905301 - (1<<30) + 1) + (1<<30)); } uint32_t primeinit = prime; diff --git a/src/neogb/Makefile.am b/src/neogb/Makefile.am index 1f3c2a4b..6904102e 100644 --- a/src/neogb/Makefile.am +++ b/src/neogb/Makefile.am @@ -3,7 +3,7 @@ libneogb_la_SOURCES = libneogb.h gb.c libneogb_ladir = $(includedir)/msolve/neogb libneogb_la_HEADERS =libneogb.h basis.h data.h engine.h f4.h sba.h hash.h \ io.h modular.h nf.h f4sat.h sort_r.h meta_data.h \ - tools.h update.h + tools.h update.h rng.h mt.h libneogb_la_LDFLAGS = -version-info $(LT_VERSION) libneogb_la_CFLAGS = $(SIMD_FLAGS) $(CPUEXT_FLAGS) $(OPENMP_CFLAGS) -Wall -Wextra -Wno-unused-parameter -Wno-sign-compare @@ -21,6 +21,8 @@ EXTRA_DIST = basis.h \ meta_data.h \ tools.h \ update.h \ + rng.h \ + mt.h \ basis.c \ convert.c \ data.c \ @@ -43,4 +45,6 @@ EXTRA_DIST = basis.h \ meta_data.c \ symbol.c \ tools.c \ + rng.c \ + mt.c \ update.c diff --git a/src/neogb/gb.c b/src/neogb/gb.c index 6d6b852e..b5a81757 100644 --- a/src/neogb/gb.c +++ b/src/neogb/gb.c @@ -29,6 +29,8 @@ #include "libneogb.h" #include "data.c" #include "meta_data.c"/* computational meta data */ +#include "mt.c" /* MT19937 from GSL 1.9 */ +#include "rng.c" /* platform independent pseudo-random numbers */ #include "tools.c" /* tools like inversion mod p, * tracer construction, timings etc. */ #include "sort_r.h" /* special quicksort implementation */ diff --git a/src/neogb/la_ff_16.c b/src/neogb/la_ff_16.c index 597f63d7..c915320f 100644 --- a/src/neogb/la_ff_16.c +++ b/src/neogb/la_ff_16.c @@ -420,7 +420,7 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_ff_16( #define LA_CF_T cf16_t #define LA_CF cf_16 #define LA_FF(name) name##_ff_16 -#define LA_RANDOM_MULTIPLIER(fc) ((int64_t)rand() % (fc)) +#define LA_RANDOM_MULTIPLIER(rng, fc) ((int64_t)msolve_rng_rand(rng) % (fc)) #define LA_SMALL(name) name##_ff_16 #include "la_ff_small_prime.c" diff --git a/src/neogb/la_ff_32.c b/src/neogb/la_ff_32.c index f34d73c3..a677bcde 100644 --- a/src/neogb/la_ff_32.c +++ b/src/neogb/la_ff_32.c @@ -189,7 +189,7 @@ static int is_kernel_trivial( /* fill random value array */ for (i = 0; i < ncr; ++i) { do { - mull[j] = (int64_t)rand() % fc; + mull[j] = (int64_t)msolve_rand() % fc; } while (mull[j] == 0); } for (i = 0; i < sat->ld; ++i) { @@ -1997,7 +1997,8 @@ static cf32_t *reduce_dense_row_by_dense_new_pivots_31_bit( #define LA_FF(name) name##_ff_32 /* restrict random multipliers to < 2^15, see * probabilistic_sparse_reduced_echelon_form_ff_32() */ -#define LA_RANDOM_MULTIPLIER(fc) ((int64_t)rand() & ((1 << 15) - 1)) +#define LA_RANDOM_MULTIPLIER(rng, fc) \ + ((int64_t)msolve_rng_rand(rng) & ((1 << 15) - 1)) #define LA_SMALL(name) name##_17_bit #include "la_ff_small_prime.c" diff --git a/src/neogb/la_ff_64.c b/src/neogb/la_ff_64.c index 35d86aae..5d874ac8 100644 --- a/src/neogb/la_ff_64.c +++ b/src/neogb/la_ff_64.c @@ -1183,8 +1183,15 @@ static cf64_t *reduce_dense_row_by_all_pivots_ff_64( #define LA_CF_T cf64_t #define LA_CF cf_64 #define LA_FF(name) name##_ff_64 -#define LA_RANDOM_MULTIPLIER(fc) \ - ((int64_t)((((uint64_t)rand() << 31) ^ (uint64_t)rand()) % (fc))) +/* the two draws are sequenced so that the multiplier does not depend on + * the evaluation order chosen by the compiler */ +static inline int64_t random_multiplier_ff_64(msolve_rng_t *rng, uint64_t fc) +{ + const uint64_t hi = (uint64_t)msolve_rng_rand(rng); + const uint64_t lo = (uint64_t)msolve_rng_rand(rng); + return (int64_t)(((hi << 31) ^ lo) % fc); +} +#define LA_RANDOM_MULTIPLIER(rng, fc) random_multiplier_ff_64((rng), (fc)) /* maximal number of rows reduced together by the known pivots, see * LA_FF(reduce_dense_rows_by_known_pivots) in la_ff_common.c */ diff --git a/src/neogb/la_ff_8.c b/src/neogb/la_ff_8.c index ce818392..43dfb323 100644 --- a/src/neogb/la_ff_8.c +++ b/src/neogb/la_ff_8.c @@ -566,7 +566,7 @@ static hm_t *reduce_dense_row_by_known_pivots_sparse_ff_8( #define LA_CF_T cf8_t #define LA_CF cf_8 #define LA_FF(name) name##_ff_8 -#define LA_RANDOM_MULTIPLIER(fc) ((int64_t)rand() % (fc)) +#define LA_RANDOM_MULTIPLIER(rng, fc) ((int64_t)msolve_rng_rand(rng) % (fc)) #define LA_SMALL(name) name##_ff_8 #include "la_ff_small_prime.c" diff --git a/src/neogb/la_ff_common.c b/src/neogb/la_ff_common.c index 919222ec..6732e0a2 100644 --- a/src/neogb/la_ff_common.c +++ b/src/neogb/la_ff_common.c @@ -32,8 +32,8 @@ * mat_t, e.g. cf_8 * LA_FF(name) appends the field suffix to name, e.g. * LA_FF(foo) expands to foo_ff_8 - * LA_RANDOM_MULTIPLIER(fc) random multiplier for the probabilistic - * linear algebra + * LA_RANDOM_MULTIPLIER(rng, fc) random multiplier for the probabilistic + * linear algebra, drawn from the generator rng * * Moreover, the following functions have to be defined: * @@ -635,6 +635,7 @@ static void LA_FF(probabilistic_sparse_reduced_echelon_form)( /* mo need to have any sharing dependencies on parallel computation, * no data to be synchronized at this step of the linear algebra */ + const uint32_t rng_base = msolve_rand_u32(); #pragma omp parallel for num_threads(st->nthrds) \ private(i, j, k, l, m) \ schedule(dynamic) @@ -649,6 +650,8 @@ static void LA_FF(probabilistic_sparse_reduced_echelon_form)( /* starting column, offset, coefficient array position in tmpcf */ hm_t sc, cfp; len_t bctr = 0; + msolve_rng_t rng; + msolve_rng_seed(&rng, msolve_rng_derive_seed(rng_base, i)); /* the combinations start at the smallest lead term of the * rows of the block. As in the exact linear algebra, they are * only reduced together if the rows fill at least 1/32 of @@ -673,7 +676,7 @@ static void LA_FF(probabilistic_sparse_reduced_echelon_form)( /* fill random value array */ for (j = 0; j < nrbl; ++j) { do { - mull[j] = LA_RANDOM_MULTIPLIER(fc); + mull[j] = LA_RANDOM_MULTIPLIER(&rng, fc); } while (mull[j] == 0); } /* generate one dense row as random linear combination @@ -758,7 +761,7 @@ static void LA_FF(probabilistic_sparse_reduced_echelon_form)( for (j = 0; j < nl; ++j) { int64_t c; do { - c = LA_RANDOM_MULTIPLIER(fc); + c = LA_RANDOM_MULTIPLIER(&rng, fc); } while (c == 0); mm[j] = ff64_mont((uint64_t)c, &f); } @@ -772,7 +775,7 @@ static void LA_FF(probabilistic_sparse_reduced_echelon_form)( int64_t mc[LA_NR]; for (j = 0; j < nl; ++j) { do { - mc[j] = LA_RANDOM_MULTIPLIER(fc); + mc[j] = LA_RANDOM_MULTIPLIER(&rng, fc); } while (mc[j] == 0); } for (l = 0; l < len; ++l) { @@ -2065,6 +2068,7 @@ static LA_CF_T **LA_FF(probabilistic_dense_linear_algebra)( (uint64_t)rpb * st->nthrds * sizeof(int64_t)); /* reduction process to get all possible pivots, no interreduction here */ + const uint32_t rng_base = msolve_rand_u32(); #pragma omp parallel for num_threads(st->nthrds) \ private(i, j, k, l) shared(nps, tbr) \ schedule(dynamic) @@ -2079,6 +2083,8 @@ static LA_CF_T **LA_FF(probabilistic_dense_linear_algebra)( hm_t os; LA_CF_T *tmp; len_t bctr = 0; + msolve_rng_t rng; + msolve_rng_seed(&rng, msolve_rng_derive_seed(rng_base, i)); while (bctr < nrbl) { npc = 0; #ifndef LA_WIDE @@ -2088,7 +2094,7 @@ static LA_CF_T **LA_FF(probabilistic_dense_linear_algebra)( /* fill random value array */ for (j = 0; j < nrbl; ++j) { do { - mull[j] = LA_RANDOM_MULTIPLIER(fc); + mull[j] = LA_RANDOM_MULTIPLIER(&rng, fc); } while (mull[j] == 0); } /* generate one dense row as random linear combination @@ -2214,6 +2220,7 @@ static LA_CF_T **LA_FF(probabilistic_sparse_dense_echelon_form)( (uint64_t)rpb * st->nthrds * sizeof(int64_t)); /* reduction process to get all possible pivots, no interreduction here */ + const uint32_t rng_base = msolve_rand_u32(); #pragma omp parallel for num_threads(st->nthrds) \ private(i, j, k, l, m) shared(nps) \ schedule(dynamic) @@ -2227,13 +2234,15 @@ static LA_CF_T **LA_FF(probabilistic_sparse_dense_echelon_form)( LA_CF_T *tmp; hm_t npc; len_t bctr = 0; + msolve_rng_t rng; + msolve_rng_seed(&rng, msolve_rng_derive_seed(rng_base, i)); while (bctr < nrbl) { npc = 0; /* fill random value array */ for (j = 0; j < nrbl; ++j) { do { - mull[j] = (int64_t)rand() % fc; + mull[j] = (int64_t)msolve_rng_rand(&rng) % fc; } while (mull[j] == 0); } /* generate one dense row as random linear combination diff --git a/src/neogb/libneogb.h b/src/neogb/libneogb.h index 031a88e9..09f29f91 100644 --- a/src/neogb/libneogb.h +++ b/src/neogb/libneogb.h @@ -36,5 +36,6 @@ #include "sort_r.h" #include "tools.h" #include "update.h" +#include "rng.h" #endif diff --git a/src/neogb/mt.c b/src/neogb/mt.c new file mode 100644 index 00000000..ea9a8adc --- /dev/null +++ b/src/neogb/mt.c @@ -0,0 +1,153 @@ +/* This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + This program is distributed in the hope that it will be useful, but + WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + General Public License for more details. You should have received + a copy of the GNU General Public License along with this program; + if not, write to the Free Foundation, Inc., 59 Temple Place, Suite + 330, Boston, MA 02111-1307 USA + + Original implementation was copyright (C) 1997 Makoto Matsumoto and + Takuji Nishimura. Coded by Takuji Nishimura, considering the + suggestions by Topher Cooper and Marc Rieffel in July-Aug. 1997, "A + C-program for MT19937: Integer version (1998/4/6)" + + This implementation copyright (C) 1998 Brian Gough. I reorganized + the code to use the module framework of GSL. The license on this + implementation was changed from LGPL to GPL, following paragraph 3 + of the LGPL, version 2. + + Update: + + The seeding procedure has been updated to match the 10/99 release + of MT19937. + + Update: + + The seeding procedure has been updated again to match the 2002 + release of MT19937 + + The original code included the comment: "When you use this, send an + email to: matumoto@math.keio.ac.jp with an appropriate reference to + your work". + + Makoto Matsumoto has a web page with more information about the + generator, http://www.math.keio.ac.jp/~matumoto/emt.html. + + The paper below has details of the algorithm. + + From: Makoto Matsumoto and Takuji Nishimura, "Mersenne Twister: A + 623-dimensionally equidistributerd uniform pseudorandom number + generator". ACM Transactions on Modeling and Computer Simulation, + Vol. 8, No. 1 (Jan. 1998), Pages 3-30 + + You can obtain the paper directly from Makoto Matsumoto's web page. + + The period of this generator is 2^{19937} - 1. + +*/ + +/* This file is derived from rng/mt.c of GSL 1.9, the last GSL release + under GPL version 2 or later. + + It was modified for msolve on 2026-10-01 as follows: + + - The GSL includes and the gsl_rng_type framework were removed. + - The 1998 and 1999 seeding procedures were removed. + - The state type mt_state_t and the functions mt_get and mt_set are + declared in mt.h. Hence, mt_get and mt_set are no longer static. + - The state words, as well as all other unsigned long variables, + parameters and return values, are uint32_t instead of unsigned long. + The UL constants are then U constants, and the masks with 0xffffffff + in mt_set were removed. + - The local macros are undefined at the end of the file. + +*/ + +#include "mt.h" + +#define N 624 /* Period parameters */ +#define M 397 + +/* most significant w-r bits */ +static const uint32_t UPPER_MASK = 0x80000000U; + +/* least significant r bits */ +static const uint32_t LOWER_MASK = 0x7fffffffU; + +uint32_t +mt_get (void *vstate) +{ + mt_state_t *state = (mt_state_t *) vstate; + + uint32_t k ; + uint32_t *const mt = state->mt; + +#define MAGIC(y) (((y)&0x1) ? 0x9908b0dfU : 0) + + if (state->mti >= N) + { /* generate N words at one time */ + int kk; + + for (kk = 0; kk < N - M; kk++) + { + uint32_t y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK); + mt[kk] = mt[kk + M] ^ (y >> 1) ^ MAGIC(y); + } + for (; kk < N - 1; kk++) + { + uint32_t y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK); + mt[kk] = mt[kk + (M - N)] ^ (y >> 1) ^ MAGIC(y); + } + + { + uint32_t y = (mt[N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK); + mt[N - 1] = mt[M - 1] ^ (y >> 1) ^ MAGIC(y); + } + + state->mti = 0; + } + + /* Tempering */ + + k = mt[state->mti]; + k ^= (k >> 11); + k ^= (k << 7) & 0x9d2c5680U; + k ^= (k << 15) & 0xefc60000U; + k ^= (k >> 18); + + state->mti++; + + return k; +} + +void +mt_set (void *vstate, uint32_t s) +{ + mt_state_t *state = (mt_state_t *) vstate; + int i; + + if (s == 0) + s = 4357; /* the default seed is 4357 */ + + state->mt[0]= s; + + for (i = 1; i < N; i++) + { + /* See Knuth's "Art of Computer Programming" Vol. 2, 3rd + Ed. p.106 for multiplier. */ + + state->mt[i] = + (1812433253U * (state->mt[i-1] ^ (state->mt[i-1] >> 30)) + i); + } + + state->mti = i; +} + +#undef MAGIC +#undef M +#undef N diff --git a/src/neogb/mt.h b/src/neogb/mt.h new file mode 100644 index 00000000..5e87ddab --- /dev/null +++ b/src/neogb/mt.h @@ -0,0 +1,78 @@ +/* This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + This program is distributed in the hope that it will be useful, but + WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + General Public License for more details. You should have received + a copy of the GNU General Public License along with this program; + if not, write to the Free Foundation, Inc., 59 Temple Place, Suite + 330, Boston, MA 02111-1307 USA + + Original implementation was copyright (C) 1997 Makoto Matsumoto and + Takuji Nishimura. Coded by Takuji Nishimura, considering the + suggestions by Topher Cooper and Marc Rieffel in July-Aug. 1997, "A + C-program for MT19937: Integer version (1998/4/6)" + + This implementation copyright (C) 1998 Brian Gough. I reorganized + the code to use the module framework of GSL. The license on this + implementation was changed from LGPL to GPL, following paragraph 3 + of the LGPL, version 2. + + Update: + + The seeding procedure has been updated to match the 10/99 release + of MT19937. + + Update: + + The seeding procedure has been updated again to match the 2002 + release of MT19937 + + The original code included the comment: "When you use this, send an + email to: matumoto@math.keio.ac.jp with an appropriate reference to + your work". + + Makoto Matsumoto has a web page with more information about the + generator, http://www.math.keio.ac.jp/~matumoto/emt.html. + + The paper below has details of the algorithm. + + From: Makoto Matsumoto and Takuji Nishimura, "Mersenne Twister: A + 623-dimensionally equidistributerd uniform pseudorandom number + generator". ACM Transactions on Modeling and Computer Simulation, + Vol. 8, No. 1 (Jan. 1998), Pages 3-30 + + You can obtain the paper directly from Makoto Matsumoto's web page. + + The period of this generator is 2^{19937} - 1. + +*/ + +/* Interface of the Mersenne Twister MT19937 in mt.c, derived from rng/mt.c + of GSL 1.9 (see mt.c for the modifications). + +*/ + +#ifndef MT19937_MT_H +#define MT19937_MT_H + +#include + +/* state of the generator */ +typedef struct + { + uint32_t mt[624]; + int32_t mti; + } +mt_state_t; + +/* seeds the generator state with s */ +void mt_set (void *state, uint32_t s); + +/* returns the next 32 bit output of the generator state */ +uint32_t mt_get (void *vstate); + +#endif diff --git a/src/neogb/rng.c b/src/neogb/rng.c new file mode 100644 index 00000000..f8d054dc --- /dev/null +++ b/src/neogb/rng.c @@ -0,0 +1,78 @@ +/* This file is part of msolve. + * + * msolve is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 2 of the License, or + * (at your option) any later version. + * + * msolve is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with msolve. If not, see + * + * Authors: + * Jérémy Berthomieu + * Christian Eder + * Mohab Safey El Din */ + + +#include "rng.h" + +#include "mt.h" + +#if defined(__STDC_VERSION__) && __STDC_VERSION__ >= 201112L +#define MSOLVE_THREAD_LOCAL _Thread_local +#else +#define MSOLVE_THREAD_LOCAL __thread +#endif + +static MSOLVE_THREAD_LOCAL msolve_rng_t msolve_rng_thread; +static MSOLVE_THREAD_LOCAL int msolve_rng_thread_seeded = 0; + +void msolve_rng_seed(msolve_rng_t *rng, uint32_t seed) +{ + mt_set(rng, seed); +} + +uint32_t msolve_rng_get(msolve_rng_t *rng) +{ + return mt_get(rng); +} + +int32_t msolve_rng_rand(msolve_rng_t *rng) +{ + return (int32_t)(msolve_rng_get(rng) >> 1); +} + +uint32_t msolve_rng_derive_seed(uint32_t base, uint32_t idx) +{ + return base + idx; +} + +msolve_rng_t *msolve_rng_current(void) +{ + if (!msolve_rng_thread_seeded) { + msolve_rng_seed(&msolve_rng_thread, 0); + msolve_rng_thread_seeded = 1; + } + return &msolve_rng_thread; +} + +void msolve_srand(uint32_t seed) +{ + msolve_rng_seed(&msolve_rng_thread, seed); + msolve_rng_thread_seeded = 1; +} + +int32_t msolve_rand(void) +{ + return msolve_rng_rand(msolve_rng_current()); +} + +uint32_t msolve_rand_u32(void) +{ + return msolve_rng_get(msolve_rng_current()); +} diff --git a/src/neogb/rng.h b/src/neogb/rng.h new file mode 100644 index 00000000..cbaa5373 --- /dev/null +++ b/src/neogb/rng.h @@ -0,0 +1,75 @@ +/* This file is part of msolve. + * + * msolve is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 2 of the License, or + * (at your option) any later version. + * + * msolve is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with msolve. If not, see + * + * Authors: + * Jérémy Berthomieu + * Christian Eder + * Mohab Safey El Din */ + + +#ifndef MSOLVE_RNG_H +#define MSOLVE_RNG_H + +#include +#include "mt.h" + +/* Platform independent pseudo-random number generation based on the + * Mersenne Twister MT19937 (M. Matsumoto and T. Nishimura, 1998, with the + * 2002 seeding procedure). The generator is GSL's gsl_rng_mt19937 taken from + * GSL 1.9 (see mt.c): for a given seed the generated sequence is + * identical on all platforms and coincides with the one of GSL (including + * the convention that seed 0 is replaced by 4357). + * + * Each thread owns its own generator (see msolve_rng_current()), the one + * of the main thread is seeded by msolve_srand(). Code drawing random + * numbers inside an OpenMP parallel region must not depend on the thread + * scheduling: it should use a local generator seeded via + * msolve_rng_derive_seed() from a base seed drawn before entering the + * parallel region and from the loop index. */ + +/* largest value returned by msolve_rand() and msolve_rng_rand() */ +#define MSOLVE_RAND_MAX 0x7fffffff + +/* state of the generator, see mt.h */ +typedef mt_state_t msolve_rng_t; + +/* seeds rng with seed, seed 0 is replaced by 4357 as in GSL */ +void msolve_rng_seed(msolve_rng_t *rng, uint32_t seed); + +/* returns the next 32 bit output of rng, i.e. the same value as + * gsl_rng_get() for gsl_rng_mt19937 */ +uint32_t msolve_rng_get(msolve_rng_t *rng); + +/* returns the next output of rng in [0, MSOLVE_RAND_MAX] (the upper 31 bits + * of msolve_rng_get()), as a drop-in replacement for rand() */ +int32_t msolve_rng_rand(msolve_rng_t *rng); + +/* returns the seed base + idx (mod 2^32) for a child generator, depending + * only on base and idx */ +uint32_t msolve_rng_derive_seed(uint32_t base, uint32_t idx); + +/* returns the generator of the calling thread; if it was never seeded it + * is seeded with 0 (i.e. 4357) */ +msolve_rng_t *msolve_rng_current(void); + +/* seeds the generator of the calling thread */ +void msolve_srand(uint32_t seed); + +/* draws from the generator of the calling thread: msolve_rand() returns a + * value in [0, MSOLVE_RAND_MAX], msolve_rand_u32() a full 32 bit value */ +int32_t msolve_rand(void); +uint32_t msolve_rand_u32(void); + +#endif diff --git a/test/diff/diff_hankel-degenerate-16.sh b/test/diff/diff_hankel-degenerate-16.sh index c8228462..24ad6952 100755 --- a/test/diff/diff_hankel-degenerate-16.sh +++ b/test/diff/diff_hankel-degenerate-16.sh @@ -1,12 +1,12 @@ #!/usr/bin/env bash -# With this seed and glibc's rand(), the Hankel matrix used to compute the +# With this seed, the Hankel matrix used to compute the # parametrizations has a singular sub-block of size dim - 1, so that the # computation must fall back to another method. file=cyclic5-16 -SEED=1790795968 +SEED=37892 source test/diff/diff_source.sh diff --git a/test/neogb/rng/mt19937.c b/test/neogb/rng/mt19937.c new file mode 100644 index 00000000..72789ece --- /dev/null +++ b/test/neogb/rng/mt19937.c @@ -0,0 +1,42 @@ +#include +#include "../../../src/neogb/rng.h" + +/* returns the n-th output of the generator seeded with seed */ +static uint32_t nth_output(uint32_t seed, uint32_t n) +{ + msolve_rng_t rng; + uint32_t i, v = 0; + msolve_rng_seed(&rng, seed); + for (i = 0; i < n; ++i) { + v = msolve_rng_get(&rng); + } + return v; +} + +int main(void) +{ + int ret = 0; + /* reference value of GSL's test suite for gsl_rng_mt19937 */ + if (nth_output(4357, 1000) != 1186927261U) { + fprintf(stderr, "seed 4357: wrong 1000th output\n"); + ret = 1; + } + /* seed 0 is replaced by 4357, as in GSL */ + if (nth_output(0, 1000) != 1186927261U) { + fprintf(stderr, "seed 0: wrong 1000th output\n"); + ret = 1; + } + /* reference value required by the C++ standard for std::mt19937 */ + if (nth_output(5489, 10000) != 4123659995U) { + fprintf(stderr, "seed 5489: wrong 10000th output\n"); + ret = 1; + } + /* the thread generator must follow msolve_srand */ + msolve_srand(5489); + if (msolve_rand_u32() != nth_output(5489, 1) + || msolve_rand() != (int32_t)(nth_output(5489, 2) >> 1)) { + fprintf(stderr, "msolve_srand/msolve_rand mismatch\n"); + ret = 1; + } + return ret; +}