Sat Dec 27 16:16:49 2003 Ben Pfaff <blp@gnu.org>
[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;
50
51   rng = xmalloc (sizeof *rng);
52   if (t == 0)
53     time (&t);
54   else
55     t++;
56   rng_seed (rng, &t, sizeof t);
57   rng->next_normal = NOT_DOUBLE;
58   return rng;
59 }
60
61 /* Destroys RNG. */
62 void
63 rng_destroy (struct rng *rng) 
64 {
65   free (rng);
66 }
67
68 /* Swap bytes. */
69 static void
70 swap_byte (uint8_t *a, uint8_t *b) 
71 {
72   uint8_t t = *a;
73   *a = *b;
74   *b = t;
75 }
76
77 /* Seeds RNG based on the SIZE bytes in BUF.
78    At most the first 256 bytes of BUF are used. */
79 void
80 rng_seed (struct rng *rng, const void *key_, size_t size) 
81 {
82   const uint8_t *key = key_;
83   size_t key_idx;
84   uint8_t *s;
85   int i, j;
86
87   assert (rng != NULL);
88
89   s = rng->s;
90   rng->i = rng->j = 0;
91   for (i = 0; i < 256; i++) 
92     s[i] = i;
93   for (key_idx = 0, i = j = 0; i < 256; i++) 
94     {
95       j = (j + s[i] + key[key_idx]) & 255;
96       swap_byte (s + i, s + j);
97       if (++key_idx >= size)
98         key_idx = 0;
99     }
100 }
101
102 /* Reads SIZE random bytes from RNG into BUF. */
103 void
104 rng_get_bytes (struct rng *rng, void *buf_, size_t size) 
105 {
106   uint8_t *buf = buf_;
107   uint8_t *s;
108   uint8_t i, j;
109
110   assert (rng != 0);
111
112   s = rng->s;
113   i = rng->i;
114   j = rng->j;
115   while (size-- > 0) 
116     {
117       i += 1;
118       j += s[i];
119       swap_byte (s + i, s + j);
120       *buf++ = s[(s[i] + s[j]) & 255];
121     }
122   rng->i = i;
123   rng->j = j;
124 }
125
126 /* Returns a random int in the range [0, INT_MAX]. */
127 int
128 rng_get_int (struct rng *rng) 
129 {
130   int value;
131
132   do 
133     {
134       rng_get_bytes (rng, &value, sizeof value);
135       value = abs (value);
136     }
137   while (value < 0);
138    
139   return value;
140 }
141
142 /* Returns a random unsigned in the range [0, UINT_MAX]. */
143 unsigned
144 rng_get_unsigned (struct rng *rng) 
145 {
146   unsigned value;
147
148   rng_get_bytes (rng, &value, sizeof value);
149   return value;
150 }
151
152 /* Returns a random number from the uniform distribution with
153    range [0,1). */
154 double
155 rng_get_double (struct rng *rng) 
156 {
157   unsigned long value;
158
159   rng_get_bytes (rng, &value, sizeof value);
160   return value / ULONG_MAX;
161 }
162
163 /* Returns a random number from the distribution with mean 0 and
164    standard deviation 1.  (Multiply the result by the desired
165    standard deviation, then add the desired mean.) */
166 double 
167 rng_get_double_normal (struct rng *rng)
168 {
169   /* Knuth, _The Art of Computer Programming_, Vol. 2, 3.4.1C,
170      Algorithm P. */
171   double this_normal;
172   
173   if (rng->next_normal != NOT_DOUBLE)
174     {
175       this_normal = rng->next_normal;
176       rng->next_normal = NOT_DOUBLE;
177     }
178   else 
179     {
180       double v1, v2, s;
181       
182       do
183         {
184           double u1 = rng_get_double (rng);
185           double u2 = rng_get_double (rng);
186           v1 = 2.0 * u1 - 1.0;
187           v2 = 2.0 * u2 - 1.0;
188           s = v1 * v1 + v2 * v2;
189         }
190       while (s >= 1);
191
192       this_normal = v1 * sqrt (-2. * log (s) / s);
193       rng->next_normal = v2 * sqrt (-2. * log (s) / s); 
194     }
195   
196   return this_normal;
197 }
198
199 /* Gets an initialized RNG for use in PSPP transformations and
200    procedures. */
201 struct rng *
202 pspp_rng (void)
203 {
204   static struct rng *rng;
205
206   if (rng == NULL)
207     rng = rng_create ();
208   return rng;
209 }