1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glprng01.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,227 @@
1.4 +/* glprng01.c */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* This code is a modified version of the module GB_FLIP, a portable
1.10 +* pseudo-random number generator. The original version of GB_FLIP is
1.11 +* a part of The Stanford GraphBase developed by Donald E. Knuth (see
1.12 +* http://www-cs-staff.stanford.edu/~knuth/sgb.html).
1.13 +*
1.14 +* Note that all changes concern only external names, so this modified
1.15 +* version produces exactly the same results as the original version.
1.16 +*
1.17 +* Changes were made by Andrew Makhorin <mao@gnu.org>.
1.18 +*
1.19 +* GLPK is free software: you can redistribute it and/or modify it
1.20 +* under the terms of the GNU General Public License as published by
1.21 +* the Free Software Foundation, either version 3 of the License, or
1.22 +* (at your option) any later version.
1.23 +*
1.24 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.25 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.26 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.27 +* License for more details.
1.28 +*
1.29 +* You should have received a copy of the GNU General Public License
1.30 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.31 +***********************************************************************/
1.32 +
1.33 +#include "glpenv.h"
1.34 +#include "glprng.h"
1.35 +
1.36 +#if 0
1.37 +int A[56] = { -1 };
1.38 +#else
1.39 +#define A (rand->A)
1.40 +#endif
1.41 +/* pseudo-random values */
1.42 +
1.43 +#if 0
1.44 +int *fptr = A;
1.45 +#else
1.46 +#define fptr (rand->fptr)
1.47 +#endif
1.48 +/* the next A value to be exported */
1.49 +
1.50 +#define mod_diff(x, y) (((x) - (y)) & 0x7FFFFFFF)
1.51 +/* difference modulo 2^31 */
1.52 +
1.53 +static int flip_cycle(RNG *rand)
1.54 +{ /* this is an auxiliary routine to do 55 more steps of the basic
1.55 + recurrence, at high speed, and to reset fptr */
1.56 + int *ii, *jj;
1.57 + for (ii = &A[1], jj = &A[32]; jj <= &A[55]; ii++, jj++)
1.58 + *ii = mod_diff(*ii, *jj);
1.59 + for (jj = &A[1]; ii <= &A[55]; ii++, jj++)
1.60 + *ii = mod_diff(*ii, *jj);
1.61 + fptr = &A[54];
1.62 + return A[55];
1.63 +}
1.64 +
1.65 +/***********************************************************************
1.66 +* NAME
1.67 +*
1.68 +* rng_create_rand - create pseudo-random number generator
1.69 +*
1.70 +* SYNOPSIS
1.71 +*
1.72 +* #include "glprng.h"
1.73 +* RNG *rng_create_rand(void);
1.74 +*
1.75 +* DESCRIPTION
1.76 +*
1.77 +* The routine rng_create_rand creates and initializes a pseudo-random
1.78 +* number generator.
1.79 +*
1.80 +* RETURNS
1.81 +*
1.82 +* The routine returns a pointer to the generator created. */
1.83 +
1.84 +RNG *rng_create_rand(void)
1.85 +{ RNG *rand;
1.86 + int i;
1.87 + rand = xmalloc(sizeof(RNG));
1.88 + A[0] = -1;
1.89 + for (i = 1; i <= 55; i++) A[i] = 0;
1.90 + fptr = A;
1.91 + rng_init_rand(rand, 1);
1.92 + return rand;
1.93 +}
1.94 +
1.95 +/***********************************************************************
1.96 +* NAME
1.97 +*
1.98 +* rng_init_rand - initialize pseudo-random number generator
1.99 +*
1.100 +* SYNOPSIS
1.101 +*
1.102 +* #include "glprng.h"
1.103 +* void rng_init_rand(RNG *rand, int seed);
1.104 +*
1.105 +* DESCRIPTION
1.106 +*
1.107 +* The routine rng_init_rand initializes the pseudo-random number
1.108 +* generator. The parameter seed may be any integer number. Note that
1.109 +* on creating the generator this routine is called with the parameter
1.110 +* seed equal to 1. */
1.111 +
1.112 +void rng_init_rand(RNG *rand, int seed)
1.113 +{ int i;
1.114 + int prev = seed, next = 1;
1.115 + seed = prev = mod_diff(prev, 0);
1.116 + A[55] = prev;
1.117 + for (i = 21; i; i = (i + 21) % 55)
1.118 + { A[i] = next;
1.119 + next = mod_diff(prev, next);
1.120 + if (seed & 1)
1.121 + seed = 0x40000000 + (seed >> 1);
1.122 + else
1.123 + seed >>= 1;
1.124 + next = mod_diff(next, seed);
1.125 + prev = A[i];
1.126 + }
1.127 + flip_cycle(rand);
1.128 + flip_cycle(rand);
1.129 + flip_cycle(rand);
1.130 + flip_cycle(rand);
1.131 + flip_cycle(rand);
1.132 + return;
1.133 +}
1.134 +
1.135 +/***********************************************************************
1.136 +* NAME
1.137 +*
1.138 +* rng_next_rand - obtain pseudo-random integer in the range [0, 2^31-1]
1.139 +*
1.140 +* SYNOPSIS
1.141 +*
1.142 +* #include "glprng.h"
1.143 +* int rng_next_rand(RNG *rand);
1.144 +*
1.145 +* RETURNS
1.146 +*
1.147 +* The routine rng_next_rand returns a next pseudo-random integer which
1.148 +* is uniformly distributed between 0 and 2^31-1, inclusive. The period
1.149 +* length of the generated numbers is 2^85 - 2^30. The low order bits of
1.150 +* the generated numbers are just as random as the high-order bits. */
1.151 +
1.152 +int rng_next_rand(RNG *rand)
1.153 +{ return
1.154 + *fptr >= 0 ? *fptr-- : flip_cycle(rand);
1.155 +}
1.156 +
1.157 +/***********************************************************************
1.158 +* NAME
1.159 +*
1.160 +* rng_unif_rand - obtain pseudo-random integer in the range [0, m-1]
1.161 +*
1.162 +* SYNOPSIS
1.163 +*
1.164 +* #include "glprng.h"
1.165 +* int rng_unif_rand(RNG *rand, int m);
1.166 +*
1.167 +* RETURNS
1.168 +*
1.169 +* The routine rng_unif_rand returns a next pseudo-random integer which
1.170 +* is uniformly distributed between 0 and m-1, inclusive, where m is any
1.171 +* positive integer less than 2^31. */
1.172 +
1.173 +#define two_to_the_31 ((unsigned int)0x80000000)
1.174 +
1.175 +int rng_unif_rand(RNG *rand, int m)
1.176 +{ unsigned int t = two_to_the_31 - (two_to_the_31 % m);
1.177 + int r;
1.178 + xassert(m > 0);
1.179 + do { r = rng_next_rand(rand); } while (t <= (unsigned int)r);
1.180 + return r % m;
1.181 +}
1.182 +
1.183 +/***********************************************************************
1.184 +* NAME
1.185 +*
1.186 +* rng_delete_rand - delete pseudo-random number generator
1.187 +*
1.188 +* SYNOPSIS
1.189 +*
1.190 +* #include "glprng.h"
1.191 +* void rng_delete_rand(RNG *rand);
1.192 +*
1.193 +* DESCRIPTION
1.194 +*
1.195 +* The routine rng_delete_rand frees all the memory allocated to the
1.196 +* specified pseudo-random number generator. */
1.197 +
1.198 +void rng_delete_rand(RNG *rand)
1.199 +{ xfree(rand);
1.200 + return;
1.201 +}
1.202 +
1.203 +/**********************************************************************/
1.204 +
1.205 +#if 0
1.206 +/* To be sure that this modified version produces the same results as
1.207 + the original version, run this validation program. */
1.208 +
1.209 +int main(void)
1.210 +{ RNG *rand;
1.211 + int j;
1.212 + rand = rng_create_rand();
1.213 + rng_init_rand(rand, -314159);
1.214 + if (rng_next_rand(rand) != 119318998)
1.215 + { fprintf(stderr, "Failure on the first try!\n");
1.216 + return -1;
1.217 + }
1.218 + for (j = 1; j <= 133; j++) rng_next_rand(rand);
1.219 + if (rng_unif_rand(rand, 0x55555555) != 748103812)
1.220 + { fprintf(stderr, "Failure on the second try!\n");
1.221 + return -2;
1.222 + }
1.223 + fprintf(stderr, "OK, the random-number generator routines seem to"
1.224 + " work!\n");
1.225 + rng_delete_rand(rand);
1.226 + return 0;
1.227 +}
1.228 +#endif
1.229 +
1.230 +/* eof */