src/random.c

   1 /* PSPP - computes sample statistics.
   2    Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
   3    Written by Ben Pfaff <blp@gnu.org>.
   4
   5    This program is free software; you can redistribute it and/or
   6    modify it under the terms of the GNU General Public License as
   7    published by the Free Software Foundation; either version 2 of the
   8    License, or (at your option) any later version.
   9
  10    This program is distributed in the hope that it will be useful, but
  11    WITHOUT ANY WARRANTY; without even the implied warranty of
  12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13    General Public License for more details.
  14
  15    You should have received a copy of the GNU General Public License
  16    along with this program; if not, write to the Free Software
  17    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  18    02111-1307, USA. */
  19
  20 #include <config.h>
  21 #include <math.h>
  22 #include <stdlib.h>
  23 #include <time.h>
  24 #include "magic.h"
  25 #include "random.h"
  26 #include "settings.h"
  27
  28 /* Deal with broken system random number generator. */
  29 #if HAVE_GOOD_RANDOM
  30 #define real_rand rand
  31 #define real_srand srand
  32 #define REAL_RAND_MAX RAND_MAX
  33 #else /* !HAVE_GOOD_RANDOM */
  34 #define REAL_RAND_MAX 32767
  35
  36 /* Some systems are so broken that they do not supply a value for
  37    RAND_MAX.  There is absolutely no reliable way to determine this
  38    value, either.  So we must supply our own.  This one is the one
  39    presented in the ANSI C standard as strictly compliant. */
  40 static unsigned long int next = 1;
  41
  42 int
  43 real_rand (void)
  44 {
  45   next = next * 1103515245 + 12345;
  46   return (unsigned int)(next / 65536) % 32768;
  47 }
  48
  49 void
  50 real_srand (unsigned int seed)
  51 {
  52   next = seed;
  53 }
  54 #endif /* !HAVE_GOOD_RANDOM */
  55
  56 /* The random number generator here is an implementation in C of
  57    Knuth's Algorithm 3.2.2B (Randomizing by Shuffling) in _The Art of
  58    Computer Programming_, Vol. 2. */
  59
  60 #define k 13
  61 static int V[k];
  62 static int Y;
  63
  64 static double X2;
  65
  66 /* Initializes the random number generator.  Should be called once by
  67    every cmd_*() that uses random numbers.  Note that this includes
  68    all procedures that use expressions since they may generate random
  69    numbers. */
  70 void
  71 setup_randomize (void)
  72 {
  73   static time_t curtime;
  74   int i;
  75
  76   if (set_seed == NOT_LONG)
  77     {
  78       if (!curtime)
  79         time (&curtime);
  80       real_srand (curtime++);
  81     }
  82   else
  83     real_srand (set_seed);
  84
  85   set_seed_used = 1;
  86
  87   for (i = 0; i < k; i++)
  88     V[i] = real_rand ();
  89   Y = real_rand ();
  90   X2 = NOT_DOUBLE;
  91 }
  92
  93 /* Standard shuffling procedure for increasing randomness of the ANSI
  94    C random number generator. Returns a random number R where 0 <= R
  95    <= RAND_MAX. */
  96 inline int
  97 shuffle (void)
  98 {
  99   int j = k * Y / RAND_MAX;
 100   Y = V[j];
 101   V[j] = real_rand ();
 102   return Y;
 103 }
 104
 105 /* Returns a random number R where 0 <= R <= X. */
 106 double
 107 rand_uniform (double x)
 108 {
 109   return ((double) shuffle ()) / (((double) RAND_MAX) / x);
 110 }
 111
 112 /* Returns a random number from the distribution with mean 0 and
 113    standard deviation X.  This uses algorithm P in section 3.4.1C of
 114    Knuth's _Art of Computer Programming_, Vol 2. */
 115 double
 116 rand_normal (double x)
 117 {
 118   double U1, U2;
 119   double V1, V2;
 120   double S;
 121   double X1;
 122
 123   if (X2 != NOT_DOUBLE)
 124     {
 125       double t = X2;
 126       X2 = NOT_DOUBLE;
 127       return t * x;
 128     }
 129   do
 130     {
 131       U1 = ((double) shuffle ()) / RAND_MAX;
 132       U2 = ((double) shuffle ()) / RAND_MAX;
 133       V1 = 2 * U1 - 1;
 134       V2 = 2 * U2 - 1;
 135       S = V1 * V1 + V2 * V2;
 136     }
 137   while (S >= 1);
 138   X1 = V1 * sqrt (-2. * log (S) / S);
 139   X2 = V2 * sqrt (-2. * log (S) / S);
 140   return X1 * x;
 141 }
 142
 143 /* Returns a random integer R, where 0 <= R < X. */
 144 int
 145 rand_simple (int x)
 146 {
 147   return shuffle () % x;
 148 }
 149