praat/external/glpk/glprng01.c

/* 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 <mao@gnu.org>.
*
*  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 <http://www.gnu.org/licenses/>.
***********************************************************************/

#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 */