diff -r d59bea55db9b -r c445c931472f src/glprng01.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/glprng01.c Mon Dec 06 13:09:21 2010 +0100 @@ -0,0 +1,227 @@ +/* glprng01.c */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* This code is a modified version of the module GB_FLIP, a portable +* pseudo-random number generator. The original version of GB_FLIP is +* a part of The Stanford GraphBase developed by Donald E. Knuth (see +* http://www-cs-staff.stanford.edu/~knuth/sgb.html). +* +* Note that all changes concern only external names, so this modified +* version produces exactly the same results as the original version. +* +* Changes were made by Andrew Makhorin . +* +* GLPK 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 3 of the License, or +* (at your option) any later version. +* +* GLPK 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 GLPK. If not, see . +***********************************************************************/ + +#include "glpenv.h" +#include "glprng.h" + +#if 0 +int A[56] = { -1 }; +#else +#define A (rand->A) +#endif +/* pseudo-random values */ + +#if 0 +int *fptr = A; +#else +#define fptr (rand->fptr) +#endif +/* the next A value to be exported */ + +#define mod_diff(x, y) (((x) - (y)) & 0x7FFFFFFF) +/* difference modulo 2^31 */ + +static int flip_cycle(RNG *rand) +{ /* this is an auxiliary routine to do 55 more steps of the basic + recurrence, at high speed, and to reset fptr */ + int *ii, *jj; + for (ii = &A[1], jj = &A[32]; jj <= &A[55]; ii++, jj++) + *ii = mod_diff(*ii, *jj); + for (jj = &A[1]; ii <= &A[55]; ii++, jj++) + *ii = mod_diff(*ii, *jj); + fptr = &A[54]; + return A[55]; +} + +/*********************************************************************** +* NAME +* +* rng_create_rand - create pseudo-random number generator +* +* SYNOPSIS +* +* #include "glprng.h" +* RNG *rng_create_rand(void); +* +* DESCRIPTION +* +* The routine rng_create_rand creates and initializes a pseudo-random +* number generator. +* +* RETURNS +* +* The routine returns a pointer to the generator created. */ + +RNG *rng_create_rand(void) +{ RNG *rand; + int i; + rand = xmalloc(sizeof(RNG)); + A[0] = -1; + for (i = 1; i <= 55; i++) A[i] = 0; + fptr = A; + rng_init_rand(rand, 1); + return rand; +} + +/*********************************************************************** +* NAME +* +* rng_init_rand - initialize pseudo-random number generator +* +* SYNOPSIS +* +* #include "glprng.h" +* void rng_init_rand(RNG *rand, int seed); +* +* DESCRIPTION +* +* The routine rng_init_rand initializes the pseudo-random number +* generator. The parameter seed may be any integer number. Note that +* on creating the generator this routine is called with the parameter +* seed equal to 1. */ + +void rng_init_rand(RNG *rand, int seed) +{ int i; + int prev = seed, next = 1; + seed = prev = mod_diff(prev, 0); + A[55] = prev; + for (i = 21; i; i = (i + 21) % 55) + { A[i] = next; + next = mod_diff(prev, next); + if (seed & 1) + seed = 0x40000000 + (seed >> 1); + else + seed >>= 1; + next = mod_diff(next, seed); + prev = A[i]; + } + flip_cycle(rand); + flip_cycle(rand); + flip_cycle(rand); + flip_cycle(rand); + flip_cycle(rand); + return; +} + +/*********************************************************************** +* NAME +* +* rng_next_rand - obtain pseudo-random integer in the range [0, 2^31-1] +* +* SYNOPSIS +* +* #include "glprng.h" +* int rng_next_rand(RNG *rand); +* +* RETURNS +* +* The routine rng_next_rand returns a next pseudo-random integer which +* is uniformly distributed between 0 and 2^31-1, inclusive. The period +* length of the generated numbers is 2^85 - 2^30. The low order bits of +* the generated numbers are just as random as the high-order bits. */ + +int rng_next_rand(RNG *rand) +{ return + *fptr >= 0 ? *fptr-- : flip_cycle(rand); +} + +/*********************************************************************** +* NAME +* +* rng_unif_rand - obtain pseudo-random integer in the range [0, m-1] +* +* SYNOPSIS +* +* #include "glprng.h" +* int rng_unif_rand(RNG *rand, int m); +* +* RETURNS +* +* The routine rng_unif_rand returns a next pseudo-random integer which +* is uniformly distributed between 0 and m-1, inclusive, where m is any +* positive integer less than 2^31. */ + +#define two_to_the_31 ((unsigned int)0x80000000) + +int rng_unif_rand(RNG *rand, int m) +{ unsigned int t = two_to_the_31 - (two_to_the_31 % m); + int r; + xassert(m > 0); + do { r = rng_next_rand(rand); } while (t <= (unsigned int)r); + return r % m; +} + +/*********************************************************************** +* NAME +* +* rng_delete_rand - delete pseudo-random number generator +* +* SYNOPSIS +* +* #include "glprng.h" +* void rng_delete_rand(RNG *rand); +* +* DESCRIPTION +* +* The routine rng_delete_rand frees all the memory allocated to the +* specified pseudo-random number generator. */ + +void rng_delete_rand(RNG *rand) +{ xfree(rand); + return; +} + +/**********************************************************************/ + +#if 0 +/* To be sure that this modified version produces the same results as + the original version, run this validation program. */ + +int main(void) +{ RNG *rand; + int j; + rand = rng_create_rand(); + rng_init_rand(rand, -314159); + if (rng_next_rand(rand) != 119318998) + { fprintf(stderr, "Failure on the first try!\n"); + return -1; + } + for (j = 1; j <= 133; j++) rng_next_rand(rand); + if (rng_unif_rand(rand, 0x55555555) != 748103812) + { fprintf(stderr, "Failure on the second try!\n"); + return -2; + } + fprintf(stderr, "OK, the random-number generator routines seem to" + " work!\n"); + rng_delete_rand(rand); + return 0; +} +#endif + +/* eof */