1 /*************************************************************************\
2 *
3 * Package: TestU01
4 * File: ugfsr.c
5 * Environment: ANSI C
6 *
7 * Copyright (c) 2002 Pierre L'Ecuyer, DIRO, Université de Montréal.
8 * e-mail: lecuyer@iro.umontreal.ca
9 * All rights reserved.
10 *
11 * Redistribution and use in source and binary forms, with or without
12 * modification, are permitted without a fee for private, research,
13 * academic, or other non-commercial purposes.
14 * Any use of this software in a commercial environment requires a
15 * written licence from the copyright owner.
16 *
17 * Any changes made to this package must be clearly identified as such.
18 *
19 * In scientific publications which used this software, a reference to it
20 * would be appreciated.
21 *
22 * Redistributions of source code must retain this copyright notice
23 * and the following disclaimer.
24 *
25 * THIS PACKAGE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS OR
26 * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
27 * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
28 *
29 \*************************************************************************/
30
31 #include "util.h"
32 #include "num.h"
33 #include "addstr.h"
34
35 #include "ugfsr.h"
36 #include "unif01.h"
37
38 #include <stdio.h>
39 #include <string.h>
40
41
42
43
44
45 /*=============================== Macros ===============================*/
46
47 #define LEN 300 /* Max length of strings */
48
49 #define MATRIX_A 0x9908b0dfUL /* constant vector a for MT19937 */
50 #define UPPER_MASK 0x80000000UL /* most significant w - r bits */
51 #define LOWER_MASK 0x7fffffff /* least significant r bits */
52
53 /* for tempering */
54 #define TEMPERING_MASK_B 0x9d2c5680UL
55 #define TEMPERING_MASK_C 0xefc60000UL
56 #define TEMPERING_SHIFT_U(y) ((y) >> 11)
57 #define TEMPERING_SHIFT_S(y) ((y) << 7)
58 #define TEMPERING_SHIFT_T(y) ((y) << 15)
59 #define TEMPERING_SHIFT_L(y) ((y) >> 18)
60
61
62
63
64
65
66 /*================================ Types ================================*/
67
68 typedef struct {
69 unsigned long Shift, Mask, mag01[2];
70 double Norm;
71 } GFSR_param;
72
73 typedef struct {
74 unsigned long *X;
75 unsigned int r, s, K;
76 } GFSR_state;
77
78 /*------------------------------------*/
79
80 typedef struct {
81 unsigned long Shift, Mask, mag01[2];
82 unsigned long b, c, S2, T2;
83 } TGFSR2_param;
84
85 typedef GFSR_state TGFSR2_state;
86
87 /*------------------------------------*/
88
89 typedef struct {
90 unsigned long mag01[2];
91 } MT19937_param;
92
93 typedef GFSR_state MT19937_state;
94
95 /*------------------------------------*/
96
97 typedef struct {
98 unsigned long Shift;
99 } GFSR5_param;
100
101 typedef struct {
102 unsigned long *X;
103 unsigned int r1, r2, r3, s, K;
104 } GFSR5_state;
105
106 /*------------------------------------*/
107
108
109
110
111
112
113 /***************************************************************************/
114
PRODUIT(int k,unsigned int F[],unsigned int G[],unsigned int H[])115 static void PRODUIT (int k, unsigned int F[], unsigned int G[],
116 unsigned int H[])
117 {
118 int i, j;
119 unsigned int Work;
120
121 for (i = 0; i < k; i++) {
122 Work = 0;
123 for (j = 0; j < k; j++) {
124 if (G[j] == 1)
125 Work += F[i + j];
126 }
127 H[i] = Work % 2;
128 }
129 }
130
131
132 /***************************************************************************/
133
134 #define Fushia 69069 /* multiplier */
135
InitFushimi(int k,int r,int s,GFSR_state * state)136 static void InitFushimi (int k, int r, int s, GFSR_state *state)
137 /*
138 * Initialization copied from Fushimi90
139 */
140 {
141 unsigned int *x, *x0, *x1, *x2, *x3;
142 unsigned int E0, E1;
143 int wi, ix, B[32];
144 int i, iMD2, k2, k3, j;
145
146 E0 = 0;
147 E1 = 1;
148 k2 = k * 2;
149 k3 = k * 3;
150 state->K = k3;
151 state->r = 3 * (k - r);
152 state->s = 0;
153
154 x = calloc ((size_t) 3 * (k + 1), sizeof (unsigned int));
155 x0 = calloc ((size_t) 2 * (k + 1), sizeof (unsigned int));
156 x1 = calloc ((size_t) 2 * (k + 1), sizeof (unsigned int));
157 x2 = calloc ((size_t) 1 * (k + 1), sizeof (unsigned int));
158 x3 = calloc ((size_t) 3 * (k + 1), sizeof (unsigned int));
159
160 /* Initialization */
161 for (j = 0; j < k2; j++) {
162 x0[j] = E0;
163 x3[j] = E0;
164 }
165
166 /* Table B contains 2^30, 2^29, ... , 1, 0. */
167 B[31] = 0;
168 B[30] = 1;
169 for (j = 0; j <= 29; j++)
170 B[29 - j] = B[30 - j] * 2;
171
172 ix = s;
173 for (j = 0; j < k; j++) {
174 ix *= Fushia;
175 if (ix > 0)
176 x0[j] = E1;
177 }
178 for (j = k; j < k2; j++)
179 x0[j] = x0[j - k] ^ x0[j - r];
180 x3[1] = E1;
181 for (i = 0; i <= k - 2; i++) {
182 iMD2 = i % 2;
183 for (j = k - 1; j >= 0; j--) {
184 x3[2 * j + iMD2] = x3[j];
185 x3[2 * j + 1 - iMD2] = E0;
186 }
187 for (j = k - 1; j >= 0; j--) {
188 x3[j] ^= x3[j + k];
189 x3[j + k - r] ^= x3[j + k];
190 x3[j + k] = E0;
191 }
192 }
193 PRODUIT (k, x0, x3, x1);
194 for (j = k; j < k2; j++)
195 x1[j] = x1[j - k] ^ x1[j - r];
196 PRODUIT (k, x1, x3, x2);
197 for (j = 0; j <= k; j++) {
198 x[3 * j] = x0[j];
199 x[3 * j + 1] = x1[j];
200 x[3 * j + 2] = x2[j];
201 }
202 for (i = 0; i < k3; i++) {
203 wi = 0;
204 for (j = 0; j < 32; j++) {
205 if (x[state->s] != E0)
206 wi += B[j];
207 x[state->s] ^= x[state->r];
208 state->s++;
209 if (state->s == state->K)
210 state->s = 0;
211 state->r++;
212 if (state->r == state->K)
213 state->r = 0;
214 }
215 state->X[i] = wi & unif01_MASK32;
216 }
217 free (x);
218 free (x0);
219 free (x1);
220 free (x2);
221 free (x3);
222 }
223
224
225 /***************************************************************************/
226
GFSR_Bits(void * vpar,void * vsta)227 static unsigned long GFSR_Bits (void *vpar, void *vsta)
228 {
229 GFSR_param *param = vpar;
230 GFSR_state *state = vsta;
231
232 if (++state->s == state->K)
233 state->s = 0; /* s = (s + 1) % K */
234 if (++state->r == state->K)
235 state->r = 0; /* r = (r + 1) % K */
236 state->X[state->s] ^= state->X[state->r];
237 return state->X[state->s] << param->Shift;
238 }
239
240 /*-----------------------------------------------------------------------*/
241
GFSR_U01(void * vpar,void * vsta)242 static double GFSR_U01 (void *vpar, void *vsta)
243 {
244 return GFSR_Bits (vpar, vsta) * unif01_INV32;
245 }
246
247 /*-----------------------------------------------------------------------*/
248
WrGFSR(void * vsta)249 static void WrGFSR (void *vsta)
250 {
251 GFSR_state *state = vsta;
252 unsigned int j;
253 unsigned int s = state->s;
254
255 if (unif01_WrLongStateFlag) {
256 printf (" S = {\n ");
257 for (j = 0; j < state->K; j++) {
258 if (++s >= state->K)
259 s = 0;
260 printf (" %12lu", state->X[s]);
261 if (j < state->K - 1)
262 printf (",");
263 if ((j % 5) == 4)
264 printf ("\n ");
265 };
266 printf (" }\n");
267 } else
268 unif01_WrLongStateDef ();
269 }
270
271 /*-----------------------------------------------------------------------*/
272
CreateGFSR0(unsigned int k,unsigned int r,unsigned int l,unsigned long S[],const char nom[])273 static unif01_Gen *CreateGFSR0 (unsigned int k, unsigned int r,
274 unsigned int l, unsigned long S[], const char nom[])
275 /*
276 * CreateGen common to many generators GFSR and TGFSR
277 */
278 {
279 unif01_Gen *gen;
280 GFSR_param *param;
281 GFSR_state *state;
282 size_t leng;
283 char name[LEN + 1];
284 unsigned int j;
285 unsigned long Mask;
286
287 strcpy (name, nom);
288 addstr_Uint (name, " k = ", k);
289 addstr_Uint (name, ", r = ", r);
290 addstr_Uint (name, ", l = ", l);
291 addstr_ArrayUlong (name, ", S = ", (int) k, S);
292 if ((l < 1) || (r >= k) || (l > 32))
293 util_Error (name);
294
295 gen = util_Malloc (sizeof (unif01_Gen));
296 param = util_Malloc (sizeof (GFSR_param));
297 state = util_Malloc (sizeof (GFSR_state));
298
299 leng = strlen (name);
300 gen->name = util_Calloc (leng + 1, sizeof (char));
301 strncpy (gen->name, name, leng);
302
303 if (l == 32)
304 Mask = unif01_MASK32;
305 else
306 Mask = num_TwoExp[l] - 1.0;
307 state->X = util_Calloc ((size_t) k, sizeof (unsigned long));
308 for (j = 0; j < k; j++)
309 state->X[j] = S[j] & Mask;
310
311 state->s = 0;
312 state->r = k - r;
313 state->K = k;
314 param->Shift = 32 - l;
315 param->Mask = Mask;
316
317 gen->param = param;
318 gen->state = state;
319 gen->GetBits = &GFSR_Bits;
320 gen->GetU01 = &GFSR_U01;
321 gen->Write = &WrGFSR;
322 return gen;
323 }
324
325 /*=========================================================================*/
326
ugfsr_CreateGFSR3(unsigned int k,unsigned int r,unsigned int l,unsigned long S[])327 unif01_Gen *ugfsr_CreateGFSR3 (unsigned int k, unsigned int r,
328 unsigned int l, unsigned long S[])
329 {
330 return CreateGFSR0 (k, r, l, S, "ugfsr_CreateGFSR3:");
331 }
332
333
334 /*=========================================================================*/
335
336 #define RIPa 16807 /* multiplier */
337 #define RIPk 521 /* Parametres for */
338 #define RIPr 32 /* Ripley90 .... */
339 #define RIPk2 (2*RIPk)
340 #define RIPc 127773
341
342
Ripley90_U01(void * vpar,void * vsta)343 static double Ripley90_U01 (void *vpar, void *vsta)
344 {
345 GFSR_param *param = vpar;
346 GFSR_state *state = vsta;
347 unsigned long temp;
348
349 state->s--;
350 state->r--;
351 temp = state->X[state->r];
352 state->X[state->r] ^= state->X[state->s];
353 if (state->s == 0)
354 state->s = RIPk;
355 if (state->r == 0)
356 state->r = RIPk;
357 return temp * param->Norm;
358 }
359
360 /*-----------------------------------------------------------------------*/
361
Ripley90_Bits(void * vpar,void * vsta)362 static unsigned long Ripley90_Bits (void *vpar, void *vsta)
363 {
364 return (unsigned long) (Ripley90_U01 (vpar, vsta) * unif01_NORM32);
365 }
366
367 /*-----------------------------------------------------------------------*/
368
WrRipley90(void * vsta)369 static void WrRipley90 (void *vsta)
370 {
371 GFSR_state *state = vsta;
372 int j;
373 int r = state->r;
374
375 if (unif01_WrLongStateFlag) {
376 printf (" S = {\n ");
377 for (j = 0; j < RIPk; j++) {
378 r--;
379 printf (" %12lu", state->X[r]);
380 if (r <= 0)
381 r = RIPk;
382 if (j < RIPk - 1)
383 printf (",");
384 if ((j % 5) == 4)
385 printf ("\n ");
386 };
387 printf (" }\n");
388 } else
389 unif01_WrLongStateDef ();
390 }
391
392 /*-----------------------------------------------------------------------*/
393
ugfsr_CreateRipley90(long s)394 unif01_Gen * ugfsr_CreateRipley90 (long s)
395 {
396 unif01_Gen *gen;
397 GFSR_param *param;
398 GFSR_state *state;
399 size_t leng;
400 char name[LEN + 1];
401 unsigned long x0[RIPk2];
402 unsigned long wi, wt;
403 long vt, ix, i, j; /* Work variables */
404
405 gen = util_Malloc (sizeof (unif01_Gen));
406 param = util_Malloc (sizeof (GFSR_param));
407 state = util_Malloc (sizeof (GFSR_state));
408
409 strcpy (name, "ugfsr_CreateRipley90:");
410 addstr_Long (name, " s = ", s);
411 leng = strlen (name);
412 gen->name = util_Calloc (leng + 1, sizeof (char));
413 strncpy (gen->name, name, leng);
414
415 state->K = RIPk;
416 state->r = RIPk - RIPr;
417 state->s = RIPk;
418 param->Norm = 1.0 / (num_TwoExp[31] - 1.0);
419 state->X = util_Calloc ((size_t) RIPk, sizeof (unsigned long));
420
421 /* Initialization */
422 ix = s;
423 for (j = 0; j < RIPk; j++) {
424 x0[j] = 0;
425 vt = ix / RIPc;
426 ix = RIPa * (ix - RIPc * vt) - 2836 * vt;
427 if (ix < 0)
428 ix += 2147483647;
429 if (ix > 1073741824)
430 x0[j] = 1;
431 }
432 for (j = RIPk; j < RIPk2; j++)
433 x0[j] = x0[j - RIPk] ^ x0[j - RIPr];
434 for (i = 0; i < RIPk; i++) {
435 wi = 0;
436 for (j = 0; j <= 30; j++) {
437 wt = x0[i + 16 * (j + 1)];
438 wt <<= j;
439 #ifndef IS_ULONG32
440 wt &= unif01_MASK32;
441 #endif
442 wi += wt;
443 #ifndef IS_ULONG32
444 wi &= unif01_MASK32;
445 #endif
446 }
447 state->X[i] = wi & unif01_MASK32;
448 }
449 gen->param = param;
450 gen->state = state;
451 gen->GetBits = &Ripley90_Bits;
452 gen->GetU01 = &Ripley90_U01;
453 gen->Write = &WrRipley90;
454 return gen;
455 }
456
457
458 /***************************************************************************/
459
ugfsr_CreateToot73(unsigned long S[])460 unif01_Gen * ugfsr_CreateToot73 (unsigned long S[])
461 {
462 unif01_Gen *gen;
463 GFSR_param *param;
464 GFSR_state *state;
465 size_t leng;
466 char name[LEN + 1];
467 long q, l, m, n;
468 int h, t, j, i, K = 607, R = 273;
469 unsigned long Q[700], Mask;
470
471 gen = util_Malloc (sizeof (unif01_Gen));
472 param = util_Malloc (sizeof (GFSR_param));
473 state = util_Malloc (sizeof (GFSR_state));
474
475 strcpy (name, "ugfsr_CreateToot73:");
476 addstr_ArrayUlong (name, " S = ", K, S);
477 leng = strlen (name);
478 gen->name = util_Calloc (leng + 1, sizeof (char));
479 strncpy (gen->name, name, leng);
480
481 state->X = util_Calloc ((size_t) 700, sizeof (unsigned long));
482 state->r = K - R;
483 state->s = 0;
484 state->K = K;
485 Mask = num_TwoExp[23] - 1.0;
486 param->Shift = 9;
487 param->Mask = Mask;
488
489 q = S[19];
490 m = S[10];
491 for (j = 1; j <= 19; j++)
492 Q[j] = S[j];
493 for (j = 19; j <= K; j++) {
494 l = (long) Q[j - 18];
495 n = (long) Q[j - 8];
496 Q[j] = (((unsigned long) q) << 1) + (((unsigned long) l) >> 31);
497 Q[j] ^= (((unsigned long) m) << 15) + (((unsigned long) n) >> 17);
498 Q[j] &= unif01_MASK32;
499 q = l;
500 m = n;
501 Q[j - 18] = l & Mask;
502 }
503
504 for (j = 590; j <= K; j++)
505 Q[j] &= Mask;
506 i = 1;
507 t = 0;
508 do {
509 for (j = i; j <= K; j += 16) {
510 state->X[t] = Q[j];
511 t++;
512 }
513 i++;
514 /* Recompute table Q so that the equation */
515 /* Q(n) =( ( x^15*Q(n- 9)+x^- 17*Q(n- 8) ) XOR */
516 /* ( x*Q(n-19)+x^- 31*Q(n-18) ) */
517 /* remains valid */
518 for (h = 1; h <= R; h++)
519 Q[h] ^= Q[h + K - R];
520 for (h = R + 1; h <= K; h++)
521 Q[h] ^= Q[h - R];
522 } while (t <= K);
523
524 gen->param = param;
525 gen->state = state;
526 gen->GetBits = &GFSR_Bits;
527 gen->GetU01 = &GFSR_U01;
528 gen->Write = &WrGFSR;
529 return gen;
530 }
531
532
533 /**************************************************************************/
534
Fushimi90_Bits(void * junk,void * vsta)535 static unsigned long Fushimi90_Bits (void *junk, void *vsta)
536 {
537 GFSR_state *state = vsta;
538 unsigned long V;
539
540 V = state->X[state->s];
541 state->X[state->s] ^= state->X[state->r];
542 if (++state->s == state->K)
543 state->s = 0;
544 if (++state->r == state->K)
545 state->r = 0;
546 return V << 1;
547 }
548
549 /*-----------------------------------------------------------------------*/
550
Fushimi90_U01(void * vpar,void * vsta)551 static double Fushimi90_U01 (void *vpar, void *vsta)
552 {
553 return Fushimi90_Bits (vpar, vsta) * unif01_INV32;
554 }
555
556 /*-----------------------------------------------------------------------*/
557
ugfsr_CreateFushimi90(int s)558 unif01_Gen *ugfsr_CreateFushimi90 (int s)
559 {
560 unif01_Gen *gen;
561 GFSR_state *state;
562 GFSR_param *param;
563 size_t leng;
564 char name[LEN + 1];
565
566 gen = util_Malloc (sizeof (unif01_Gen));
567 param = util_Malloc (sizeof (GFSR_param));
568 state = util_Malloc (sizeof (GFSR_state));
569
570 state->K = 521;
571 state->r = 521 - 32;
572 state->s = 0;
573
574 strcpy (name, "ugfsr_CreateFushimi90:");
575 addstr_Int (name, " s = ", s);
576 leng = strlen (name);
577 gen->name = util_Calloc (leng + 1, sizeof (char));
578 strncpy (gen->name, name, leng);
579
580 state->X = util_Calloc ((size_t) 3 * 522, sizeof (unsigned long));
581 InitFushimi (521, 32, s, state);
582 gen->param = param;
583 gen->state = state;
584 gen->GetBits = &Fushimi90_Bits;
585 gen->GetU01 = &Fushimi90_U01;
586 gen->Write = &WrGFSR;
587 return gen;
588 }
589
590 /*-----------------------------------------------------------------------*/
591
ugfsr_CreateFushimi(int k,int r,int s)592 unif01_Gen *ugfsr_CreateFushimi (int k, int r, int s)
593 {
594 unif01_Gen *gen;
595 GFSR_state *state;
596 GFSR_param *param;
597 size_t leng;
598 char name[LEN + 1];
599
600 gen = util_Malloc (sizeof (unif01_Gen));
601 param = util_Malloc (sizeof (GFSR_param));
602 state = util_Malloc (sizeof (GFSR_state));
603
604 state->K = k;
605 state->r = k - r;
606 state->s = 0;
607
608 strcpy (name, "ugfsr_CreateFushimi:");
609 addstr_Int (name, " k = ", k);
610 addstr_Int (name, ", r = ", r);
611 addstr_Int (name, ", s = ", s);
612 leng = strlen (name);
613 gen->name = util_Calloc (leng + 1, sizeof (char));
614 strncpy (gen->name, name, leng);
615
616 state->X = util_Calloc ((size_t) k * 3, sizeof (unsigned long));
617 InitFushimi (k, r, s, state);
618 gen->param = param;
619 gen->state = state;
620 gen->GetBits = &Fushimi90_Bits;
621 gen->GetU01 = &Fushimi90_U01;
622 gen->Write = &WrGFSR;
623 return gen;
624 }
625
626
627 /**************************************************************************/
628
ugfsr_CreateKirk81(long s)629 unif01_Gen *ugfsr_CreateKirk81 (long s)
630 {
631 unif01_Gen *gen;
632 GFSR_state *state;
633 GFSR_param *param;
634 size_t leng;
635 char name[LEN + 1];
636 unsigned long mask, msb;
637 unsigned int j, k;
638 long vt;
639
640 gen = util_Malloc (sizeof (unif01_Gen));
641 param = util_Malloc (sizeof (GFSR_param));
642 state = util_Malloc (sizeof (GFSR_state));
643
644 strcpy (name, "ugfsr_CreateKirk81:");
645 addstr_Long (name, " s = ", s);
646 leng = strlen (name);
647 gen->name = util_Calloc (leng + 1, sizeof (char));
648 strncpy (gen->name, name, leng);
649
650 state->K = 250;
651 state->s = 0;
652 state->r = 147;
653 state->X = util_Calloc ((size_t) state->K, sizeof (unsigned long));
654
655 msb = 2147483648UL; /* TwoExp31 */
656 mask = 4294967295UL; /* TwoExp32 - 1 */
657 for (j = 0; j < state->K; j++) {
658 vt = s / 127773;
659 s = 16807 * (s - 127773 * vt) - 2836 * vt;
660 if (s < 0)
661 s += 2147483647;
662 state->X[j] = 2 * ((unsigned long) s);
663 if (s > 1000000000)
664 state->X[j] += 1;
665 /* When s > 1000000000, then set last bit to 1 */
666 }
667 for (j = 1; j <= 31; j++) {
668 k = 7 * j + 3;
669 state->X[k] = (state->X[k] & mask) | msb;
670 mask >>= 1;
671 msb >>= 1;
672 }
673 param->Shift = 0;
674 gen->param = param;
675 gen->state = state;
676 gen->GetBits = &GFSR_Bits;
677 gen->GetU01 = &GFSR_U01;
678 gen->Write = &WrGFSR;
679 return gen;
680 }
681
682
683 /**************************************************************************/
684
TGFSR_Bits(void * vpar,void * vsta)685 static unsigned long TGFSR_Bits (void *vpar, void *vsta)
686 {
687 GFSR_param *param = vpar;
688 GFSR_state *state = vsta;
689 unsigned long v;
690
691 state->X[state->s] = state->X[state->r] ^ (state->X[state->s] >> 1)
692 ^ param->mag01[state->X[state->s] % 2];
693 v = state->X[state->s] & param->Mask;
694 if (++state->s == state->K)
695 state->s = 0; /* s = (s + 1) % K */
696 if (++state->r == state->K)
697 state->r = 0; /* r = (r + 1) % K */
698 return v << param->Shift;
699 }
700
701 /*-----------------------------------------------------------------------*/
702
TGFSR_U01(void * vpar,void * vsta)703 static double TGFSR_U01 (void *vpar, void *vsta)
704 {
705 return TGFSR_Bits (vpar, vsta) * unif01_INV32;
706 }
707
708 /*-----------------------------------------------------------------------*/
709
ugfsr_CreateTGFSR(unsigned int k,unsigned int r,unsigned int l,unsigned long Av,unsigned long S[])710 unif01_Gen *ugfsr_CreateTGFSR (unsigned int k, unsigned int r,
711 unsigned int l, unsigned long Av, unsigned long S[])
712 {
713 unif01_Gen *gen;
714 GFSR_param *param;
715 char name[LEN + 1] = "";
716 size_t leng;
717
718 gen = CreateGFSR0 (k, r, l, S, "ugfsr_CreateTGFSR:");
719 addstr_Ulong (name, ", Av = ", Av);
720 leng = strlen (gen->name) + strlen (name);
721 gen->name = util_Realloc (gen->name, (leng + 1) * sizeof (char));
722 strncat (gen->name, name, leng);
723
724 param = gen->param;
725 param->mag01[0] = 0x0;
726 param->mag01[1] = Av;
727 gen->GetBits = &TGFSR_Bits;
728 gen->GetU01 = &TGFSR_U01;
729 return gen;
730 }
731
732
733 /**************************************************************************/
734
T800_Bits(void * vpar,void * vsta)735 static unsigned long T800_Bits (void *vpar, void *vsta)
736 {
737 GFSR_state *state = vsta;
738 GFSR_param *param = vpar;
739 unsigned long v;
740
741 state->X[state->s] = state->X[state->r] ^ (state->X[state->s] >> 1)
742 ^ param->mag01[state->X[state->s] % 2];
743 v = state->X[state->s];
744 #ifndef IS_ULONG32
745 v &= 0xffffffffUL; /* you may delete this line if word size = 32 */
746 #endif
747 if (++state->s == state->K)
748 state->s = 0; /* s = (s + 1) % K */
749 if (++state->r == state->K)
750 state->r = 0; /* r = (r + 1) % K */
751 return v;
752 }
753
754 /*-----------------------------------------------------------------------*/
755
T800_U01(void * vpar,void * vsta)756 static double T800_U01 (void *vpar, void *vsta)
757 {
758 return T800_Bits (vpar, vsta) * unif01_INV32;
759 }
760
761 /*-----------------------------------------------------------------------*/
762
ugfsr_CreateT800(unsigned long S[])763 unif01_Gen *ugfsr_CreateT800 (unsigned long S[])
764 {
765 unif01_Gen *gen;
766 GFSR_param *param;
767 char name[LEN + 1] = "";
768 size_t leng;
769
770 gen = CreateGFSR0 (25, 18, 32, S, "ugfsr_CreateT800:");
771 addstr_Ulong (name, ", Av = ", 2394935336UL);
772 leng = strlen (gen->name) + strlen (name);
773 gen->name = util_Realloc (gen->name, (leng + 1) * sizeof (char));
774 strncat (gen->name, name, leng);
775
776 param = gen->param;
777 param->mag01[0] = 0x0;
778 param->mag01[1] = 0x8ebfd028UL; /* this is magic vector 'a' */
779 gen->GetBits = &T800_Bits;
780 gen->GetU01 = &T800_U01;
781 return gen;
782 }
783
784
785 /***************************************************************************/
786
TGFSR2_Bits(void * vpar,void * vsta)787 static unsigned long TGFSR2_Bits (void *vpar, void *vsta)
788 {
789 TGFSR2_param *param = vpar;
790 TGFSR2_state *state = vsta;
791 unsigned long v;
792
793 v = state->X[state->s];
794 state->X[state->s] =
795 state->X[state->r] ^ (v >> 1) ^ param->mag01[v % 2];
796 if (++state->s == state->K)
797 state->s = 0; /* s = (s + 1) % K */
798 if (++state->r == state->K)
799 state->r = 0; /* r = (r + 1) % K */
800 v ^= (v << param->S2) & param->b; /* s and b, magic vectors */
801 v ^= (v << param->T2) & param->c; /* t and c, magic vectors */
802 v &= param->Mask; /* l least significant bits */
803 return v << param->Shift;
804 }
805
806 /*-----------------------------------------------------------------------*/
807
TGFSR2_U01(void * vpar,void * vsta)808 static double TGFSR2_U01 (void *vpar, void *vsta)
809 {
810 return TGFSR2_Bits (vpar, vsta) * unif01_INV32;
811 }
812
813 /*-----------------------------------------------------------------------*/
814
ugfsr_CreateTGFSR2(unsigned int k,unsigned int r,unsigned int l,unsigned int s,unsigned int t,unsigned long Av,unsigned long Bv,unsigned long Cv,unsigned long S[])815 unif01_Gen *ugfsr_CreateTGFSR2 (unsigned int k, unsigned int r,
816 unsigned int l, unsigned int s, unsigned int t, unsigned long Av,
817 unsigned long Bv, unsigned long Cv, unsigned long S[])
818 {
819 unif01_Gen *gen;
820 TGFSR2_param *param;
821 char name[LEN + 1];
822 size_t leng;
823
824 gen = CreateGFSR0 (k, r, l, S, "");
825 util_Free (gen->name);
826 strcpy (name, "ugfsr_CreateTGFSR2:");
827 addstr_Uint (name, " k = ", k);
828 addstr_Uint (name, ", r = ", r);
829 addstr_Uint (name, ", l = ", l);
830 addstr_Ulong (name, ", Av = ", Av);
831 addstr_Ulong (name, ", Bv = ", Bv);
832 addstr_Ulong (name, ", Cv = ", Cv);
833 addstr_Uint (name, ", s = ", s);
834 addstr_Uint (name, ", t = ", t);
835 addstr_ArrayUlong (name, ", S", (int) k, S);
836 leng = strlen (name);
837 gen->name = util_Calloc (leng + 1, sizeof (char));
838 strncpy (gen->name, name, leng);
839
840 util_Free (gen->param);
841 gen->param = param = util_Malloc (sizeof (TGFSR2_param));
842 param->S2 = s;
843 param->T2 = t;
844 param->mag01[0] = 0x0;
845 param->mag01[1] = Av; /* this is magic vector 'a' */
846 param->b = Bv;
847 param->c = Cv;
848 param->Shift = 32 - l;
849 param->Mask = num_TwoExp[l] - 1.0;
850 if (l == 32)
851 param->Mask = unif01_MASK32;
852 gen->GetBits = &TGFSR2_Bits;
853 gen->GetU01 = &TGFSR2_U01;
854 gen->Write = &WrGFSR;
855 return gen;
856 }
857
858
859 /***************************************************************************/
860
TT400_U01(void * vpar,void * vsta)861 static double TT400_U01 (void *vpar, void *vsta)
862 {
863 GFSR_param *param = vpar;
864 GFSR_state *state = vsta;
865 unsigned long v;
866
867 v = state->X[state->s];
868 state->X[state->s] = state->X[state->r] ^ (v >> 1) ^ param->mag01[v % 2];
869 v ^= (v << 2) & 0x6A68; /* s and b, magic vectors */
870 v ^= (v << 7) & 0x7500; /* t and c, magic vectors */
871 v &= 0xffff;
872 if (++state->s == state->K)
873 state->s = 0;
874 if (++state->r == state->K)
875 state->r = 0;
876 return (double) v / 0xffff;
877 }
878
879 /*-----------------------------------------------------------------------*/
880
TT400_Bits(void * vpar,void * vsta)881 static unsigned long TT400_Bits (void *vpar, void *vsta)
882 {
883 return (unsigned long) (TT400_U01 (vpar, vsta) * unif01_NORM32);
884 }
885
886 /*-----------------------------------------------------------------------*/
887
ugfsr_CreateTT400(unsigned long S[])888 unif01_Gen *ugfsr_CreateTT400 (unsigned long S[])
889 {
890 unif01_Gen *gen;
891 GFSR_param *param;
892 unsigned int N, M;
893
894 N = 25;
895 M = 25 - 11;
896 gen = CreateGFSR0 (N, M, 16, S, "ugfsr_CreateTT400:");
897
898 param = gen->param;
899 param->mag01[0] = 0x0;
900 param->mag01[1] = 0xA875; /* this is magic vector 'a' */
901 gen->GetBits = &TT400_Bits;
902 gen->GetU01 = &TT400_U01;
903 return gen;
904 }
905
906
907 /***************************************************************************/
908
TT403_U01(void * vpar,void * vsta)909 static double TT403_U01 (void *vpar, void *vsta)
910 {
911 GFSR_param *param = vpar;
912 GFSR_state *state = vsta;
913 unsigned long v;
914
915 v = state->X[state->s];
916 state->X[state->s] = state->X[state->r] ^ (v >> 1) ^ param->mag01[v % 2];
917 v ^= (v << 8) & 0x102D1200; /* s and b, magic vectors */
918 v ^= (v << 14) & 0x66E50000; /* t and c, magic vectors */
919 v &= 0x7fffffff;
920 if (++state->s == state->K)
921 state->s = 0;
922 if (++state->r == state->K)
923 state->r = 0;
924 return (double) v / 0x7fffffff;
925 }
926
927 /*-----------------------------------------------------------------------*/
928
TT403_Bits(void * vpar,void * vsta)929 static unsigned long TT403_Bits (void *vpar, void *vsta)
930 {
931 return (unsigned long) (TT403_U01 (vpar, vsta) * unif01_NORM32);
932 }
933
934 /*-----------------------------------------------------------------------*/
935
ugfsr_CreateTT403(unsigned long S[])936 unif01_Gen *ugfsr_CreateTT403 (unsigned long S[])
937 {
938 unif01_Gen *gen;
939 GFSR_param *param;
940 unsigned int N, M;
941
942 N = 13;
943 M = 13 - 2;
944 gen = CreateGFSR0 (N, M, 31, S, "ugfsr_CreateTT403:");
945
946 param = gen->param;
947 param->mag01[0] = 0x0;
948 param->mag01[1] = 0x6B5ECCF6; /* this is magic vector 'a' */
949 gen->GetBits = &TT403_Bits;
950 gen->GetU01 = &TT403_U01;
951 return gen;
952 }
953
954
955 /***************************************************************************/
956
TT775_U01(void * vpar,void * vsta)957 static double TT775_U01 (void *vpar, void *vsta)
958 {
959 GFSR_param *param = vpar;
960 GFSR_state *state = vsta;
961 unsigned long v;
962
963 v = state->X[state->s];
964 state->X[state->s] = state->X[state->r] ^ (v >> 1) ^ param->mag01[v % 2];
965 v ^= (v << 6) & 0x1ABD5900; /* s and b, magic vectors */
966 v ^= (v << 14) & 0x776A0000; /* t and c, magic vectors */
967 v &= 0x7fffffff;
968 if (++state->s == state->K)
969 state->s = 0;
970 if (++state->r == state->K)
971 state->r = 0;
972 return (double) v / 0x7fffffff;
973 }
974
975 /*-----------------------------------------------------------------------*/
976
TT775_Bits(void * vpar,void * vsta)977 static unsigned long TT775_Bits (void *vpar, void *vsta)
978 {
979 return (unsigned long) (TT775_U01 (vpar, vsta) * unif01_NORM32);
980 }
981
982 /*-----------------------------------------------------------------------*/
983
ugfsr_CreateTT775(unsigned long S[])984 unif01_Gen *ugfsr_CreateTT775 (unsigned long S[])
985 {
986 unif01_Gen *gen;
987 GFSR_param *param;
988 unsigned int N, M;
989
990 N = 25;
991 M = 25 - 8;
992 gen = CreateGFSR0 (N, M, 31, S, "ugfsr_CreateTT775:");
993
994 param = gen->param;
995 param->mag01[0] = 0x0;
996 param->mag01[1] = 0x6C6CB38C; /* this is magic vector 'a' */
997 gen->GetBits = &TT775_Bits;
998 gen->GetU01 = &TT775_U01;
999 return gen;
1000 }
1001
1002
1003 /**************************************************************************
1004 *
1005 * Our version of TT800
1006 *
1007 **************************************************************************/
1008
TT800_Bits(void * vpar,void * vsta)1009 static unsigned long TT800_Bits (void *vpar, void *vsta)
1010 {
1011 GFSR_param *param = vpar;
1012 GFSR_state *state = vsta;
1013 unsigned long v;
1014
1015 v = state->X[state->s];
1016 state->X[state->s] = state->X[state->r] ^ (v >> 1) ^ param->mag01[v % 2];
1017 v ^= (v << 7) & 0x2b5b2500; /* s and b, magic vectors */
1018 v ^= (v << 15) & 0xdb8b0000UL; /* t and c, magic vectors */
1019 #ifndef IS_ULONG32
1020 v &= 0xffffffffUL; /* you may delete this line if word size = 32 */
1021 #endif
1022 if (++state->s == state->K)
1023 state->s = 0;
1024 if (++state->r == state->K)
1025 state->r = 0;
1026 return v;
1027 }
1028
1029 /*-----------------------------------------------------------------------*/
1030
TT800_U01(void * vpar,void * vsta)1031 static double TT800_U01 (void *vpar, void *vsta)
1032 {
1033 return TT800_Bits (vpar, vsta) * unif01_INV32;
1034 }
1035
1036 /*-----------------------------------------------------------------------*/
1037
ugfsr_CreateTT800(unsigned long S[])1038 unif01_Gen *ugfsr_CreateTT800 (unsigned long S[])
1039 {
1040 unif01_Gen *gen;
1041 GFSR_param *param;
1042 unsigned int N, M;
1043
1044 N = 25;
1045 M = 25 - 7;
1046 gen = CreateGFSR0 (N, M, 32, S, "ugfsr_CreateTT800:");
1047
1048 param = gen->param;
1049 param->mag01[0] = 0x0;
1050 param->mag01[1] = 0x8ebfd028UL; /* this is magic vector 'a' */
1051 gen->GetBits = &TT800_Bits;
1052 gen->GetU01 = &TT800_U01;
1053 return gen;
1054 }
1055
1056
1057 /***************************************************************************/
1058
1059 /* A C-program for TT800 : July 8th 1994 Version */
1060 /* by M. Matsumoto, email: matumoto@math.keio.ac.jp */
1061 /* genrand() generate one pseudorandom number with double precision */
1062 /* which is uniformly distributed on [0,1]-interval */
1063 /* for each call. One may choose any initial 25 seeds */
1064 /* except all zeros. */
1065
1066 #define NN 25
1067 #define MM 7
1068
TT800M94_U01(void * vpar,void * vsta)1069 static double TT800M94_U01 (void *vpar, void *vsta)
1070 {
1071 GFSR_param *param = vpar;
1072 GFSR_state *state = vsta;
1073 unsigned long y;
1074 unsigned int j;
1075
1076 if (state->s == NN) {
1077 /* generate NN words at one time */
1078 for (j = 0; j < NN - MM; j++) {
1079 state->X[j] = state->X[j + MM] ^ (state->X[j] >> 1) ^
1080 param->mag01[state->X[j] % 2];
1081 }
1082 for (; j < NN; j++) {
1083 state->X[j] = state->X[j + (MM - NN)] ^ (state->X[j] >> 1) ^
1084 param->mag01[state->X[j] % 2];
1085 }
1086 state->s = 0;
1087 }
1088 y = state->X[state->s++];
1089 y ^= (y << 7) & 0x2b5b2500; /* s and b, magic vectors */
1090 y ^= (y << 15) & 0xdb8b0000UL; /* t and c, magic vectors */
1091 #ifndef IS_ULONG32
1092 y &= 0xffffffffUL; /* you may delete this line if word size = 32 */
1093 #endif
1094
1095 return (double) y / 0xffffffffUL;
1096 }
1097
1098 /*-----------------------------------------------------------------------*/
1099
TT800M94_Bits(void * vpar,void * vsta)1100 static unsigned long TT800M94_Bits (void *vpar, void *vsta)
1101 {
1102 return (unsigned long) (TT800M94_U01 (vpar, vsta) * unif01_NORM32);
1103 }
1104
1105 /*-----------------------------------------------------------------------*/
1106
1107
1108 /* A C-program for TT800 : July 8th 1996 Version */
1109 /* by M. Matsumoto, email: matumoto@math.keio.ac.jp */
1110 /* genrand() generate one pseudorandom number with double precision */
1111 /* which is uniformly distributed on [0,1]-interval */
1112 /* for each call. One may choose any initial 25 seeds */
1113 /* except all zeros. */
1114
TT800M96_U01(void * vpar,void * vsta)1115 static double TT800M96_U01 (void *vpar, void *vsta)
1116 {
1117 GFSR_param *param = vpar;
1118 GFSR_state *state = vsta;
1119 unsigned long y;
1120 unsigned int j;
1121
1122 if (state->s == NN) {
1123 /* generate NN words at one time */
1124 for (j = 0; j < NN - MM; j++) {
1125 state->X[j] = state->X[j + MM] ^ (state->X[j] >> 1) ^
1126 param->mag01[state->X[j] % 2];
1127 }
1128 for (; j < NN; j++) {
1129 state->X[j] = state->X[j + (MM - NN)] ^ (state->X[j] >> 1) ^
1130 param->mag01[state->X[j] % 2];
1131 }
1132 state->s = 0;
1133 }
1134 y = state->X[state->s++];
1135 y ^= (y << 7) & 0x2b5b2500; /* s and b, magic vectors */
1136 y ^= (y << 15) & 0xdb8b0000UL; /* t and c, magic vectors */
1137 #ifndef IS_ULONG32
1138 y &= 0xffffffffUL; /* you may delete this line if word
1139 size = 32 */
1140 #endif
1141
1142 /* the following line was added by Makoto Matsumoto in the 1996 version
1143 to improve lower bit's correlation. This line differs from the code
1144 published in 1994. */
1145 y ^= (y >> 16); /* added to the 1994 version */
1146
1147 return (double) y / 0xffffffffUL;
1148 }
1149
1150 /*-----------------------------------------------------------------------*/
1151
TT800M96_Bits(void * vpar,void * vsta)1152 static unsigned long TT800M96_Bits (void *vpar, void *vsta)
1153 {
1154 return (unsigned long) (TT800M96_U01 (vpar, vsta) * unif01_NORM32);
1155 }
1156
1157 /*-----------------------------------------------------------------------*/
1158
ugfsr_CreateTT800M94(unsigned long S[])1159 unif01_Gen *ugfsr_CreateTT800M94 (unsigned long S[])
1160 {
1161 unif01_Gen *gen;
1162 GFSR_param *param;
1163
1164 gen = CreateGFSR0 (NN, MM, 32, S, "ugfsr_CreateTT800M94:");
1165 param = gen->param;
1166 param->mag01[0] = 0x0;
1167 param->mag01[1] = 0x8ebfd028UL; /* this is magic vector 'a' */
1168 gen->GetBits = &TT800M94_Bits;
1169 gen->GetU01 = &TT800M94_U01;
1170 return gen;
1171 }
1172
1173 /*-----------------------------------------------------------------------*/
1174
ugfsr_CreateTT800M96(unsigned long S[])1175 unif01_Gen *ugfsr_CreateTT800M96 (unsigned long S[])
1176 {
1177 unif01_Gen *gen;
1178 GFSR_param *param;
1179
1180 gen = CreateGFSR0 (NN, MM, 32, S, "ugfsr_CreateTT800M96:");
1181 param = gen->param;
1182 param->mag01[0] = 0x0;
1183 param->mag01[1] = 0x8ebfd028UL; /* this is magic vector 'a' */
1184 gen->GetBits = &TT800M96_Bits;
1185 gen->GetU01 = &TT800M96_U01;
1186 return gen;
1187 }
1188
1189
1190 /**************************************************************************
1191 *
1192 * The code for the two versions of MT19937 is copyrighted by Makoto
1193 * Matsumoto and Takuji Nishimura. I have changed the name of some variables
1194 * to make it consistent with TestU01. (R. Simard)
1195 */
1196
1197 /*
1198 A C-program for MT19937, with initialization improved 2002/1/26.
1199 Coded by Takuji Nishimura and Makoto Matsumoto.
1200
1201 Before using, initialize the state by using init_genrand(seed)
1202 or init_by_array(init_key, key_length).
1203
1204 Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
1205 All rights reserved.
1206
1207 Redistribution and use in source and binary forms, with or without
1208 modification, are permitted provided that the following conditions
1209 are met:
1210
1211 1. Redistributions of source code must retain the above copyright
1212 notice, this list of conditions and the following disclaimer.
1213
1214 2. Redistributions in binary form must reproduce the above copyright
1215 notice, this list of conditions and the following disclaimer in the
1216 documentation and/or other materials provided with the distribution.
1217
1218 3. The names of its contributors may not be used to endorse or promote
1219 products derived from this software without specific prior written
1220 permission.
1221
1222 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
1223 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
1224 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
1225 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
1226 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
1227 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
1228 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
1229 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
1230 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
1231 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
1232 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
1233
1234
1235 Any feedback is very welcome.
1236 http://www.math.keio.ac.jp/matumoto/emt.html
1237 email: matumoto@math.keio.ac.jp
1238 */
1239
1240 /* Period parameters */
1241 #undef NN
1242 #undef MM
1243 #define NN 624
1244 #define MM 397
1245 #define MATRIX_A 0x9908b0dfUL /* constant vector a */
1246 #define UPPER_MASK 0x80000000UL /* most significant w-r bits */
1247 #undef LOWER_MASK
1248 #define LOWER_MASK 0x7fffffffUL /* least significant r bits */
1249
1250
1251 /*-----------------------------------------------------------------------*/
1252
1253 /* initializes mt[NN] with a seed */
init_genrand(void * vsta,unsigned long s)1254 static void init_genrand (void *vsta, unsigned long s)
1255 {
1256 MT19937_state *state = vsta;
1257 int j;
1258
1259 state->X[0] = s & 0xffffffffUL;
1260 for (j = 1; j < NN; j++) {
1261 state->X[j] =
1262 (1812433253UL * (state->X[j - 1] ^ (state->X[j - 1] >> 30)) + j);
1263 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
1264 /* In the previous versions, MSBs of the seed affect */
1265 /* only MSBs of the array state->X[]. */
1266 /* 2002/01/09 modified by Makoto Matsumoto */
1267 state->X[j] &= 0xffffffffUL;
1268 /* for >32 bit machines */
1269 }
1270 state->s = NN;
1271 }
1272
1273
1274 /*-----------------------------------------------------------------------*/
1275
1276 /* initialize by an array with array-length */
1277 /* init_key is the array for initializing keys */
1278 /* key_length is its length */
init_by_array(void * vsta,unsigned long init_key[],int key_length)1279 static void init_by_array (void *vsta, unsigned long init_key[],
1280 int key_length)
1281 {
1282 MT19937_state *state = vsta;
1283 int i, j, k;
1284 init_genrand (vsta, 19650218UL);
1285 i = 1;
1286 j = 0;
1287 k = (NN > key_length ? NN : key_length);
1288 for (; k; k--) {
1289 state->X[i] = (state->X[i] ^ ((state->X[i - 1] ^ (state->X[i - 1] >>
1290 30)) * 1664525UL)) + init_key[j] + j; /* non linear */
1291 state->X[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
1292 i++;
1293 j++;
1294 if (i >= NN) {
1295 state->X[0] = state->X[NN - 1];
1296 i = 1;
1297 }
1298 if (j >= key_length)
1299 j = 0;
1300 }
1301 for (k = NN - 1; k; k--) {
1302 state->X[i] = (state->X[i] ^ ((state->X[i - 1] ^ (state->X[i - 1] >>
1303 30)) * 1566083941UL)) - i; /* non linear */
1304 state->X[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
1305 i++;
1306 if (i >= NN) {
1307 state->X[0] = state->X[NN - 1];
1308 i = 1;
1309 }
1310 }
1311
1312 state->X[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial
1313 array */
1314 }
1315
1316
1317 /*-----------------------------------------------------------------------*/
1318
1319 /* generates a random number on [0,0xffffffff]-interval
1320 unsigned long genrand_int32(void) */
MT19937_02_Bits(void * vpar,void * vsta)1321 static unsigned long MT19937_02_Bits (void *vpar, void *vsta)
1322 {
1323 MT19937_param *param = vpar;
1324 MT19937_state *state = vsta;
1325 unsigned long y;
1326
1327 if (state->s >= NN) { /* generate NN words at one time */
1328 int kk;
1329
1330 if (state->s == NN + 1) /* if init_genrand() has not been called, */
1331 init_genrand (state, 5489UL); /* a default initial seed is used */
1332
1333 for (kk = 0; kk < NN - MM; kk++) {
1334 y = (state->X[kk] & UPPER_MASK) | (state->X[kk + 1] & LOWER_MASK);
1335 state->X[kk] =
1336 state->X[kk + MM] ^ (y >> 1) ^ param->mag01[y & 0x1UL];
1337 }
1338 for (; kk < NN - 1; kk++) {
1339 y = (state->X[kk] & UPPER_MASK) | (state->X[kk + 1] & LOWER_MASK);
1340 state->X[kk] = state->X[kk + (MM - NN)] ^ (y >> 1) ^
1341 param->mag01[y & 0x1UL];
1342 }
1343 y = (state->X[NN - 1] & UPPER_MASK) | (state->X[0] & LOWER_MASK);
1344 state->X[NN - 1] =
1345 state->X[MM - 1] ^ (y >> 1) ^ param->mag01[y & 0x1UL];
1346
1347 state->s = 0;
1348 }
1349
1350 y = state->X[state->s++];
1351
1352 /* Tempering */
1353 y ^= (y >> 11);
1354 y ^= (y << 7) & 0x9d2c5680UL;
1355 y ^= (y << 15) & 0xefc60000UL;
1356 y ^= (y >> 18);
1357
1358 return y;
1359 }
1360
1361
1362 /*-----------------------------------------------------------------------*/
1363
1364 /* generates a random number on (0,1)-real-interval
1365 double genrand_real3(void) */
MT19937_02_U01(void * vpar,void * vsta)1366 static double MT19937_02_U01 (void *vpar, void *vsta)
1367 {
1368 return ((double) MT19937_02_Bits (vpar, vsta) + 0.5) *
1369 (1.0 / 4294967296.0); /* divided by 2^32 */
1370 }
1371
1372 /* These real versions are due to Isaku Wada, 2002/01/09 added */
1373
1374
1375 /*-----------------------------------------------------------------------*/
1376
1377 static void WrMT19937 (void *);
1378
1379 /*-----------------------------------------------------------------------*/
1380
ugfsr_CreateMT19937_02(unsigned long seed,unsigned long Key[],int len)1381 unif01_Gen *ugfsr_CreateMT19937_02 (unsigned long seed,
1382 unsigned long Key[], int len)
1383 {
1384 unif01_Gen *gen;
1385 MT19937_param *param;
1386 MT19937_state *state;
1387 unsigned long S[NN];
1388 char name[LEN + 1];
1389 size_t leng;
1390
1391 gen = CreateGFSR0 (NN, MM, 32, S, "");
1392 param = gen->param;
1393 state = gen->state;
1394 param->mag01[0] = 0x0UL;
1395 param->mag01[1] = MATRIX_A;
1396 gen->GetBits = &MT19937_02_Bits;
1397 gen->GetU01 = &MT19937_02_U01;
1398 gen->Write = &WrMT19937;
1399 strcpy (name, "ugfsr_CreateMT19937_02:");
1400 if ((len <= 0) || (NULL == Key)) {
1401 init_genrand (state, seed);
1402 addstr_Ulong (name, " seed = ", seed);
1403 } else {
1404 state->s = NN + 1;
1405 init_by_array (state, Key, len);
1406 addstr_ArrayUlong (name, " Key = ", len, Key);
1407 }
1408 leng = strlen (name);
1409 gen->name = util_Realloc (gen->name, (leng + 1) * sizeof (char));
1410 strncpy (gen->name, name, leng);
1411 gen->name[leng] = '\0';
1412 return gen;
1413 }
1414
1415
1416
1417 /**************************************************************************/
1418 /* A C-program for MT19937: Real number version - */
1419 /* generates one pseudorandom number with double precision */
1420 /* which is uniformly distributed on [0,1]-interval for each call. */
1421 /* CreateMT19937(seed) set initial values to the working area of 624 words. */
1422 /* (seed is any integer except for 0). */
1423
1424 #undef NN
1425 #undef MM
1426 #define NN 624
1427 #define MM 397
1428
MT19937_98_U01(void * vpar,void * vsta)1429 static double MT19937_98_U01 (void *vpar, void *vsta)
1430 {
1431 MT19937_param *param = vpar;
1432 MT19937_state *state = vsta;
1433 unsigned long y;
1434 unsigned int j;
1435
1436 if (state->s == NN) { /* generate NN words at one time */
1437 for (j = 0; j < NN - MM; j++) {
1438 y = (state->X[j] & UPPER_MASK) | (state->X[j + 1] & LOWER_MASK);
1439 state->X[j] = state->X[j + MM] ^ (y >> 1) ^ param->mag01[y & 0x1];
1440 }
1441 for (; j < NN - 1; j++) {
1442 y = (state->X[j] & UPPER_MASK) | (state->X[j + 1] & LOWER_MASK);
1443 state->X[j] = state->X[j + (MM - NN)] ^ (y >> 1) ^
1444 param->mag01[y & 0x1];
1445 }
1446 y = (state->X[NN - 1] & UPPER_MASK) | (state->X[0] & LOWER_MASK);
1447 state->X[NN - 1] = state->X[MM - 1] ^ (y >> 1) ^ param->mag01[y & 0x1];
1448
1449 state->s = 0;
1450 }
1451 y = state->X[state->s++];
1452 y ^= TEMPERING_SHIFT_U (y);
1453 y ^= TEMPERING_SHIFT_S (y) & TEMPERING_MASK_B;
1454 y ^= TEMPERING_SHIFT_T (y) & TEMPERING_MASK_C;
1455 #ifndef IS_ULONG32
1456 y &= 0xffffffffUL; /* you may delete this line if word size = 32 */
1457 #endif
1458 y ^= TEMPERING_SHIFT_L (y);
1459 return (double) y / 0xffffffffUL;
1460 }
1461
1462 /*-----------------------------------------------------------------------*/
1463
MT19937_98_Bits(void * vpar,void * vsta)1464 static unsigned long MT19937_98_Bits (void *vpar, void *vsta)
1465 {
1466 return (unsigned long) (MT19937_98_U01 (vpar, vsta) * unif01_NORM32);
1467 }
1468
1469 /*-----------------------------------------------------------------------*/
1470
WrMT19937(void * vsta)1471 static void WrMT19937 (void *vsta)
1472 {
1473 GFSR_state *state = vsta;
1474 unsigned int j;
1475 if (unif01_WrLongStateFlag) {
1476 printf ("S = {\n ");
1477 for (j = 0; j < state->K; j++) {
1478 printf (" %12lu", state->X[j]);
1479 if (j < state->K - 1)
1480 printf (",");
1481 if ((j % 5) == 4)
1482 printf ("\n ");
1483 };
1484 printf (" }\n");
1485 } else
1486 unif01_WrLongStateDef ();
1487 }
1488
1489 /*-----------------------------------------------------------------------*/
1490
ugfsr_CreateMT19937_98(unsigned long seed)1491 unif01_Gen *ugfsr_CreateMT19937_98 (unsigned long seed)
1492 /*
1493 * Setting initial seeds to X[NN] using
1494 * the generator Line 25 of Table 1 in
1495 * [KNUTH 1981, The Art of Computer Programming
1496 * Vol. 2 (2nd Ed.), pp102]
1497 */
1498 {
1499 unif01_Gen *gen;
1500 MT19937_param *param;
1501 unsigned int j;
1502 unsigned long S[NN];
1503 char name[LEN + 1];
1504 size_t leng;
1505
1506 S[0] = seed & 0xffffffffUL;
1507 for (j = 1; j < NN; j++)
1508 S[j] = (69069 * S[j - 1]) & 0xffffffffUL;
1509 gen = CreateGFSR0 (NN, MM, 32, S, "");
1510 param = gen->param;
1511 param->mag01[0] = 0x0;
1512 param->mag01[1] = MATRIX_A;
1513 gen->GetBits = &MT19937_98_Bits;
1514 gen->GetU01 = &MT19937_98_U01;
1515 gen->Write = &WrMT19937;
1516 strcpy (name, "ugfsr_CreateMT19937_98:");
1517 addstr_Ulong (name, " seed = ", seed);
1518 leng = strlen (name);
1519 gen->name = util_Realloc (gen->name, (leng + 1) * sizeof (char));
1520 strncpy (gen->name, name, leng);
1521 gen->name[leng] = '\0';
1522 return gen;
1523 }
1524
1525
1526 /**************************************************************************/
1527
GFSR5_Bits(void * vpar,void * vsta)1528 static unsigned long GFSR5_Bits (void *vpar, void *vsta)
1529 {
1530 GFSR5_param *param = vpar;
1531 GFSR5_state *state = vsta;
1532
1533 if (++state->s == state->K)
1534 state->s = 0;
1535 if (++state->r1 == state->K)
1536 state->r1 = 0;
1537 if (++state->r2 == state->K)
1538 state->r2 = 0;
1539 if (++state->r3 == state->K)
1540 state->r3 = 0;
1541 state->X[state->s] ^= state->X[state->r1] ^ state->X[state->r2] ^
1542 state->X[state->r3];
1543 return state->X[state->s] << param->Shift;
1544 }
1545
1546 /*-----------------------------------------------------------------------*/
1547
GFSR5_U01(void * vpar,void * vsta)1548 static double GFSR5_U01 (void *vpar, void *vsta)
1549 {
1550 return GFSR5_Bits (vpar, vsta) * unif01_INV32;
1551 }
1552
1553 /*-----------------------------------------------------------------------*/
1554
ugfsr_CreateGFSR5(unsigned int k,unsigned int r1,unsigned int r2,unsigned int r3,unsigned int l,unsigned long S[])1555 unif01_Gen *ugfsr_CreateGFSR5 (unsigned int k, unsigned int r1,
1556 unsigned int r2, unsigned int r3, unsigned int l, unsigned long S[])
1557 {
1558 unif01_Gen *gen;
1559 GFSR5_param *param;
1560 GFSR5_state *state;
1561 size_t leng;
1562 char name[LEN + 1];
1563 unsigned int j;
1564 unsigned long Mask;
1565
1566 util_Assert ((l >= 1) && (l <= 32),
1567 "ugfsr_CreateGFSR5: l must be in [1..32]");
1568 util_Assert ((r3 > 0) && (r3 < r2),
1569 "ugfsr_CreateGFSR5: we must have 0 < r3 < r2");
1570 util_Assert ((r1 > r2) && (r1 < k),
1571 "ugfsr_CreateGFSR5: we must have r2 < r1 < k");
1572
1573 gen = util_Malloc (sizeof (unif01_Gen));
1574 param = util_Malloc (sizeof (GFSR5_param));
1575 state = util_Malloc (sizeof (GFSR5_state));
1576
1577 strcpy (name, "ugfsr_CreateGFSR5:");
1578 addstr_Uint (name, " k = ", k);
1579 addstr_Uint (name, ", r1 = ", r1);
1580 addstr_Uint (name, ", r2 = ", r2);
1581 addstr_Uint (name, ", r3 = ", r3);
1582 addstr_Uint (name, ", l = ", l);
1583 addstr_ArrayUlong (name, ", S = ", (int) k, S);
1584 leng = strlen (name);
1585 gen->name = util_Calloc (leng + 1, sizeof (char));
1586 strncpy (gen->name, name, leng);
1587
1588 if (l == 32)
1589 Mask = unif01_MASK32;
1590 else
1591 Mask = num_TwoExp[l] - 1.0;
1592 state->X = util_Calloc ((size_t) k, sizeof (unsigned long));
1593 for (j = 0; j < k; j++)
1594 state->X[j] = S[j] & Mask;
1595
1596 state->r1 = k - r1;
1597 state->r2 = k - r2;
1598 state->r3 = k - r3;
1599 state->s = 0;
1600 state->K = k;
1601 param->Shift = 32 - l;
1602
1603 gen->param = param;
1604 gen->state = state;
1605 gen->GetBits = &GFSR5_Bits;
1606 gen->GetU01 = &GFSR5_U01;
1607 gen->Write = &WrGFSR;
1608 return gen;
1609 }
1610
1611
1612 /**************************************************************************/
1613
1614 #define A 471
1615 #define B 1586
1616 #define C 6988
1617 #define D 9689
1618 #define M 16383
1619
1620 #define RandomInteger (++state->s, state->X[state->s & M] = \
1621 state->X[(state->s-A) & M] ^ state->X[(state->s-B) & M] ^ \
1622 state->X[(state->s-C) & M] ^ state->X[(state->s-D) & M])
1623
WrZiff98(void * vsta)1624 static void WrZiff98 (void *vsta)
1625 {
1626 GFSR_state *state = vsta;
1627 unsigned int j;
1628 int s = state->s;
1629
1630 s = (s - D) % (M + 1);
1631 if (unif01_WrLongStateFlag) {
1632 printf (" S = {\n ");
1633 for (j = 0; j < state->K; j++) {
1634 if (++s > M)
1635 s = 0;
1636 printf (" %12lu", state->X[s]);
1637 if (j < state->K - 1)
1638 printf (",");
1639 if ((j % 5) == 4)
1640 printf ("\n ");
1641 };
1642 printf (" }\n");
1643 } else
1644 unif01_WrLongStateDef ();
1645 }
1646
1647 /*-----------------------------------------------------------------------*/
1648
Ziff98_Bits(void * junk,void * vsta)1649 static unsigned long Ziff98_Bits (void *junk, void *vsta)
1650 {
1651 GFSR_state *state = vsta;
1652 return RandomInteger;
1653 }
1654
1655 /*-----------------------------------------------------------------------*/
1656
Ziff98_U01(void * junk,void * vsta)1657 static double Ziff98_U01 (void *junk, void *vsta)
1658 {
1659 GFSR_state *state = vsta;
1660 return RandomInteger * unif01_INV32;
1661 }
1662
1663 /*-----------------------------------------------------------------------*/
1664
ugfsr_CreateZiff98(unsigned long S[])1665 unif01_Gen *ugfsr_CreateZiff98 (unsigned long S[])
1666 {
1667 unif01_Gen *gen;
1668 GFSR_state *state;
1669 size_t leng;
1670 char name[LEN + 1];
1671 unsigned int j, l = 32, Mask;
1672
1673 gen = util_Malloc (sizeof (unif01_Gen));
1674 state = util_Malloc (sizeof (GFSR_state));
1675
1676 Mask = unif01_MASK32;
1677 state->X = util_Calloc ((size_t) M + 1, sizeof (unsigned long));
1678 for (j = 0; j < D; j++)
1679 state->X[j] = S[j] & Mask;
1680 state->s = D;
1681 state->K = D;
1682
1683 strcpy (name, "ugfsr_CreateZiff98:");
1684 addstr_Uint (name, " k = ", D);
1685 addstr_Uint (name, ", r1 = ", C);
1686 addstr_Uint (name, ", r2 = ", B);
1687 addstr_Uint (name, ", r3 = ", A);
1688 addstr_Uint (name, ", l = ", l);
1689 addstr_ArrayUlong (name, ", S = ", D, S);
1690 leng = strlen (name);
1691 gen->name = util_Calloc (leng + 1, sizeof (char));
1692 strncpy (gen->name, name, leng);
1693
1694 gen->param = NULL;
1695 gen->state = state;
1696 gen->GetBits = &Ziff98_Bits;
1697 gen->GetU01 = &Ziff98_U01;
1698 gen->Write = &WrZiff98;
1699 return gen;
1700 }
1701
1702 /**************************************************************************/
1703
ugfsr_DeleteGFSR5(unif01_Gen * gen)1704 void ugfsr_DeleteGFSR5 (unif01_Gen *gen)
1705 {
1706 GFSR5_state *state;
1707 if (NULL == gen)
1708 return;
1709 state = gen->state;
1710 util_Free (state->X);
1711 gen->state = util_Free (gen->state);
1712 gen->param = util_Free (gen->param);
1713 gen->name = util_Free (gen->name);
1714 util_Free (gen);
1715 }
1716
1717
1718 /**************************************************************************/
1719
ugfsr_DeleteGen(unif01_Gen * gen)1720 void ugfsr_DeleteGen (unif01_Gen *gen)
1721 {
1722 GFSR_state *state;
1723 if (NULL == gen)
1724 return;
1725 state = gen->state;
1726 util_Free (state->X);
1727 gen->state = util_Free (gen->state);
1728 gen->param = util_Free (gen->param);
1729 gen->name = util_Free (gen->name);
1730 util_Free (gen);
1731 }
1732