(rng_get_double) Fix always-returning-zero bug in my preferred way,
[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 "random.h"
22 #include <assert.h>
23 #include <inttypes.h>
24 #include <limits.h>
25 #include <math.h>
26 #include <stdlib.h>
27 #include <time.h>
28 #include "alloc.h"
29 #include "magic.h"
30 #include "settings.h"
31
32 /* Random number generator. */
33 struct rng 
34   {
35     /* RC4-based random bytes. */
36     uint8_t s[256];
37     uint8_t i, j;
38
39     /* Normal distribution. */
40     double next_normal;
41   };
42
43 /* Creates a new random number generator, seeds it based on
44    the current time, and returns it. */
45 struct rng *
46 rng_create (void) 
47 {
48   struct rng *rng;
49   static time_t t=0;
50
51   rng = xmalloc (sizeof *rng);
52   if (t == 0 || set_seed == NOT_LONG)
53   {
54     if (set_seed == NOT_LONG) 
55       time (&t);
56     else
57       t = set_seed;
58   }
59   else
60     t++;
61   rng_seed (rng, &t, sizeof t);
62   rng->next_normal = NOT_DOUBLE;
63   return rng;
64 }
65
66 /* Destroys RNG. */
67 void
68 rng_destroy (struct rng *rng) 
69 {
70   free (rng);
71 }
72
73 /* Swap bytes. */
74 static void
75 swap_byte (uint8_t *a, uint8_t *b) 
76 {
77   uint8_t t = *a;
78   *a = *b;
79   *b = t;
80 }
81
82 /* Seeds RNG based on the SIZE bytes in BUF.
83    At most the first 256 bytes of BUF are used. */
84 void
85 rng_seed (struct rng *rng, const void *key_, size_t size) 
86 {
87   const uint8_t *key = key_;
88   size_t key_idx;
89   uint8_t *s;
90   int i, j;
91
92   assert (rng != NULL);
93
94   s = rng->s;
95   rng->i = rng->j = 0;
96   for (i = 0; i < 256; i++) 
97     s[i] = i;
98   for (key_idx = 0, i = j = 0; i < 256; i++) 
99     {
100       j = (j + s[i] + key[key_idx]) & 255;
101       swap_byte (s + i, s + j);
102       if (++key_idx >= size)
103         key_idx = 0;
104     }
105 }
106
107 /* Reads SIZE random bytes from RNG into BUF. */
108 void
109 rng_get_bytes (struct rng *rng, void *buf_, size_t size) 
110 {
111   uint8_t *buf = buf_;
112   uint8_t *s;
113   uint8_t i, j;
114
115   assert (rng != 0);
116
117   s = rng->s;
118   i = rng->i;
119   j = rng->j;
120   while (size-- > 0) 
121     {
122       i += 1;
123       j += s[i];
124       swap_byte (s + i, s + j);
125       *buf++ = s[(s[i] + s[j]) & 255];
126     }
127   rng->i = i;
128   rng->j = j;
129 }
130
131 /* Returns a random int in the range [0, INT_MAX]. */
132 int
133 rng_get_int (struct rng *rng) 
134 {
135   int value;
136
137   do 
138     {
139       rng_get_bytes (rng, &value, sizeof value);
140       value = abs (value);
141     }
142   while (value < 0);
143    
144   return value;
145 }
146
147 /* Returns a random unsigned in the range [0, UINT_MAX]. */
148 unsigned
149 rng_get_unsigned (struct rng *rng) 
150 {
151   unsigned value;
152
153   rng_get_bytes (rng, &value, sizeof value);
154   return value;
155 }
156
157 /* Returns a random number from the uniform distribution with
158    range [0,1). */
159 double
160 rng_get_double (struct rng *rng) 
161 {
162   for (;;) 
163     {
164       unsigned long ulng;
165       double dbl;
166   
167       rng_get_bytes (rng, &ulng, sizeof ulng);
168       dbl = ulng / (ULONG_MAX + 1.0);
169       if (dbl >= 0 && dbl < 1)
170         return dbl;
171     }
172 }
173
174 /* Returns a random number from the distribution with mean 0 and
175    standard deviation 1.  (Multiply the result by the desired
176    standard deviation, then add the desired mean.) */
177 double 
178 rng_get_double_normal (struct rng *rng)
179 {
180   /* Knuth, _The Art of Computer Programming_, Vol. 2, 3.4.1C,
181      Algorithm P. */
182   double this_normal;
183   
184   if (rng->next_normal != NOT_DOUBLE)
185     {
186       this_normal = rng->next_normal;
187       rng->next_normal = NOT_DOUBLE;
188     }
189   else 
190     {
191       double v1, v2, s;
192       
193       do
194         {
195           double u1 = rng_get_double (rng);
196           double u2 = rng_get_double (rng);
197           v1 = 2.0 * u1 - 1.0;
198           v2 = 2.0 * u2 - 1.0;
199           s = v1 * v1 + v2 * v2;
200         }
201       while (s >= 1);
202
203       this_normal = v1 * sqrt (-2. * log (s) / s);
204       rng->next_normal = v2 * sqrt (-2. * log (s) / s); 
205     }
206   
207   return this_normal;
208 }
209
210 /* Gets an initialized RNG for use in PSPP transformations and
211    procedures. */
212 struct rng *
213 pspp_rng (void)
214 {
215   static struct rng *rng;
216
217   if (rng == NULL)
218     rng = rng_create ();
219   return rng;
220 }