1 /* rng/random.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 #include "gsl__config.h"
21 #include <stdlib.h>
22 #include "gsl_rng.h"
23
24 /* This file provides support for random() generators. There are three
25 versions in widespread use today,
26
27 - The original BSD version, e.g. on SunOS 4.1 and FreeBSD.
28
29 - The Linux libc5 version, which is differs from the BSD version in
30 its seeding procedure, possibly due to the introduction of a typo
31 in the multiplier.
32
33 - The GNU glibc2 version, which has a new (and better) seeding
34 procedure.
35
36 They all produce different numbers, due to the different seeding
37 algorithms, but the algorithm for the generator is the same in each
38 case.
39
40 */
41
42 static inline long int random_get (int * i, int * j, int n, long int * x);
43
44 static inline unsigned long int random8_get (void *vstate);
45 static inline unsigned long int random32_get (void *vstate);
46 static inline unsigned long int random64_get (void *vstate);
47 static inline unsigned long int random128_get (void *vstate);
48 static inline unsigned long int random256_get (void *vstate);
49
50 static double random8_get_double (void *vstate);
51 static double random32_get_double (void *vstate);
52 static double random64_get_double (void *vstate);
53 static double random128_get_double (void *vstate);
54 static double random256_get_double (void *vstate);
55
56 static void random8_glibc2_set (void *state, unsigned long int s);
57 static void random32_glibc2_set (void *state, unsigned long int s);
58 static void random64_glibc2_set (void *state, unsigned long int s);
59 static void random128_glibc2_set (void *state, unsigned long int s);
60 static void random256_glibc2_set (void *state, unsigned long int s);
61
62 static void random8_libc5_set (void *state, unsigned long int s);
63 static void random32_libc5_set (void *state, unsigned long int s);
64 static void random64_libc5_set (void *state, unsigned long int s);
65 static void random128_libc5_set (void *state, unsigned long int s);
66 static void random256_libc5_set (void *state, unsigned long int s);
67
68 static void random8_bsd_set (void *state, unsigned long int s);
69 static void random32_bsd_set (void *state, unsigned long int s);
70 static void random64_bsd_set (void *state, unsigned long int s);
71 static void random128_bsd_set (void *state, unsigned long int s);
72 static void random256_bsd_set (void *state, unsigned long int s);
73
74 static void bsd_initialize (long int * x, int n, unsigned long int s);
75 static void libc5_initialize (long int * x, int n, unsigned long int s);
76 static void glibc2_initialize (long int * x, int n, unsigned long int s);
77
78 typedef struct
79 {
80 long int x;
81 }
82 random8_state_t;
83
84 typedef struct
85 {
86 int i, j;
87 long int x[7];
88 }
89 random32_state_t;
90
91 typedef struct
92 {
93 int i, j;
94 long int x[15];
95 }
96 random64_state_t;
97
98 typedef struct
99 {
100 int i, j;
101 long int x[31];
102 }
103 random128_state_t;
104
105 typedef struct
106 {
107 int i, j;
108 long int x[63];
109 }
110 random256_state_t;
111
112 static inline unsigned long int
random8_get(void * vstate)113 random8_get (void *vstate)
114 {
115 random8_state_t *state = (random8_state_t *) vstate;
116
117 state->x = (1103515245 * state->x + 12345) & 0x7fffffffUL;
118 return state->x;
119 }
120
121 static inline long int
random_get(int * i,int * j,int n,long int * x)122 random_get (int * i, int * j, int n, long int * x)
123 {
124 long int k ;
125
126 x[*i] += x[*j] ;
127 k = (x[*i] >> 1) & 0x7FFFFFFF ;
128
129 (*i)++ ;
130 if (*i == n)
131 *i = 0 ;
132
133 (*j)++ ;
134 if (*j == n)
135 *j = 0 ;
136
137 return k ;
138 }
139
140 static inline unsigned long int
random32_get(void * vstate)141 random32_get (void *vstate)
142 {
143 random32_state_t *state = (random32_state_t *) vstate;
144 unsigned long int k = random_get (&state->i, &state->j, 7, state->x) ;
145 return k ;
146 }
147
148 static inline unsigned long int
random64_get(void * vstate)149 random64_get (void *vstate)
150 {
151 random64_state_t *state = (random64_state_t *) vstate;
152 long int k = random_get (&state->i, &state->j, 15, state->x) ;
153 return k ;
154 }
155
156 static inline unsigned long int
random128_get(void * vstate)157 random128_get (void *vstate)
158 {
159 random128_state_t *state = (random128_state_t *) vstate;
160 unsigned long int k = random_get (&state->i, &state->j, 31, state->x) ;
161 return k ;
162 }
163
164 static inline unsigned long int
random256_get(void * vstate)165 random256_get (void *vstate)
166 {
167 random256_state_t *state = (random256_state_t *) vstate;
168 long int k = random_get (&state->i, &state->j, 63, state->x) ;
169 return k ;
170 }
171
172 static double
random8_get_double(void * vstate)173 random8_get_double (void *vstate)
174 {
175 return random8_get (vstate) / 2147483648.0 ;
176 }
177
178 static double
random32_get_double(void * vstate)179 random32_get_double (void *vstate)
180 {
181 return random32_get (vstate) / 2147483648.0 ;
182 }
183
184 static double
random64_get_double(void * vstate)185 random64_get_double (void *vstate)
186 {
187 return random64_get (vstate) / 2147483648.0 ;
188 }
189
190 static double
random128_get_double(void * vstate)191 random128_get_double (void *vstate)
192 {
193 return random128_get (vstate) / 2147483648.0 ;
194 }
195
196 static double
random256_get_double(void * vstate)197 random256_get_double (void *vstate)
198 {
199 return random256_get (vstate) / 2147483648.0 ;
200 }
201
202 static void
random8_bsd_set(void * vstate,unsigned long int s)203 random8_bsd_set (void *vstate, unsigned long int s)
204 {
205 random8_state_t *state = (random8_state_t *) vstate;
206
207 if (s == 0)
208 s = 1;
209
210 state->x = s;
211 }
212
213 static void
random32_bsd_set(void * vstate,unsigned long int s)214 random32_bsd_set (void *vstate, unsigned long int s)
215 {
216 random32_state_t *state = (random32_state_t *) vstate;
217 int i;
218
219 bsd_initialize (state->x, 7, s) ;
220
221 state->i = 3;
222 state->j = 0;
223
224 for (i = 0 ; i < 10 * 7 ; i++)
225 random32_get (state) ;
226 }
227
228 static void
random64_bsd_set(void * vstate,unsigned long int s)229 random64_bsd_set (void *vstate, unsigned long int s)
230 {
231 random64_state_t *state = (random64_state_t *) vstate;
232 int i;
233
234 bsd_initialize (state->x, 15, s) ;
235
236 state->i = 1;
237 state->j = 0;
238
239 for (i = 0 ; i < 10 * 15 ; i++)
240 random64_get (state) ;
241 }
242
243 static void
random128_bsd_set(void * vstate,unsigned long int s)244 random128_bsd_set (void *vstate, unsigned long int s)
245 {
246 random128_state_t *state = (random128_state_t *) vstate;
247 int i;
248
249 bsd_initialize (state->x, 31, s) ;
250
251 state->i = 3;
252 state->j = 0;
253
254 for (i = 0 ; i < 10 * 31 ; i++)
255 random128_get (state) ;
256 }
257
258 static void
random256_bsd_set(void * vstate,unsigned long int s)259 random256_bsd_set (void *vstate, unsigned long int s)
260 {
261 random256_state_t *state = (random256_state_t *) vstate;
262 int i;
263
264 bsd_initialize (state->x, 63, s) ;
265
266 state->i = 1;
267 state->j = 0;
268
269 for (i = 0 ; i < 10 * 63 ; i++)
270 random256_get (state) ;
271 }
272
273 static void
bsd_initialize(long int * x,int n,unsigned long int s)274 bsd_initialize (long int * x, int n, unsigned long int s)
275 {
276 int i;
277
278 if (s == 0)
279 s = 1 ;
280
281 x[0] = s;
282
283 for (i = 1 ; i < n ; i++)
284 x[i] = 1103515245 * x[i-1] + 12345 ;
285 }
286
287 static void
libc5_initialize(long int * x,int n,unsigned long int s)288 libc5_initialize (long int * x, int n, unsigned long int s)
289 {
290 int i;
291
292 if (s == 0)
293 s = 1 ;
294
295 x[0] = s;
296
297 for (i = 1 ; i < n ; i++)
298 x[i] = 1103515145 * x[i-1] + 12345 ;
299 }
300
301 static void
glibc2_initialize(long int * x,int n,unsigned long int s)302 glibc2_initialize (long int * x, int n, unsigned long int s)
303 {
304 int i;
305
306 if (s == 0)
307 s = 1 ;
308
309 x[0] = s;
310
311 for (i = 1 ; i < n ; i++)
312 {
313 const long int h = s / 127773;
314 const long int t = 16807 * (s - h * 127773) - h * 2836;
315 if (t < 0)
316 {
317 s = t + 2147483647 ;
318 }
319 else
320 {
321 s = t ;
322 }
323
324 x[i] = s ;
325 }
326 }
327
328 static void
random8_glibc2_set(void * vstate,unsigned long int s)329 random8_glibc2_set (void *vstate, unsigned long int s)
330 {
331 random8_state_t *state = (random8_state_t *) vstate;
332
333 if (s == 0)
334 s = 1;
335
336 state->x = s;
337 }
338
339 static void
random32_glibc2_set(void * vstate,unsigned long int s)340 random32_glibc2_set (void *vstate, unsigned long int s)
341 {
342 random32_state_t *state = (random32_state_t *) vstate;
343 int i;
344
345 glibc2_initialize (state->x, 7, s) ;
346
347 state->i = 3;
348 state->j = 0;
349
350 for (i = 0 ; i < 10 * 7 ; i++)
351 random32_get (state) ;
352 }
353
354 static void
random64_glibc2_set(void * vstate,unsigned long int s)355 random64_glibc2_set (void *vstate, unsigned long int s)
356 {
357 random64_state_t *state = (random64_state_t *) vstate;
358 int i;
359
360 glibc2_initialize (state->x, 15, s) ;
361
362 state->i = 1;
363 state->j = 0;
364
365 for (i = 0 ; i < 10 * 15 ; i++)
366 random64_get (state) ;
367 }
368
369 static void
random128_glibc2_set(void * vstate,unsigned long int s)370 random128_glibc2_set (void *vstate, unsigned long int s)
371 {
372 random128_state_t *state = (random128_state_t *) vstate;
373 int i;
374
375 glibc2_initialize (state->x, 31, s) ;
376
377 state->i = 3;
378 state->j = 0;
379
380 for (i = 0 ; i < 10 * 31 ; i++)
381 random128_get (state) ;
382 }
383
384 static void
random256_glibc2_set(void * vstate,unsigned long int s)385 random256_glibc2_set (void *vstate, unsigned long int s)
386 {
387 random256_state_t *state = (random256_state_t *) vstate;
388 int i;
389
390 glibc2_initialize (state->x, 63, s) ;
391
392 state->i = 1;
393 state->j = 0;
394
395 for (i = 0 ; i < 10 * 63 ; i++)
396 random256_get (state) ;
397 }
398
399
400 static void
random8_libc5_set(void * vstate,unsigned long int s)401 random8_libc5_set (void *vstate, unsigned long int s)
402 {
403 random8_state_t *state = (random8_state_t *) vstate;
404
405 if (s == 0)
406 s = 1;
407
408 state->x = s;
409 }
410
411 static void
random32_libc5_set(void * vstate,unsigned long int s)412 random32_libc5_set (void *vstate, unsigned long int s)
413 {
414 random32_state_t *state = (random32_state_t *) vstate;
415 int i;
416
417 libc5_initialize (state->x, 7, s) ;
418
419 state->i = 3;
420 state->j = 0;
421
422 for (i = 0 ; i < 10 * 7 ; i++)
423 random32_get (state) ;
424 }
425
426 static void
random64_libc5_set(void * vstate,unsigned long int s)427 random64_libc5_set (void *vstate, unsigned long int s)
428 {
429 random64_state_t *state = (random64_state_t *) vstate;
430 int i;
431
432 libc5_initialize (state->x, 15, s) ;
433
434 state->i = 1;
435 state->j = 0;
436
437 for (i = 0 ; i < 10 * 15 ; i++)
438 random64_get (state) ;
439 }
440
441 static void
random128_libc5_set(void * vstate,unsigned long int s)442 random128_libc5_set (void *vstate, unsigned long int s)
443 {
444 random128_state_t *state = (random128_state_t *) vstate;
445 int i;
446
447 libc5_initialize (state->x, 31, s) ;
448
449 state->i = 3;
450 state->j = 0;
451
452 for (i = 0 ; i < 10 * 31 ; i++)
453 random128_get (state) ;
454 }
455
456 static void
random256_libc5_set(void * vstate,unsigned long int s)457 random256_libc5_set (void *vstate, unsigned long int s)
458 {
459 random256_state_t *state = (random256_state_t *) vstate;
460 int i;
461
462 libc5_initialize (state->x, 63, s) ;
463
464 state->i = 1;
465 state->j = 0;
466
467 for (i = 0 ; i < 10 * 63 ; i++)
468 random256_get (state) ;
469 }
470
471 static const gsl_rng_type random_glibc2_type =
472 {"random-glibc2", /* name */
473 0x7fffffffUL, /* RAND_MAX */
474 0, /* RAND_MIN */
475 sizeof (random128_state_t),
476 &random128_glibc2_set,
477 &random128_get,
478 &random128_get_double};
479
480 static const gsl_rng_type random8_glibc2_type =
481 {"random8-glibc2", /* name */
482 0x7fffffffUL, /* RAND_MAX */
483 0, /* RAND_MIN */
484 sizeof (random8_state_t),
485 &random8_glibc2_set,
486 &random8_get,
487 &random8_get_double};
488
489 static const gsl_rng_type random32_glibc2_type =
490 {"random32-glibc2", /* name */
491 0x7fffffffUL, /* RAND_MAX */
492 0, /* RAND_MIN */
493 sizeof (random32_state_t),
494 &random32_glibc2_set,
495 &random32_get,
496 &random32_get_double};
497
498 static const gsl_rng_type random64_glibc2_type =
499 {"random64-glibc2", /* name */
500 0x7fffffffUL, /* RAND_MAX */
501 0, /* RAND_MIN */
502 sizeof (random64_state_t),
503 &random64_glibc2_set,
504 &random64_get,
505 &random64_get_double};
506
507 static const gsl_rng_type random128_glibc2_type =
508 {"random128-glibc2", /* name */
509 0x7fffffffUL, /* RAND_MAX */
510 0, /* RAND_MIN */
511 sizeof (random128_state_t),
512 &random128_glibc2_set,
513 &random128_get,
514 &random128_get_double};
515
516 static const gsl_rng_type random256_glibc2_type =
517 {"random256-glibc2", /* name */
518 0x7fffffffUL, /* RAND_MAX */
519 0, /* RAND_MIN */
520 sizeof (random256_state_t),
521 &random256_glibc2_set,
522 &random256_get,
523 &random256_get_double};
524
525 static const gsl_rng_type random_libc5_type =
526 {"random-libc5", /* name */
527 0x7fffffffUL, /* RAND_MAX */
528 0, /* RAND_MIN */
529 sizeof (random128_state_t),
530 &random128_libc5_set,
531 &random128_get,
532 &random128_get_double};
533
534 static const gsl_rng_type random8_libc5_type =
535 {"random8-libc5", /* name */
536 0x7fffffffUL, /* RAND_MAX */
537 0, /* RAND_MIN */
538 sizeof (random8_state_t),
539 &random8_libc5_set,
540 &random8_get,
541 &random8_get_double};
542
543 static const gsl_rng_type random32_libc5_type =
544 {"random32-libc5", /* name */
545 0x7fffffffUL, /* RAND_MAX */
546 0, /* RAND_MIN */
547 sizeof (random32_state_t),
548 &random32_libc5_set,
549 &random32_get,
550 &random32_get_double};
551
552 static const gsl_rng_type random64_libc5_type =
553 {"random64-libc5", /* name */
554 0x7fffffffUL, /* RAND_MAX */
555 0, /* RAND_MIN */
556 sizeof (random64_state_t),
557 &random64_libc5_set,
558 &random64_get,
559 &random64_get_double};
560
561 static const gsl_rng_type random128_libc5_type =
562 {"random128-libc5", /* name */
563 0x7fffffffUL, /* RAND_MAX */
564 0, /* RAND_MIN */
565 sizeof (random128_state_t),
566 &random128_libc5_set,
567 &random128_get,
568 &random128_get_double};
569
570 static const gsl_rng_type random256_libc5_type =
571 {"random256-libc5", /* name */
572 0x7fffffffUL, /* RAND_MAX */
573 0, /* RAND_MIN */
574 sizeof (random256_state_t),
575 &random256_libc5_set,
576 &random256_get,
577 &random256_get_double};
578
579 static const gsl_rng_type random_bsd_type =
580 {"random-bsd", /* name */
581 0x7fffffffUL, /* RAND_MAX */
582 0, /* RAND_MIN */
583 sizeof (random128_state_t),
584 &random128_bsd_set,
585 &random128_get,
586 &random128_get_double};
587
588 static const gsl_rng_type random8_bsd_type =
589 {"random8-bsd", /* name */
590 0x7fffffffUL, /* RAND_MAX */
591 0, /* RAND_MIN */
592 sizeof (random8_state_t),
593 &random8_bsd_set,
594 &random8_get,
595 &random8_get_double};
596
597 static const gsl_rng_type random32_bsd_type =
598 {"random32-bsd", /* name */
599 0x7fffffffUL, /* RAND_MAX */
600 0, /* RAND_MIN */
601 sizeof (random32_state_t),
602 &random32_bsd_set,
603 &random32_get,
604 &random32_get_double};
605
606 static const gsl_rng_type random64_bsd_type =
607 {"random64-bsd", /* name */
608 0x7fffffffUL, /* RAND_MAX */
609 0, /* RAND_MIN */
610 sizeof (random64_state_t),
611 &random64_bsd_set,
612 &random64_get,
613 &random64_get_double};
614
615 static const gsl_rng_type random128_bsd_type =
616 {"random128-bsd", /* name */
617 0x7fffffffUL, /* RAND_MAX */
618 0, /* RAND_MIN */
619 sizeof (random128_state_t),
620 &random128_bsd_set,
621 &random128_get,
622 &random128_get_double};
623
624 static const gsl_rng_type random256_bsd_type =
625 {"random256-bsd", /* name */
626 0x7fffffffUL, /* RAND_MAX */
627 0, /* RAND_MIN */
628 sizeof (random256_state_t),
629 &random256_bsd_set,
630 &random256_get,
631 &random256_get_double};
632
633 const gsl_rng_type *gsl_rng_random_libc5 = &random_libc5_type;
634 const gsl_rng_type *gsl_rng_random8_libc5 = &random8_libc5_type;
635 const gsl_rng_type *gsl_rng_random32_libc5 = &random32_libc5_type;
636 const gsl_rng_type *gsl_rng_random64_libc5 = &random64_libc5_type;
637 const gsl_rng_type *gsl_rng_random128_libc5 = &random128_libc5_type;
638 const gsl_rng_type *gsl_rng_random256_libc5 = &random256_libc5_type;
639
640 const gsl_rng_type *gsl_rng_random_glibc2 = &random_glibc2_type;
641 const gsl_rng_type *gsl_rng_random8_glibc2 = &random8_glibc2_type;
642 const gsl_rng_type *gsl_rng_random32_glibc2 = &random32_glibc2_type;
643 const gsl_rng_type *gsl_rng_random64_glibc2 = &random64_glibc2_type;
644 const gsl_rng_type *gsl_rng_random128_glibc2 = &random128_glibc2_type;
645 const gsl_rng_type *gsl_rng_random256_glibc2 = &random256_glibc2_type;
646
647 const gsl_rng_type *gsl_rng_random_bsd = &random_bsd_type;
648 const gsl_rng_type *gsl_rng_random8_bsd = &random8_bsd_type;
649 const gsl_rng_type *gsl_rng_random32_bsd = &random32_bsd_type;
650 const gsl_rng_type *gsl_rng_random64_bsd = &random64_bsd_type;
651 const gsl_rng_type *gsl_rng_random128_bsd = &random128_bsd_type;
652 const gsl_rng_type *gsl_rng_random256_bsd = &random256_bsd_type;
653
654
655
656
657