From 1c66d7a6a1f1821cc17ab3489ae3e8974fcfc14a 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. 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 | 15 +++- src/neogb/la_ff_32.c | 15 +++- src/neogb/la_ff_8.c | 15 +++- 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/neogb/rng/mt19937.c | 42 +++++++++++ 17 files changed, 501 insertions(+), 33 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 e2942621..8b764751 100644 --- a/.gitignore +++ b/.gitignore @@ -51,6 +51,7 @@ fglm_build_matrixn_nonradical-shape-31 fglm_build_matrixn_radical-shape-31 line_endings_support neogb_io +neogb_rng msolve-*.tar.gz msolve-*/ src/msolve/msolve-version.h diff --git a/Makefile.am b/Makefile.am index 41023691..a44348c1 100644 --- a/Makefile.am +++ b/Makefile.am @@ -8,6 +8,7 @@ bin_PROGRAMS = msolve msolve_SOURCES = src/msolve/main.c check_PROGRAMS = neogb_io \ + neogb_rng \ fglm_build_matrixn_radical-shape-31 \ fglm_build_matrixn_nonradical-shape-31 \ fglm_build_matrixn_nonradical-radicalshape-31 \ @@ -82,6 +83,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_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 85e8c66e..109cdf52 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 bda35a8b..3ba01e98 100644 --- a/src/msolve/lifting-gb.c +++ b/src/msolve/lifting-gb.c @@ -1309,9 +1309,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; @@ -1426,7 +1426,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 134e930c..501c0e8e 100644 --- a/src/msolve/main.c +++ b/src/msolve/main.c @@ -543,14 +543,14 @@ int main(int argc, char **argv){ &precision, &refine, &isolate, &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 cb82a156..c7026e94 100644 --- a/src/msolve/msolve.c +++ b/src/msolve/msolve.c @@ -735,9 +735,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; @@ -745,9 +745,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); @@ -1855,9 +1855,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; @@ -1945,6 +1948,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; } @@ -2248,9 +2252,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; @@ -4840,9 +4844,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; @@ -5024,9 +5028,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 f9800338..f6e1c767 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 \ @@ -40,4 +42,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 fee1cf0f..e70b886b 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 9b08e8df..f4d10757 100644 --- a/src/neogb/la_ff_16.c +++ b/src/neogb/la_ff_16.c @@ -758,6 +758,7 @@ static void probabilistic_sparse_reduced_echelon_form_ff_16( /* 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) @@ -772,6 +773,8 @@ static void probabilistic_sparse_reduced_echelon_form_ff_16( /* 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)); while (bctr < nrbl) { cfp = bctr + i*rpb; sc = 0; @@ -779,7 +782,7 @@ static void probabilistic_sparse_reduced_echelon_form_ff_16( /* 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 @@ -1657,6 +1660,7 @@ static cf16_t **probabilistic_dense_linear_algebra_ff_16( (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) @@ -1671,6 +1675,8 @@ static cf16_t **probabilistic_dense_linear_algebra_ff_16( hm_t os; cf16_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; os = ncr % UNROLL; @@ -1678,7 +1684,7 @@ static cf16_t **probabilistic_dense_linear_algebra_ff_16( /* 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 @@ -1793,6 +1799,7 @@ static cf16_t **probabilistic_sparse_dense_echelon_form_ff_16( (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) @@ -1806,13 +1813,15 @@ static cf16_t **probabilistic_sparse_dense_echelon_form_ff_16( cf16_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/la_ff_32.c b/src/neogb/la_ff_32.c index 52c912b6..6008bfc7 100644 --- a/src/neogb/la_ff_32.c +++ b/src/neogb/la_ff_32.c @@ -2364,6 +2364,7 @@ static void probabilistic_sparse_reduced_echelon_form_ff_32( /* 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) @@ -2378,6 +2379,8 @@ static void probabilistic_sparse_reduced_echelon_form_ff_32( /* 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)); while (bctr < nrbl) { cfp = bctr + i*rpb; sc = 0; @@ -2385,7 +2388,7 @@ static void probabilistic_sparse_reduced_echelon_form_ff_32( /* fill random value array */ for (j = 0; j < nrbl; ++j) { do { - mull[j] = (int64_t)rand() & mask; + mull[j] = (int64_t)msolve_rng_rand(&rng) & mask; } while (mull[j] == 0); } /* generate one dense row as random linear combination @@ -3583,6 +3586,7 @@ static cf32_t **probabilistic_dense_linear_algebra_ff_32( (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) @@ -3597,6 +3601,8 @@ static cf32_t **probabilistic_dense_linear_algebra_ff_32( hm_t os; cf32_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; os = ncr % UNROLL; @@ -3604,7 +3610,7 @@ static cf32_t **probabilistic_dense_linear_algebra_ff_32( /* fill random value array */ for (j = 0; j < nrbl; ++j) { do { - mull[j] = (int64_t)rand() & mask; + mull[j] = (int64_t)msolve_rng_rand(&rng) & mask; } while (mull[j] == 0); } /* generate one dense row as random linear combination @@ -3719,6 +3725,7 @@ static cf32_t **probabilistic_sparse_dense_echelon_form_ff_32( (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) @@ -3732,13 +3739,15 @@ static cf32_t **probabilistic_sparse_dense_echelon_form_ff_32( cf32_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/la_ff_8.c b/src/neogb/la_ff_8.c index 89af2c54..17096cdd 100644 --- a/src/neogb/la_ff_8.c +++ b/src/neogb/la_ff_8.c @@ -907,6 +907,7 @@ static void probabilistic_sparse_reduced_echelon_form_ff_8( /* 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) @@ -921,6 +922,8 @@ static void probabilistic_sparse_reduced_echelon_form_ff_8( /* 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)); while (bctr < nrbl) { cfp = bctr + i*rpb; sc = 0; @@ -928,7 +931,7 @@ static void probabilistic_sparse_reduced_echelon_form_ff_8( /* 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 @@ -1786,6 +1789,7 @@ static cf8_t **probabilistic_dense_linear_algebra_ff_8( (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) @@ -1800,6 +1804,8 @@ static cf8_t **probabilistic_dense_linear_algebra_ff_8( hm_t os; cf8_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; os = ncr % UNROLL; @@ -1807,7 +1813,7 @@ static cf8_t **probabilistic_dense_linear_algebra_ff_8( /* 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 @@ -1922,6 +1928,7 @@ static cf8_t **probabilistic_sparse_dense_echelon_form_ff_8( (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) @@ -1935,13 +1942,15 @@ static cf8_t **probabilistic_sparse_dense_echelon_form_ff_8( cf8_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/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; +}