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