checkin of 0.3.0
[pspp] / 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