1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glprng01.c	Sun Nov 06 20:59:10 2011 +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 */