alpar@9: /* glprng01.c */ alpar@9: alpar@9: /*********************************************************************** alpar@9: * This code is part of GLPK (GNU Linear Programming Kit). alpar@9: * alpar@9: * This code is a modified version of the module GB_FLIP, a portable alpar@9: * pseudo-random number generator. The original version of GB_FLIP is alpar@9: * a part of The Stanford GraphBase developed by Donald E. Knuth (see alpar@9: * http://www-cs-staff.stanford.edu/~knuth/sgb.html). alpar@9: * alpar@9: * Note that all changes concern only external names, so this modified alpar@9: * version produces exactly the same results as the original version. alpar@9: * alpar@9: * Changes were made by Andrew Makhorin . alpar@9: * alpar@9: * GLPK is free software: you can redistribute it and/or modify it alpar@9: * under the terms of the GNU General Public License as published by alpar@9: * the Free Software Foundation, either version 3 of the License, or alpar@9: * (at your option) any later version. alpar@9: * alpar@9: * GLPK is distributed in the hope that it will be useful, but WITHOUT alpar@9: * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY alpar@9: * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public alpar@9: * License for more details. alpar@9: * alpar@9: * You should have received a copy of the GNU General Public License alpar@9: * along with GLPK. If not, see . alpar@9: ***********************************************************************/ alpar@9: alpar@9: #include "glpenv.h" alpar@9: #include "glprng.h" alpar@9: alpar@9: #if 0 alpar@9: int A[56] = { -1 }; alpar@9: #else alpar@9: #define A (rand->A) alpar@9: #endif alpar@9: /* pseudo-random values */ alpar@9: alpar@9: #if 0 alpar@9: int *fptr = A; alpar@9: #else alpar@9: #define fptr (rand->fptr) alpar@9: #endif alpar@9: /* the next A value to be exported */ alpar@9: alpar@9: #define mod_diff(x, y) (((x) - (y)) & 0x7FFFFFFF) alpar@9: /* difference modulo 2^31 */ alpar@9: alpar@9: static int flip_cycle(RNG *rand) alpar@9: { /* this is an auxiliary routine to do 55 more steps of the basic alpar@9: recurrence, at high speed, and to reset fptr */ alpar@9: int *ii, *jj; alpar@9: for (ii = &A[1], jj = &A[32]; jj <= &A[55]; ii++, jj++) alpar@9: *ii = mod_diff(*ii, *jj); alpar@9: for (jj = &A[1]; ii <= &A[55]; ii++, jj++) alpar@9: *ii = mod_diff(*ii, *jj); alpar@9: fptr = &A[54]; alpar@9: return A[55]; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * rng_create_rand - create pseudo-random number generator alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glprng.h" alpar@9: * RNG *rng_create_rand(void); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine rng_create_rand creates and initializes a pseudo-random alpar@9: * number generator. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * The routine returns a pointer to the generator created. */ alpar@9: alpar@9: RNG *rng_create_rand(void) alpar@9: { RNG *rand; alpar@9: int i; alpar@9: rand = xmalloc(sizeof(RNG)); alpar@9: A[0] = -1; alpar@9: for (i = 1; i <= 55; i++) A[i] = 0; alpar@9: fptr = A; alpar@9: rng_init_rand(rand, 1); alpar@9: return rand; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * rng_init_rand - initialize pseudo-random number generator alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glprng.h" alpar@9: * void rng_init_rand(RNG *rand, int seed); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine rng_init_rand initializes the pseudo-random number alpar@9: * generator. The parameter seed may be any integer number. Note that alpar@9: * on creating the generator this routine is called with the parameter alpar@9: * seed equal to 1. */ alpar@9: alpar@9: void rng_init_rand(RNG *rand, int seed) alpar@9: { int i; alpar@9: int prev = seed, next = 1; alpar@9: seed = prev = mod_diff(prev, 0); alpar@9: A[55] = prev; alpar@9: for (i = 21; i; i = (i + 21) % 55) alpar@9: { A[i] = next; alpar@9: next = mod_diff(prev, next); alpar@9: if (seed & 1) alpar@9: seed = 0x40000000 + (seed >> 1); alpar@9: else alpar@9: seed >>= 1; alpar@9: next = mod_diff(next, seed); alpar@9: prev = A[i]; alpar@9: } alpar@9: flip_cycle(rand); alpar@9: flip_cycle(rand); alpar@9: flip_cycle(rand); alpar@9: flip_cycle(rand); alpar@9: flip_cycle(rand); alpar@9: return; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * rng_next_rand - obtain pseudo-random integer in the range [0, 2^31-1] alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glprng.h" alpar@9: * int rng_next_rand(RNG *rand); alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * The routine rng_next_rand returns a next pseudo-random integer which alpar@9: * is uniformly distributed between 0 and 2^31-1, inclusive. The period alpar@9: * length of the generated numbers is 2^85 - 2^30. The low order bits of alpar@9: * the generated numbers are just as random as the high-order bits. */ alpar@9: alpar@9: int rng_next_rand(RNG *rand) alpar@9: { return alpar@9: *fptr >= 0 ? *fptr-- : flip_cycle(rand); alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * rng_unif_rand - obtain pseudo-random integer in the range [0, m-1] alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glprng.h" alpar@9: * int rng_unif_rand(RNG *rand, int m); alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * The routine rng_unif_rand returns a next pseudo-random integer which alpar@9: * is uniformly distributed between 0 and m-1, inclusive, where m is any alpar@9: * positive integer less than 2^31. */ alpar@9: alpar@9: #define two_to_the_31 ((unsigned int)0x80000000) alpar@9: alpar@9: int rng_unif_rand(RNG *rand, int m) alpar@9: { unsigned int t = two_to_the_31 - (two_to_the_31 % m); alpar@9: int r; alpar@9: xassert(m > 0); alpar@9: do { r = rng_next_rand(rand); } while (t <= (unsigned int)r); alpar@9: return r % m; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * rng_delete_rand - delete pseudo-random number generator alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glprng.h" alpar@9: * void rng_delete_rand(RNG *rand); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine rng_delete_rand frees all the memory allocated to the alpar@9: * specified pseudo-random number generator. */ alpar@9: alpar@9: void rng_delete_rand(RNG *rand) alpar@9: { xfree(rand); alpar@9: return; alpar@9: } alpar@9: alpar@9: /**********************************************************************/ alpar@9: alpar@9: #if 0 alpar@9: /* To be sure that this modified version produces the same results as alpar@9: the original version, run this validation program. */ alpar@9: alpar@9: int main(void) alpar@9: { RNG *rand; alpar@9: int j; alpar@9: rand = rng_create_rand(); alpar@9: rng_init_rand(rand, -314159); alpar@9: if (rng_next_rand(rand) != 119318998) alpar@9: { fprintf(stderr, "Failure on the first try!\n"); alpar@9: return -1; alpar@9: } alpar@9: for (j = 1; j <= 133; j++) rng_next_rand(rand); alpar@9: if (rng_unif_rand(rand, 0x55555555) != 748103812) alpar@9: { fprintf(stderr, "Failure on the second try!\n"); alpar@9: return -2; alpar@9: } alpar@9: fprintf(stderr, "OK, the random-number generator routines seem to" alpar@9: " work!\n"); alpar@9: rng_delete_rand(rand); alpar@9: return 0; alpar@9: } alpar@9: #endif alpar@9: alpar@9: /* eof */