1 /*
2 * prng.c - Portable, ISO C90 and C99 compliant high-quality
3 * pseudo-random number generator based on the alleged RC4
4 * cipher. This PRNG should be suitable for most general-purpose
5 * uses. Not recommended for cryptographic or financial
6 * purposes. Not thread-safe.
7 */
8
9 /*
10 * Copyright (c) 2004 Ben Pfaff <blp@cs.stanford.edu>.
11 * All rights reserved.
12 *
13 * Redistribution and use in source and binary forms, with or
14 * without modification, are permitted provided that the
15 * following conditions are met:
16 *
17 * 1. Redistributions of source code must retain the above
18 * copyright notice, this list of conditions and the following
19 * disclaimer.
20 *
21 * 2. Redistributions in binary form must reproduce the above
22 * copyright notice, this list of conditions and the following
23 * disclaimer in the documentation and/or other materials
24 * provided with the distribution.
25 *
26 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS
27 * IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
28 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
29 * FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
30 * SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
31 * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
32 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
33 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
34 * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
36 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
37 * THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
38 * OF SUCH DAMAGE.
39 *
40 */
41
42 #include "prng.h"
43 #include <assert.h>
44 #include <float.h>
45 #include <limits.h>
46 #include <math.h>
47 #include <time.h>
48 #include <stdio.h>
49
50 /* RC4-based pseudo-random state. */
51 static unsigned char s[256];
52 static int s_i, s_j;
53
54 /* Nonzero if PRNG has been seeded. */
55 static int seeded;
56
57 /* Swap bytes that A and B point to. */
58 #define SWAP_BYTE(A, B) \
59 do { \
60 unsigned char swap_temp = *(A); \
61 *(A) = *(B); \
62 *(B) = swap_temp; \
63 } while (0)
64
65 /* Seeds the pseudo-random number generator based on the current
66 time.
67
68 If the user calls neither this function nor prng_seed_bytes()
69 before any prng_get*() function, this function is called
70 automatically to obtain a time-based seed. */
71 long
prng_seed_time(void)72 prng_seed_time (void)
73 {
74 static time_t t;
75 if (t == 0)
76 t = time (NULL);
77 else
78 t++;
79
80 prng_seed_bytes (&t, sizeof t);
81 return((long)t);
82 }
83
84 /* Retrieves one octet from the array BYTES, which is N_BYTES in
85 size, starting at an offset of OCTET_IDX octets. BYTES is
86 treated as a circular array, so that accesses past the first
87 N_BYTES bytes wrap around to the beginning. */
88 static unsigned char
get_octet(const void * bytes_,size_t n_bytes,size_t octet_idx)89 get_octet (const void *bytes_, size_t n_bytes, size_t octet_idx)
90 {
91 const unsigned char *bytes = bytes_;
92 if (CHAR_BIT == 8)
93 return bytes[octet_idx % n_bytes];
94 else
95 {
96 size_t first_byte = octet_idx * 8 / CHAR_BIT % n_bytes;
97 size_t start_bit = octet_idx * 8 % CHAR_BIT;
98 unsigned char c = (bytes[first_byte] >> start_bit) & 255;
99
100 size_t bits_filled = CHAR_BIT - start_bit;
101 if (CHAR_BIT % 8 != 0 && bits_filled < 8)
102 {
103 size_t bits_left = 8 - bits_filled;
104 unsigned char bits_left_mask = (1u << bits_left) - 1;
105 size_t second_byte = first_byte + 1 < n_bytes ? first_byte + 1 : 0;
106
107 c |= (bytes[second_byte] & bits_left_mask) << bits_filled;
108 }
109
110 return c;
111 }
112 }
113
114 /* Seeds the pseudo-random number based on the SIZE bytes in
115 KEY. At most the first 2048 bits in KEY are used. */
116 void
prng_seed_bytes(const void * key,size_t size)117 prng_seed_bytes (const void *key, size_t size)
118 {
119 int i, j;
120
121 assert (key != NULL && size > 0);
122
123 for (i = 0; i < 256; i++)
124 s[i] = i;
125 for (i = j = 0; i < 256; i++)
126 {
127 j = (j + s[i] + get_octet (key, size, i)) & 255;
128 SWAP_BYTE (s + i, s + j);
129 }
130
131 s_i = s_j = 0;
132 seeded = 1;
133 }
134
135 /* Returns a pseudo-random integer in the range [0, 255]. */
136 unsigned char
prng_get_octet(void)137 prng_get_octet (void)
138 {
139 if (!seeded)
140 prng_seed_time ();
141
142 s_i = (s_i + 1) & 255;
143 s_j = (s_j + s[s_i]) & 255;
144 SWAP_BYTE (s + s_i, s + s_j);
145
146 return s[(s[s_i] + s[s_j]) & 255];
147 }
148
149 /* Returns a pseudo-random integer in the range [0, UCHAR_MAX]. */
150 unsigned char
prng_get_byte(void)151 prng_get_byte (void)
152 {
153 unsigned byte;
154 int bits;
155
156 byte = prng_get_octet ();
157 for (bits = 8; bits < CHAR_BIT; bits += 8)
158 byte = (byte << 8) | prng_get_octet ();
159 return byte;
160 }
161
162 /* Fills BUF with SIZE pseudo-random bytes. */
163 void
prng_get_bytes(void * buf_,size_t size)164 prng_get_bytes (void *buf_, size_t size)
165 {
166 unsigned char *buf;
167
168 for (buf = buf_; size-- > 0; buf++)
169 *buf = prng_get_byte ();
170 }
171
172 /* Returns a pseudo-random unsigned long in the range [0,
173 ULONG_MAX]. */
174 unsigned long
prng_get_ulong(void)175 prng_get_ulong (void)
176 {
177 unsigned long ulng;
178 size_t bits;
179
180 ulng = prng_get_octet ();
181 for (bits = 8; bits < CHAR_BIT * sizeof ulng; bits += 8)
182 ulng = (ulng << 8) | prng_get_octet ();
183 return ulng;
184 }
185
186 /* Returns a pseudo-random long in the range [0, LONG_MAX]. */
187 long
prng_get_long(void)188 prng_get_long (void)
189 {
190 return prng_get_ulong () & LONG_MAX;
191 }
192
193 /* Returns a pseudo-random unsigned int in the range [0,
194 UINT_MAX]. */
195 unsigned
prng_get_uint(void)196 prng_get_uint (void)
197 {
198 unsigned uint;
199 size_t bits;
200
201 uint = prng_get_octet ();
202 for (bits = 8; bits < CHAR_BIT * sizeof uint; bits += 8)
203 uint = (uint << 8) | prng_get_octet ();
204 return uint;
205 }
206
207 /* Returns a pseudo-random int in the range [0, INT_MAX]. */
208 int
prng_get_int(void)209 prng_get_int (void)
210 {
211 return prng_get_uint () & INT_MAX;
212 }
213
214 /* Returns a pseudo-random floating-point number from the uniform
215 distribution with range [0,1). */
216 double
prng_get_double(void)217 prng_get_double (void)
218 {
219 for (;;)
220 {
221 double dbl = prng_get_ulong () / (ULONG_MAX + 1.0);
222 if (dbl >= 0.0 && dbl < 1.0)
223 return dbl;
224 }
225 }
226
227 /* Returns a pseudo-random floating-point number from the
228 distribution with mean 0 and standard deviation 1. (Multiply
229 the result by the desired standard deviation, then add the
230 desired mean.) */
231 double
prng_get_double_normal(void)232 prng_get_double_normal (void)
233 {
234 /* Knuth, _The Art of Computer Programming_, Vol. 2, 3.4.1C,
235 Algorithm P. */
236 static int has_next = 0;
237 static double next_normal;
238 double this_normal;
239
240 if (has_next)
241 {
242 this_normal = next_normal;
243 has_next = 0;
244 }
245 else
246 {
247 static double limit;
248 double v1, v2, s;
249
250 if (limit == 0.0)
251 limit = log (DBL_MAX / 2) / (DBL_MAX / 2);
252
253 for (;;)
254 {
255 double u1 = prng_get_double ();
256 double u2 = prng_get_double ();
257 v1 = 2.0 * u1 - 1.0;
258 v2 = 2.0 * u2 - 1.0;
259 s = v1 * v1 + v2 * v2;
260 if (s > limit && s < 1)
261 break;
262 }
263
264 this_normal = v1 * sqrt (-2. * log (s) / s);
265 next_normal = v2 * sqrt (-2. * log (s) / s);
266 has_next = 1;
267 }
268
269 return this_normal;
270 }
271