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; +}