1 /*
2 * Copyright (c) 1996 Otmar Lendl (lendl@cosy.sbg.ac.at)
3 * All rights reserved.
4 *
5 * This file contains code written by others. Any notices of ownership
6 * have been retained.
7 *
8 * Redistribution and use in source and binary forms, with or without
9 * modification, are permitted without a fee provided that the following
10 * conditions are met:
11 *
12 * 1. This software is only used for private, research, or academic
13 * purposes.
14 *
15 * 2. Redistributions of source code must retain the above copyright
16 * notice, this list of conditions and the following disclaimer.
17 *
18 * 3. Redistributions in binary form must reproduce the above copyright
19 * notice, this list of conditions and the following disclaimer in the
20 * documentation and/or other materials provided with the distribution.
21 *
22 * 4. Any changes made to this package must be submitted to the author.
23 * The legal status of the submitted changes must allow their inclusion
24 * into this package under this license.
25 *
26 * 5. Publications in the field of pseudorandom number generation, which
27 * made use of this package must include a reference to this package.
28 *
29 * Any use of this software in a commercial environment requires a
30 * written licence from the author. Contact Otmar Lendl
31 * (lendl@cosy.sbg.ac.at) to negotiate the terms.
32 *
33 * THIS PACKAGE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS OR
34 * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
35 * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
36 *
37 * IN NO EVENT SHALL OTMAR LENDL BE LIABLE FOR ANY SPECIAL, INCIDENTAL,
38 * INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY KIND, OR ANY DAMAGES
39 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER OR
40 * NOT ADVISED OF THE POSSIBILITY OF DAMAGE, AND ON ANY THEORY OF
41 * LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE
42 * OF THIS SOFTWARE.
43 *
44 */
45 /*
46 *
47 * external.c External generator library
48 *
49 * Author: Otmar Lendl (lendl@cosy.sbg.ac.at)
50 *
51 * Last Modification: Sat Jan 11 21:28:14 MET 1997
52 *
53 */
54
55 /* Modification History:
56
57 21.11.96: First version: tt800
58 28.12.96: redesign
59 11.01.96: CTG/MRG/CMRG
60 17.10.00: use GNU automake and autoconf
61 01/11/07: replaced prng.free by prng.destroy
62 */
63
64 #include "prng.h"
65 #include <config.h>
66
67 #include <stdio.h>
68 #include <stdlib.h>
69 #include <string.h>
70 #include <errno.h>
71
72 /*
73 * We rely on int being 32 bits.
74 *
75 * for tt800 this is easy to fix, but L'Ecuyer's code depends heavily on it.
76 *
77 */
78
79 #if UINT_MAX != 4294967295U
80 #error "Wordsize Problems: int is not 32 bits."
81 #endif
82
83 /******************************************************************
84
85 TT800
86
87 *******************************************************************/
88
89 #define TT800_N 25
90 #define TT800_M 7
91 #define TT800_INV_MOD 2.3283064370807974e-10 /* 1.0 / (2^32-1) */
92
93 struct tt800_state
94 {
95 unsigned int x[TT800_N];
96 int k;
97 };
98
99 /*
100
101 Original Code by M. Matsumoto, taken from
102 ftp://random.mat.sbg.ac.at/pub/data/tt800.c, enhanced (reset, more than
103 one stream) and integrated into the libprng framework by Otmar Lendl.
104
105 */
106
107 /* A C-program for TT800 : July 8th 1996 Version */
108 /* by M. Matsumoto, email: matumoto@math.keio.ac.jp */
109 /* genrand() generate one pseudorandom number with double precision */
110 /* which is uniformly distributed on [0,1]-interval */
111 /* for each call. One may choose any initial 25 seeds */
112 /* except all zeros. */
113
114 /* See: ACM Transactions on Modelling and Computer Simulation, */
115 /* Vol. 4, No. 3, 1994, pages 254-266. */
116
117 /* we need 32 bits ore more for this numbers. 64 bits do not hurt. */
118
119 static struct tt800_state tt800_initial_state = {
120 { /* initial 25 seeds, */
121 0x95f24dab, 0x0b685215, 0xe76ccae7, 0xaf3ec239, 0x715fad23,
122 0x24a590ad, 0x69e4b5ef, 0xbf456141, 0x96bc1b7b, 0xa7bdf825,
123 0xc1de75b7, 0x8858a9c9, 0x2da87693, 0xb657f9dd, 0xffdc8a9f,
124 0x8121da71, 0x8b823ecb, 0x885d05f5, 0x4e20cd47, 0x5a9ad5d9,
125 0x512c0c03, 0xea857ccd, 0x4cc1d30f, 0x8891a8a1, 0xa6b7aadb
126 },
127 0 /* initial k */
128 };
129
130 static unsigned int tt800_mag01[2]=
131 { /* this is magic vector `a', don't change */
132 0x0, 0x8ebfd028
133 };
134
135 /*
136 * prng_tt800_get_next_int: Return next TT800 number (unscaled)
137 *
138 * Input:
139 * gen: Pointer to a struct prng.
140 *
141 */
prng_tt800_get_next_int(struct prng * gen)142 inline prng_num prng_tt800_get_next_int(struct prng *gen)
143 {
144 unsigned int y;
145 struct tt800_state *g;
146
147 g = (struct tt800_state *) gen->data.external_data.state;
148
149 if (g->k == TT800_N) /* generate TT800_N words at one time */
150 {
151 int kk;
152
153 for (kk=0; kk < TT800_N - TT800_M; kk++)
154 {
155 g->x[kk] = g->x[kk+TT800_M] ^
156 (g->x[kk] >> 1) ^ tt800_mag01[g->x[kk] % 2];
157 }
158
159 for (; kk<TT800_N;kk++)
160 {
161 g->x[kk] = g->x[kk+(TT800_M-TT800_N)] ^
162 (g->x[kk] >> 1) ^ tt800_mag01[g->x[kk] % 2];
163 }
164 g->k=0;
165 }
166
167 y = g->x[g->k];
168 y ^= (y << 7) & 0x2b5b2500; /* s and b, magic vectors */
169 y ^= (y << 15) & 0xdb8b0000; /* t and c, magic vectors */
170 /* y &= 0xffffffff; you may delete this line if word size = 32 */
171
172 /*
173 the following line was added by Makoto Matsumoto in the 1996 version
174 to improve lower bit's corellation.
175 Delete this line to o use the code published in 1994.
176 */
177
178 g->k++;
179
180 y ^= (y >> 16); /* added to the 1994 version */
181
182 return(y);
183 }
184
185 /*
186 * prng_tt800_get_next: Get next (scaled) TT800 number
187 *
188 * Input:
189 * gen: Pointer to a struct prng.
190 *
191 */
prng_tt800_get_next(struct prng * gen)192 inline double prng_tt800_get_next(struct prng *gen)
193 {
194 return(prng_tt800_get_next_int(gen) * TT800_INV_MOD);
195 }
196
197 /*-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=*/
198
199 /******************************************************************
200
201 CTG
202
203 *******************************************************************/
204
205 /*
206 * This generator was proposed by Pierre L'Ecuyer in 1995 in an
207 * article titled "Maximally equidistributed combined Tausworthe generators".
208 * (Not published yet, all I have is a preprint.)
209 *
210 * The following code is an adaption of the code printed there.
211 *
212 */
213
214 struct ctg_state
215 {
216 unsigned int s1,s2,s3;
217 };
218
219 /*
220 * prng_ctg_get_next_int: Return next ctg number (unscaled)
221 *
222 * Input:
223 * gen: Pointer to a struct prng.
224 *
225 */
prng_ctg_get_next_int(struct prng * gen)226 inline prng_num prng_ctg_get_next_int(struct prng *gen)
227 {
228 unsigned int b;
229 struct ctg_state *g;
230
231 g = (struct ctg_state *) gen->data.external_data.state;
232
233 b = (((g->s1 << 13) ^ g->s1) >> 19);
234 g->s1 = (((g->s1 & 4294967294U) << 12) ^ b);
235 b = (((g->s2 << 2) ^ g->s2) >> 25);
236 g->s2 = (((g->s2 & 4294967288U) << 4) ^ b);
237 b = (((g->s3 << 3) ^ g->s3) >> 11);
238 g->s3 = (((g->s3 & 4294967280U) << 17) ^ b);
239
240 return(g->s1 ^ g->s2 ^ g->s3);
241 }
242
243 /*
244 * prng_ctg_get_next: Get next (scaled) CTG number
245 *
246 * Input:
247 * gen: Pointer to a struct prng.
248 *
249 */
prng_ctg_get_next(struct prng * gen)250 inline double prng_ctg_get_next(struct prng *gen)
251 {
252 return(prng_ctg_get_next_int(gen) * 2.3283064365e-10);
253 }
254
255 /*-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=*/
256
257
258 /******************************************************************
259
260 MRG
261
262 *******************************************************************/
263
264 /*
265 * This generator was proposed by L'Ecuyer, Blouin, and Couture.
266 *
267 * See:
268 * "A search for good multiple recursive generators", ACM Trans. Model.
269 * Comput. Simul. 3:87-98 (1993)
270 *
271 * The following code is an adaption of the code printed there, which is
272 * optimized for 32 bit computers.
273 *
274 */
275
276 struct mrg_state
277 {
278 int x1,x2,x3,x4,x5;
279 };
280
281 /*
282 * prng_mrg_get_next_int: Return next mrg number (unscaled)
283 *
284 * Input:
285 * gen: Pointer to a struct prng.
286 *
287 */
prng_mrg_get_next_int(struct prng * gen)288 inline prng_num prng_mrg_get_next_int(struct prng *gen)
289 {
290 int h,p1,p5;
291 struct mrg_state *g;
292
293 g = (struct mrg_state *) gen->data.external_data.state;
294
295 h = g->x5 / 20554; p5 = 104480 * (g->x5 - h * 20554) - h * 1727;
296 g->x5 = g->x4; g->x4 = g->x3; g->x3 = g->x2; g->x2 = g->x1;
297 h = g->x1 / 20; p1 = 107374182 * (g->x1 - h * 20) - h * 7;
298
299 if (p1 < 0) p1 = p1 + 2147483647U;
300 if (p5 > 0) p5 = p5 - 2147483647U;
301
302 g->x1 = p1 + p5;
303 if (g->x1 < 0) g->x1 = g->x1 + 2147483647U;
304
305 return(g->x1);
306 }
307
308 /*
309 * prng_mrg_get_next: Get next (scaled) MRG number
310 *
311 * Input:
312 * gen: Pointer to a struct prng.
313 *
314 */
prng_mrg_get_next(struct prng * gen)315 inline double prng_mrg_get_next(struct prng *gen)
316 {
317 return(prng_mrg_get_next_int(gen) * 4.656612873077393e-10);
318 }
319
320 /*-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=*/
321
322
323
324 /******************************************************************
325
326 MRG
327
328 *******************************************************************/
329
330 /*
331 * This generator was proposed by Pierre L'Ecuyer.
332 *
333 * See:
334 * "Combined multiple recursive random number generators"
335 * Operations Research, 44, 5, 1996
336 *
337 * The following code is an adaption of the code printed in this paper.
338 *
339 */
340
341 struct cmrg_state
342 {
343 int x10, x11, x12, x20, x21, x22;
344 };
345
346 /*
347 * prng_cmrg_get_next_int: Return next cmrg number (unscaled)
348 *
349 * Input:
350 * gen: Pointer to a struct prng.
351 *
352 */
prng_cmrg_get_next_int(struct prng * gen)353 inline prng_num prng_cmrg_get_next_int(struct prng *gen)
354 {
355 int h, p12, p13, p21, p23;
356 struct cmrg_state *g;
357
358 g = (struct cmrg_state *) gen->data.external_data.state;
359
360 /* Component 1 */
361 h = g->x10 / 11714; p13 = 183326 * (g->x10 - h * 11714) - h * 2883;
362 h = g->x11 / 33921; p12 = 63308 * (g->x11 - h * 33921) - h * 12979;
363 if (p13 < 0) p13 = p13 + 2147483647U; if (p12 < 0) p12 = p12 + 2147483647U;
364 g->x10 = g->x11; g->x11 = g->x12; g->x12 = p12 - p13;
365 if (g->x12 < 0) g->x12 = g->x12+2147483647U;
366
367 /* Component 2 */
368 h = g->x20 / 3976; p23 = 539608 * (g->x20 - h * 3976) - h * 2071;
369 h = g->x22 / 24919; p21 = 86098 * (g->x22 - h * 24919) - h * 7417;
370 if (p23 < 0) p23 = p23 + 2145483479; if (p21 < 0) p21 = p21 + 2145483479U;
371 g->x20 = g->x21; g->x21 = g->x22; g->x22 = p21 - p23;
372 if (g->x22 < 0) g->x22 = g->x22+2145483479U;
373
374 /* Combination */
375 if (g->x12 < g->x22)
376 return(g->x12 - g->x22 + 2147483647U);
377 else
378 return(g->x12 - g->x22);
379 }
380
381 /*
382 * prng_cmrg_get_next: Get next (scaled) CMRG number
383 *
384 * Input:
385 * gen: Pointer to a struct prng.
386 *
387 */
prng_cmrg_get_next(struct prng * gen)388 inline double prng_cmrg_get_next(struct prng *gen)
389 {
390 return(prng_cmrg_get_next_int(gen) * 4.656612873077393e-10);
391 }
392
393 /*-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=*/
394
395
396 /*
397 * prng_external_get_array: Fill an Array with numbers, scaled to [0,1)
398 *
399 * Input:
400 * gen: Pointer to a struct prng.
401 *
402 * array: Where to put the numbers. ( double *array)
403 * n: Size of the array
404 *
405 */
prng_external_get_array(struct prng * gen,double * array,int n)406 void prng_external_get_array(struct prng *gen,double *array,int n)
407 {
408 int i;
409
410 for(i=0;i<n;array[i++] = prng_get_next(gen));
411 }
412
413
414 /*
415 * prng_external_reset: Reset an external generator.
416 *
417 * Input:
418 * gen: Pointer to a struct prng.
419 *
420 */
prng_external_reset(struct prng * gen)421 void prng_external_reset(struct prng *gen)
422 {
423 bcopy( (char *)gen->data.external_data.initial_state,
424 (char *)gen->data.external_data.state,
425 gen->data.external_data.state_size);
426 }
427
428 /*
429 * prng_external_free: Free all memory allocated by a PRNG
430 *
431 * Input:
432 * gen: Pointer to a struct prng.
433 *
434 */
prng_external_free(struct prng * gen)435 void prng_external_free(struct prng *gen)
436 {
437 if (gen->data.external_data.state != NULL)
438 free(gen->data.external_data.state);
439 prng_generic_free(gen);
440 }
441
442
443 /*
444 * Initialize External generator.
445 *
446 * Input:
447 *
448 * def: Pointer to struct prng_definition, as returned by
449 * prng_split_def.
450 *
451 * Returncode:
452 *
453 * struct prng * on success, NULL else
454 */
455
prng_external_init(struct prng_definition * def)456 struct prng *prng_external_init(struct prng_definition *def)
457 {
458 struct prng *gen;
459
460 if (strcasecmp("external",def->type) != 0 ) /* hmm. type seems to be wrong. */
461 {
462 return(NULL);
463 }
464
465 if (def->num_parameters < 1) /* right number of parameters */
466 {
467 return(NULL);
468 }
469
470 gen = prng_allocate();
471
472 /* default values */
473
474 gen->reset = prng_external_reset;
475 gen->destroy = prng_external_free;
476 gen->get_array = prng_external_get_array;
477 gen->is_congruential = FALSE;
478 gen->get_next_int = (prng_num (*)(struct prng *)) prng_undefined;
479 gen->modulus = 0;
480 gen->can_seed = FALSE;
481 gen->seed = (void (*)(struct prng *,prng_num seed)) prng_undefined;
482 gen->can_fast_sub = FALSE;
483 gen->get_sub_def = (char *(*)(struct prng *,prng_num s, prng_num i))
484 prng_undefined;
485 gen->can_fast_con = FALSE;
486 gen->get_con_def = (char *(*)(struct prng *,prng_num l, prng_num i))
487 prng_undefined;
488
489 /*************************** Now prepare the generator itself */
490
491 if (strcasecmp("tt800",def->parameter[0]) == 0 ) /* ctg ? */
492 {
493 if (def->num_parameters != 1) /* tt800 takes no parameters */
494 {
495 free(gen);
496 return(NULL);
497 }
498
499 gen->data.external_data.state_size = sizeof(tt800_initial_state);
500 gen->data.external_data.state =
501 malloc(gen->data.external_data.state_size);
502 if (gen->data.external_data.state == NULL)
503 {
504 fprintf(stderr,"Out of Memory.");
505 free(gen);
506 return(NULL);
507 }
508 gen->data.external_data.initial_state = &tt800_initial_state;
509
510 gen->get_next = prng_tt800_get_next;
511 gen->is_congruential = TRUE;
512 gen->get_next_int = prng_tt800_get_next_int;
513 gen->modulus = 0xffffffff;
514 gen->long_name = (char *) malloc(30);
515 if (gen->long_name != NULL)
516 strcpy(gen->long_name,"external(tt800)");
517 }
518
519 else if (strcasecmp("ctg",def->parameter[0]) == 0 ) /* CTG ? */
520 {
521 struct ctg_state *states;
522
523 if (def->num_parameters != 4) /* ctg takes 3 parameters */
524 {
525 free(gen);
526 return(NULL);
527 }
528
529 gen->data.external_data.state_size = sizeof(struct ctg_state);
530 states = (struct ctg_state *)
531 calloc(2,gen->data.external_data.state_size);
532
533 /* We allocate both initial state and state with one alloc.
534 state MUST be the first one, as prng_external_free will
535 call free() with that adress, so both structs will be
536 free with one call, too. */
537
538 if (states == NULL)
539 {
540 fprintf(stderr,"Out of Memory.");
541 free(gen);
542 return(NULL);
543 }
544
545 errno = 0;
546 states[1].s1 = strtoprng_num(def->parameter[1]);
547 states[1].s2 = strtoprng_num(def->parameter[2]);
548 states[1].s3 = strtoprng_num(def->parameter[3]);
549
550 if (errno != 0) /* errors while converting the numbers .. */
551 {
552 return(NULL);
553 }
554
555 gen->data.external_data.initial_state = &(states[1]);
556 gen->data.external_data.state = states;
557
558 gen->get_next = prng_ctg_get_next;
559 gen->is_congruential = TRUE;
560 gen->get_next_int = prng_ctg_get_next_int;
561 gen->modulus = 0xffffffff;
562
563 gen->long_name = (char *) malloc(3*PRNG_MAX_NUMBER_LEN + 15);
564 if (gen->long_name != NULL)
565 {
566 /* snprintf would be better, but it's not ubiquitous :( */
567 sprintf(gen->long_name,"external(ctg,%u,%u,%u)",
568 states[1].s1,states[1].s2,states[1].s3);
569 }
570 }
571
572 else if (strcasecmp("mrg",def->parameter[0]) == 0 ) /* MRG ? */
573 {
574 struct mrg_state *states;
575
576 if (def->num_parameters != 6) /* mrg takes 5 parameters */
577 {
578 free(gen);
579 return(NULL);
580 }
581
582 gen->data.external_data.state_size = sizeof(struct mrg_state);
583 states = (struct mrg_state *)
584 calloc(2,gen->data.external_data.state_size);
585
586 /* We allocate both initial state and state with one alloc.
587 state MUST be the first one, as prng_external_free will
588 call free() with that adress, so both structs will be
589 free with one call, too. */
590
591 if (states == NULL)
592 {
593 fprintf(stderr,"Out of Memory.");
594 free(gen);
595 return(NULL);
596 }
597
598 errno = 0;
599 states[1].x1 = strtoprng_num(def->parameter[1]);
600 states[1].x2 = strtoprng_num(def->parameter[2]);
601 states[1].x3 = strtoprng_num(def->parameter[3]);
602 states[1].x4 = strtoprng_num(def->parameter[4]);
603 states[1].x5 = strtoprng_num(def->parameter[5]);
604
605 if (errno != 0) /* errors while converting the numbers .. */
606 {
607 return(NULL);
608 }
609
610 gen->data.external_data.initial_state = &(states[1]);
611 gen->data.external_data.state = states;
612
613 gen->get_next = prng_mrg_get_next;
614 gen->is_congruential = TRUE;
615 gen->get_next_int = prng_mrg_get_next_int;
616 gen->modulus = 2147483647U;
617
618 gen->long_name = (char *) malloc(5*PRNG_MAX_NUMBER_LEN + 15);
619 if (gen->long_name != NULL)
620 {
621 /* snprintf would be better, but it's not ubiquitous :( */
622 sprintf(gen->long_name,"external(mrg,%u,%u,%u,%u,%u)",
623 states[1].x1,states[1].x2,states[1].x3,
624 states[1].x4,states[1].x5);
625 }
626 }
627
628 else if (strcasecmp("cmrg",def->parameter[0]) == 0 ) /* CMRG ? */
629 {
630 struct cmrg_state *states;
631
632 if (def->num_parameters != 7) /* cmrg takes 6 parameters */
633 {
634 free(gen);
635 return(NULL);
636 }
637
638 gen->data.external_data.state_size = sizeof(struct cmrg_state);
639 states = (struct cmrg_state *)
640 calloc(2,gen->data.external_data.state_size);
641
642 /* We allocate both initial state and state with one alloc.
643 state MUST be the first one, as prng_external_free will
644 call free() with that adress, so both structs will be
645 free with one call, too. */
646
647 if (states == NULL)
648 {
649 fprintf(stderr,"Out of Memory.");
650 free(gen);
651 return(NULL);
652 }
653
654 errno = 0;
655 states[1].x10 = strtoprng_num(def->parameter[1]);
656 states[1].x11 = strtoprng_num(def->parameter[2]);
657 states[1].x12 = strtoprng_num(def->parameter[3]);
658 states[1].x20 = strtoprng_num(def->parameter[4]);
659 states[1].x21 = strtoprng_num(def->parameter[5]);
660 states[1].x22 = strtoprng_num(def->parameter[6]);
661
662 if (errno != 0) /* errors while converting the numbers .. */
663 {
664 return(NULL);
665 }
666
667 gen->data.external_data.initial_state = &(states[1]);
668 gen->data.external_data.state = states;
669
670 gen->get_next = prng_cmrg_get_next;
671 gen->is_congruential = TRUE;
672 gen->get_next_int = prng_cmrg_get_next_int;
673 gen->modulus = 2147483647U;
674
675 gen->long_name = (char *) malloc(6*PRNG_MAX_NUMBER_LEN + 15);
676 if (gen->long_name != NULL)
677 {
678 /* snprintf would be better, but it's not ubiquitous :( */
679 sprintf(gen->long_name,"external(cmrg,%u,%u,%u,%u,%u,%u)",
680 states[1].x10,states[1].x11,states[1].x12,
681 states[1].x20,states[1].x21,states[1].x22);
682 }
683 }
684 else
685 {
686 free(gen);
687 return(NULL);
688 }
689
690 prng_reset(gen);
691
692 return(gen);
693 }
694
695