1 /* PSPP - computes sample statistics.
2 Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
3 Written by Ben Pfaff <blp@gnu.org>.
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.
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.
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
32 /* Random number generator. */
35 /* RC4-based random bytes. */
39 /* Normal distribution. */
44 /* Return a `random' seed by using the real time clock */
52 return (unsigned long) t;
55 /* Creates a new random number generator, seeds it based on
56 the current time, and returns it. */
61 static unsigned long seed=0;
64 rng = xmalloc (sizeof *rng);
67 if ( seed_is_set(&s) )
77 if (t == 0 || set_seed_used)
79 if (set_seed == NOT_LONG)
88 rng_seed (rng, &seed, sizeof seed);
89 rng->next_normal = NOT_DOUBLE;
95 rng_destroy (struct rng *rng)
102 swap_byte (uint8_t *a, uint8_t *b)
109 /* Seeds RNG based on the SIZE bytes in BUF.
110 At most the first 256 bytes of BUF are used. */
112 rng_seed (struct rng *rng, const void *key_, size_t size)
114 const uint8_t *key = key_;
119 assert (rng != NULL);
123 for (i = 0; i < 256; i++)
125 for (key_idx = 0, i = j = 0; i < 256; i++)
127 j = (j + s[i] + key[key_idx]) & 255;
128 swap_byte (s + i, s + j);
129 if (++key_idx >= size)
134 /* Reads SIZE random bytes from RNG into BUF. */
136 rng_get_bytes (struct rng *rng, void *buf_, size_t size)
151 swap_byte (s + i, s + j);
152 *buf++ = s[(s[i] + s[j]) & 255];
158 /* Returns a random int in the range [0, INT_MAX]. */
160 rng_get_int (struct rng *rng)
166 rng_get_bytes (rng, &value, sizeof value);
174 /* Returns a random unsigned in the range [0, UINT_MAX]. */
176 rng_get_unsigned (struct rng *rng)
180 rng_get_bytes (rng, &value, sizeof value);
184 /* Returns a random number from the uniform distribution with
187 rng_get_double (struct rng *rng)
194 rng_get_bytes (rng, &ulng, sizeof ulng);
195 dbl = ulng / (ULONG_MAX + 1.0);
196 if (dbl >= 0 && dbl < 1)
201 /* Returns a random number from the distribution with mean 0 and
202 standard deviation 1. (Multiply the result by the desired
203 standard deviation, then add the desired mean.) */
205 rng_get_double_normal (struct rng *rng)
207 /* Knuth, _The Art of Computer Programming_, Vol. 2, 3.4.1C,
211 if (rng->next_normal != NOT_DOUBLE)
213 this_normal = rng->next_normal;
214 rng->next_normal = NOT_DOUBLE;
222 double u1 = rng_get_double (rng);
223 double u2 = rng_get_double (rng);
226 s = v1 * v1 + v2 * v2;
230 this_normal = v1 * sqrt (-2. * log (s) / s);
231 rng->next_normal = v2 * sqrt (-2. * log (s) / s);
237 /* Gets an initialized RNG for use in PSPP transformations and
242 static struct rng *rng;