1 /*************************************************************************\
2 *
3 * Package: TestU01
4 * File: umrg.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 "gdef.h"
32 #include "util.h"
33 #include "num.h"
34 #include "addstr.h"
35
36 #include "umrg.h"
37 #include "unif01.h"
38
39 #include <stdio.h>
40 #include <stdlib.h>
41 #include <string.h>
42 #include <float.h>
43
44
45 #ifdef USE_GMP
46 #include <gmp.h>
47 #endif
48
49
50
51 #define LEN 300 /* Max length of strings */
52 #define MASK32 0xffffffffUL /* 2^32 - 1 */
53
54
55
56
57
58 /*================================= Types =================================*/
59
60 typedef struct {
61 int kind;
62 long a1, q1, r1, a2, q2, r2;
63 long M;
64 double Norm;
65 } MRG2_param;
66
67 typedef struct {
68 long x1, x2;
69 } MRG2_state;
70
71 /*-------------------------------------------------------------------------*/
72
73 typedef struct {
74 int kind;
75 long a1, q1, r1, a3, q3, r3;
76 long M;
77 double Norm;
78 } MRG3_param;
79
80 typedef struct {
81 long x1, x2, x3;
82 } MRG3_state;
83
84 /*-------------------------------------------------------------------------*/
85
86 typedef struct {
87 int kind;
88 long a1, q1, r1, a5, q5, r5;
89 long M;
90 double Norm;
91 } MRG5_param;
92
93 typedef struct {
94 long x1, x2, x3, x4, x5;
95 } MRG5_state;
96
97 /*-------------------------------------------------------------------------*/
98
99 typedef struct {
100 int kind;
101 long a1, q1, r1, a7, q7, r7;
102 long M;
103 double Norm;
104 } MRG7_param;
105
106 typedef struct {
107 long x1, x2, x3, x4, x5, x6, x7;
108 } MRG7_state;
109
110 /*-------------------------------------------------------------------------*/
111
112 typedef struct {
113 int kind;
114 long *A, *Q, *R;
115 long M;
116 double Norm;
117 } MRG_param;
118
119 typedef struct {
120 long *S;
121 int k;
122 } MRG_state;
123
124 /*-------------------------------------------------------------------------*/
125
126 typedef struct {
127 double *A;
128 double M;
129 double Norm;
130 } MRGFloat_param;
131
132 typedef struct {
133 double *S;
134 int k;
135 } MRGFloat_state;
136
137 /*-------------------------------------------------------------------------*/
138
139 typedef struct {
140 long a12, a13, a21, a23; /* Multipliers */
141 long q12, q13, q21, q23; /* qij = mi DIV aj */
142 long r12, r13, r21, r23; /* rij = mi MOD aj */
143 long M1, M2; /* Modules */
144 double Norm;
145 } CombMRG3_param;
146
147 typedef struct {
148 long x10, x11, x12, x20, x21, x22; /* State */
149 } CombMRG3_state;
150
151 /*-------------------------------------------------------------------------*/
152 #ifdef USE_GMP
153
154 typedef struct {
155 mpz_t M, *A, W, T;
156 mpf_t F, Norm;
157 lebool *AnonZero;
158 } BigMRG_param;
159
160 typedef struct {
161 mpz_t *S;
162 int k;
163 } BigMRG_state;
164
165 typedef struct {
166 mpz_t M1, M2, *A1, *A2, W, T1, T2;
167 mpf_t F1, F2, Norm1, Norm2;
168 lebool *A1nonZero, *A2nonZero;
169 } BigC2MRG_param;
170
171 typedef struct {
172 mpz_t *S1, *S2;
173 int k;
174 } BigC2MRG_state;
175
176 #endif
177 /*-------------------------------------------------------------------------*/
178
179 typedef struct {
180 lebool Flag; /* TRUE if k > r, FALSE if k < r */
181 int Skip;
182 } LagFibFloat_param;
183
184 typedef struct {
185 double *X;
186 int r, s, RR;
187 int Lag;
188 } LagFibFloat_state;
189
190 /*-------------------------------------------------------------------------*/
191
192 typedef struct {
193 unsigned long Mask; /* 2^t - 1 */
194 int b; /* shift |t - 32| bits, right or left */
195 lebool LeftShift; /* TRUE for left shift, FALSE for right shift */
196 lebool Flag; /* TRUE if k > r, FALSE if k < r */
197 int Skip;
198 } LagFib_param;
199
200 typedef struct {
201 unsigned long *X;
202 int r, s, RR;
203 int Lag;
204 } LagFib_state;
205
206 /*-------------------------------------------------------------------------*/
207
208
209
210
211
212
213 /*============================ functions ==================================*/
214
215
AddArrayString(char * name,const char * add,int high,char * A[])216 static void AddArrayString (char *name, const char *add, int high, char *A[])
217 /*
218 * Add the string array A of size = high to name
219 */
220 {
221 int j;
222 strcat (name, add);
223 strcat (name, "(");
224 strcat (name, A[0]);
225 for (j = 1; (j < high) && (j < 7); j++) {
226 strcat (name, ", ");
227 if (A[j])
228 strcat (name, A[j]);
229 else
230 strcat (name, "0");
231 }
232 if (high > 7)
233 strcat (name, ", ... )");
234 else
235 strcat (name, ")");
236 }
237
238
239 /**************************************************************************/
240
MRG2_U01(void * vpar,void * vsta)241 static double MRG2_U01 (void *vpar, void *vsta)
242 /*
243 * Implementation used for k = 2, a1 > 0, a2 > 0
244 */
245 {
246 MRG2_param *param = vpar;
247 MRG2_state *state = vsta;
248 long h;
249 long p;
250
251 h = state->x2 / param->q2;
252 p = param->a2 * (state->x2 - h * param->q2) - h * param->r2;
253 if (p < 0)
254 p += param->M;
255 state->x2 = state->x1;
256
257 h = state->x1 / param->q1;
258 state->x1 = param->a1 * (state->x1 - h * param->q1) - h * param->r1;
259 if (state->x1 <= 0)
260 state->x1 += p;
261 else
262 state->x1 = (state->x1 - param->M) + p;
263 if (state->x1 < 0)
264 state->x1 += param->M;
265 return state->x1 * param->Norm;
266 }
267
268 /*-------------------------------------------------------------------------*/
269
MRG2_Bits(void * vpar,void * vsta)270 static unsigned long MRG2_Bits (void *vpar, void *vsta)
271 {
272 return (unsigned long) (unif01_NORM32 * MRG2_U01 (vpar, vsta));
273 }
274
275 /*-------------------------------------------------------------------------*/
276
WrMRG2(void * vsta)277 static void WrMRG2 (void *vsta)
278 {
279 MRG2_state *state = vsta;
280 printf (" S[1] = %10ld, S[2] = %10ld\n\n", state->x1, state->x2);
281 }
282
283
284 /**************************************************************************/
285
MRG3_U01(void * vpar,void * vsta)286 static double MRG3_U01 (void *vpar, void *vsta)
287 /*
288 * Implementation used for k = 3, a1 > 0, a3 > 0 and a2 = 0
289 */
290 {
291 MRG3_param *param = vpar;
292 MRG3_state *state = vsta;
293 long h;
294 long p;
295
296 h = state->x3 / param->q3;
297 p = param->a3 * (state->x3 - h * param->q3) - h * param->r3;
298 if (p < 0)
299 p += param->M;
300 state->x3 = state->x2;
301 state->x2 = state->x1;
302 h = state->x1 / param->q1;
303 state->x1 = param->a1 * (state->x1 - h * param->q1) - h * param->r1;
304 if (state->x1 <= 0)
305 state->x1 += p;
306 else
307 state->x1 = (state->x1 - param->M) + p;
308 if (state->x1 < 0)
309 state->x1 += param->M;
310 return state->x1 * param->Norm;
311 }
312
313 /*-------------------------------------------------------------------------*/
314
MRG3_Bits(void * vpar,void * vsta)315 static unsigned long MRG3_Bits (void *vpar, void *vsta)
316 {
317 return (unsigned long) (unif01_NORM32 * MRG3_U01 (vpar, vsta));
318 }
319
320 /*-------------------------------------------------------------------------*/
321
WrMRG3(void * vsta)322 static void WrMRG3 (void *vsta)
323 {
324 MRG3_state *state = vsta;
325 printf (" S[1] = %10ld, S[2] = %10ld, S[3] = %10ld\n\n",
326 state->x1, state->x2, state->x3);
327 }
328
329
330 /**************************************************************************/
331
MRG5_U01(void * vpar,void * vsta)332 static double MRG5_U01 (void *vpar, void *vsta)
333 /*
334 * Implementation used for k = 5, a1 > 0, a5 > 0 and a2 = a3 = a4 = 0
335 */
336 {
337 MRG5_param *param = vpar;
338 MRG5_state *state = vsta;
339 long h;
340 long p;
341
342 h = state->x5 / param->q5;
343 p = param->a5 * (state->x5 - h * param->q5) - h * param->r5;
344 if (p < 0)
345 p += param->M;
346 state->x5 = state->x4;
347 state->x4 = state->x3;
348 state->x3 = state->x2;
349 state->x2 = state->x1;
350 h = state->x1 / param->q1;
351 state->x1 = param->a1 * (state->x1 - h * param->q1) - h * param->r1;
352 if (state->x1 <= 0)
353 state->x1 += p;
354 else
355 state->x1 = (state->x1 - param->M) + p;
356 if (state->x1 < 0)
357 state->x1 += param->M;
358 return state->x1 * param->Norm;
359 }
360
361 /*-------------------------------------------------------------------------*/
362
MRG5_Bits(void * vpar,void * vsta)363 static unsigned long MRG5_Bits (void *vpar, void *vsta)
364 {
365 return (unsigned long) (unif01_NORM32 * MRG5_U01 (vpar, vsta));
366 }
367
368 /*-------------------------------------------------------------------------*/
369
WrMRG5(void * vsta)370 static void WrMRG5 (void *vsta)
371 {
372 MRG5_state *state = vsta;
373 printf (" S[1] = %10ld, S[2] = %10ld, S[3] = %10ld,\n",
374 state->x1, state->x2, state->x3);
375 printf (" S[4] = %10ld, S[5] = %10ld\n\n", state->x4, state->x5);
376 }
377
378
379 /**************************************************************************/
380
MRG7_U01(void * vpar,void * vsta)381 static double MRG7_U01 (void *vpar, void *vsta)
382 /*
383 * Implementation used for k = 7, a1 > 0, a7 > 0, a2 = a3 = a4 = a5 = a6 = 0
384 */
385 {
386 MRG7_param *param = vpar;
387 MRG7_state *state = vsta;
388 long h;
389 long p;
390
391 h = state->x7 / param->q7;
392 p = param->a7 * (state->x7 - h * param->q7) - h * param->r7;
393 if (p < 0)
394 p += param->M;
395 state->x7 = state->x6;
396 state->x6 = state->x5;
397 state->x5 = state->x4;
398 state->x4 = state->x3;
399 state->x3 = state->x2;
400 state->x2 = state->x1;
401 h = state->x1 / param->q1;
402 state->x1 = param->a1 * (state->x1 - h * param->q1) - h * param->r1;
403 if (state->x1 <= 0)
404 state->x1 += p;
405 else
406 state->x1 = (state->x1 - param->M) + p;
407 if (state->x1 < 0)
408 state->x1 += param->M;
409 return state->x1 * param->Norm;
410 }
411
412 /*-------------------------------------------------------------------------*/
413
MRG7_Bits(void * vpar,void * vsta)414 static unsigned long MRG7_Bits (void *vpar, void *vsta)
415 {
416 return (unsigned long) (unif01_NORM32 * MRG7_U01 (vpar, vsta));
417 }
418
419 /*-------------------------------------------------------------------------*/
420
WrMRG7(void * vsta)421 static void WrMRG7 (void *vsta)
422 {
423 MRG7_state *state = vsta;
424 printf (" S[1] = %10ld, S[2] = %10ld, S[3] = %10ld,\n",
425 state->x1, state->x2, state->x3);
426 printf (" S[4] = %10ld, S[5] = %10ld, S[6] = %10ld,\n",
427 state->x4, state->x5, state->x6);
428 printf (" S[7] = %10ld\n\n", state->x7);
429 }
430
431
432 /**************************************************************************/
433
MRG_U01(void * vpar,void * vsta)434 static double MRG_U01 (void *vpar, void *vsta)
435 /*
436 * Generator MRG of order k when R[i] < Q[i] for all i. The values returned
437 * by the generator MRG are in the set { 0, 1/m, 2/m, ..., (m-1)/m }
438 */
439 {
440 MRG_param *param = vpar;
441 MRG_state *state = vsta;
442 long i, p, h, t;
443
444 p = 0;
445 for (i = state->k; i > 0; i--) {
446 if (param->A[i] != 0) {
447 h = state->S[i] / param->Q[i];
448 t = labs (param->A[i]) * (state->S[i] - h * param->Q[i]) -
449 h * param->R[i];
450 if (t < 0)
451 t += param->M;
452 if (param->A[i] < 0)
453 p -= t;
454 else
455 p += (t - param->M);
456 if (p < 0)
457 p += param->M;
458 }
459 if (i > 1)
460 state->S[i] = state->S[i - 1];
461 else
462 state->S[i] = p;
463 }
464 return (p * param->Norm);
465 }
466
467 /*-------------------------------------------------------------------------*/
468
MRG_Bits(void * vpar,void * vsta)469 static unsigned long MRG_Bits (void *vpar, void *vsta)
470 {
471 return (unsigned long) (unif01_NORM32 * MRG_U01 (vpar, vsta));
472 }
473
474 /*-------------------------------------------------------------------------*/
475
WrMRG(void * vsta)476 static void WrMRG (void *vsta)
477 {
478 MRG_state *state = vsta;
479 int i;
480 if (unif01_WrLongStateFlag || (state->k < 8)) {
481 printf (" S = {\n ");
482 for (i = 1; i <= state->k; i++) {
483 printf (" %12ld", state->S[i]);
484 if (i < state->k)
485 printf (",");
486 if ((i % 4) == 0)
487 printf ("\n ");
488 }
489 printf (" }\n");
490 } else
491 unif01_WrLongStateDef ();
492 }
493
494 /*-------------------------------------------------------------------------*/
495
496 enum {MRG_ALL = 987654321};
497
CreateMRG_all(long m,int k,long AA[],long SS[])498 static unif01_Gen * CreateMRG_all (long m, int k, long AA[], long SS[])
499 {
500 unif01_Gen *gen;
501 MRG_param *param;
502 MRG_state *state;
503 size_t leng;
504 char name[LEN + 1];
505 int i;
506 long *A, *S, *Q, *R;
507
508 if ((k < 2) || (m <= 1))
509 util_Error ("umrg_CreateMRG: k < 2 or m < 2");
510
511 gen = util_Malloc (sizeof (unif01_Gen));
512 param = util_Malloc (sizeof (MRG_param));
513 state = util_Malloc (sizeof (MRG_state));
514
515 strncpy (name, "umrg_CreateMRG:", (size_t) LEN);
516 addstr_Long (name, " m = ", m);
517 addstr_Long (name, ", k = ", k);
518 addstr_ArrayLong (name, ", A = ", k, AA);
519 addstr_ArrayLong (name, ", S = ", k, SS);
520 leng = strlen (name);
521 gen->name = util_Calloc (leng + 1, sizeof (char));
522 strncpy (gen->name, name, leng);
523
524 A = util_Calloc ((size_t) k + 1, sizeof (long));
525 R = util_Calloc ((size_t) k + 1, sizeof (long));
526 Q = util_Calloc ((size_t) k + 1, sizeof (long));
527 S = util_Calloc ((size_t) k + 1, sizeof (long));
528
529 for (i = 1; i <= k; i++) {
530 A[i] = AA[i - 1];
531 S[i] = SS[i - 1];
532 if (A[i] != 0) {
533 R[i] = m % labs (A[i]);
534 Q[i] = m / labs (A[i]);
535 }
536 }
537
538 param->kind = MRG_ALL;
539 param->M = m;
540 param->Norm = 1.0 / m;
541 param->A = A;
542 param->R = R;
543 param->Q = Q;
544 state->k = k;
545 state->S = S;
546
547 gen->param = param;
548 gen->state = state;
549 gen->GetBits = &MRG_Bits;
550 gen->GetU01 = &MRG_U01;
551 gen->Write = &WrMRG;
552
553 return gen;
554 }
555
556 /*-------------------------------------------------------------------------*/
557
DeleteMRG_all(unif01_Gen * gen)558 static void DeleteMRG_all (unif01_Gen * gen)
559 {
560 MRG_param *param;
561 MRG_state *state;
562
563 if (NULL == gen)
564 return;
565 param = gen->param;
566 state = gen->state;
567 util_Free (state->S);
568 util_Free (param->A);
569 util_Free (param->Q);
570 util_Free (param->R);
571 gen->state = util_Free (gen->state);
572 gen->param = util_Free (gen->param);
573 gen->name = util_Free (gen->name);
574 util_Free (gen);
575 }
576
577
578 /**************************************************************************/
579
umrg_CreateMRG(long m,int k,long A[],long S[])580 unif01_Gen * umrg_CreateMRG (long m, int k, long A[], long S[])
581 {
582 unif01_Gen *gen;
583 size_t leng;
584 char name[LEN + 1];
585 int i, n;
586 long r, q;
587 long R[8] = { 0 };
588 long Q[8] = { 0 };
589
590 if ((k < 2) || (m <= 1))
591 util_Error ("umrg_CreateMRG: k < 2 or m < 2");
592
593 n = 0;
594 for (i = 0; i < k; i++) {
595 util_Assert (labs (A[i]) < m, "umrg_CreateMRG: |A[i]| >= m");
596 util_Assert (S[i] < m, "umrg_CreateMRG: S[i] >= m");
597 util_Assert (S[i] >= 0, "umrg_CreateMRG: S[i] < 0");
598 if (A[i] != 0) {
599 r = m % labs (A[i]);
600 q = m / labs (A[i]);
601 util_Assert (r <= q, "umrg_CreateMRG: r[i] > q[i]");
602 if (i < 7) {
603 R[i] = r;
604 Q[i] = q;
605 }
606 }
607 if (S[i] != 0)
608 n++;
609 }
610 util_Assert (n > 0, "umrg_CreateMRG: all S[i] are 0");
611
612
613 if ((k == 2) && (A[0] > 0) && (A[1] > 0)) {
614 MRG2_param *param;
615 MRG2_state *state;
616
617 gen = util_Malloc (sizeof (unif01_Gen));
618 param = util_Malloc (sizeof (MRG2_param));
619 state = util_Malloc (sizeof (MRG2_state));
620 param->kind = 2;
621 param->M = m;
622 param->Norm = 1.0 / m;
623 param->a1 = A[0];
624 param->r1 = R[0];
625 param->q1 = Q[0];
626 param->a2 = A[1];
627 param->r2 = R[1];
628 param->q2 = Q[1];
629 state->x1 = S[0];
630 state->x2 = S[1];
631
632 gen->param = param;
633 gen->state = state;
634 gen->GetBits = &MRG2_Bits;
635 gen->GetU01 = &MRG2_U01;
636 gen->Write = &WrMRG2;
637
638 } else if ((k == 3) && (A[0] > 0) && (A[1] == 0) && (A[2] > 0)) {
639 MRG3_param *param;
640 MRG3_state *state;
641
642 gen = util_Malloc (sizeof (unif01_Gen));
643 param = util_Malloc (sizeof (MRG3_param));
644 state = util_Malloc (sizeof (MRG3_state));
645 param->kind = 3;
646 param->M = m;
647 param->Norm = 1.0 / m;
648 param->a1 = A[0];
649 param->r1 = R[0];
650 param->q1 = Q[0];
651 param->a3 = A[2];
652 param->r3 = R[2];
653 param->q3 = Q[2];
654 state->x1 = S[0];
655 state->x2 = S[1];
656 state->x3 = S[2];
657
658 gen->param = param;
659 gen->state = state;
660 gen->GetBits = &MRG3_Bits;
661 gen->GetU01 = &MRG3_U01;
662 gen->Write = &WrMRG3;
663
664 } else if ((k == 5) && (A[0] > 0) && (A[1] == 0) && (A[2] == 0) &&
665 (A[3] == 0) && (A[4] > 0)) {
666 MRG5_param *param;
667 MRG5_state *state;
668
669 gen = util_Malloc (sizeof (unif01_Gen));
670 param = util_Malloc (sizeof (MRG5_param));
671 state = util_Malloc (sizeof (MRG5_state));
672 param->kind = 5;
673 param->M = m;
674 param->Norm = 1.0 / m;
675 param->a1 = A[0];
676 param->r1 = R[0];
677 param->q1 = Q[0];
678 param->a5 = A[4];
679 param->r5 = R[4];
680 param->q5 = Q[4];
681 state->x1 = S[0];
682 state->x2 = S[1];
683 state->x3 = S[2];
684 state->x4 = S[3];
685 state->x5 = S[4];
686
687 gen->param = param;
688 gen->state = state;
689 gen->GetBits = &MRG5_Bits;
690 gen->GetU01 = &MRG5_U01;
691 gen->Write = &WrMRG5;
692
693 } else if ((k == 7) && (A[0] > 0) && (A[1] == 0) && (A[2] == 0) &&
694 (A[3] == 0) && (A[4] == 0) && (A[5] == 0) && (A[6] > 0)) {
695 MRG7_param *param;
696 MRG7_state *state;
697
698 gen = util_Malloc (sizeof (unif01_Gen));
699 param = util_Malloc (sizeof (MRG7_param));
700 state = util_Malloc (sizeof (MRG7_state));
701 param->kind = 7;
702 param->M = m;
703 param->Norm = 1.0 / m;
704 param->a1 = A[0];
705 param->r1 = R[0];
706 param->q1 = Q[0];
707 param->a7 = A[6];
708 param->r7 = R[6];
709 param->q7 = Q[6];
710 state->x1 = S[0];
711 state->x2 = S[1];
712 state->x3 = S[2];
713 state->x4 = S[3];
714 state->x5 = S[4];
715 state->x6 = S[5];
716 state->x7 = S[6];
717
718 gen->param = param;
719 gen->state = state;
720 gen->GetBits = &MRG7_Bits;
721 gen->GetU01 = &MRG7_U01;
722 gen->Write = &WrMRG7;
723
724 } else {
725 return CreateMRG_all (m, k, A, S);
726 }
727
728 strncpy (name, "umrg_CreateMRG:", (size_t) LEN);
729 addstr_Long (name, " m = ", m);
730 addstr_Long (name, ", k = ", k);
731 addstr_ArrayLong (name, ", A = ", k, A);
732 addstr_ArrayLong (name, ", S = ", k, S);
733 leng = strlen (name);
734 gen->name = util_Calloc (leng + 1, sizeof (char));
735 strncpy (gen->name, name, leng);
736
737 return gen;
738 }
739
740 /*-------------------------------------------------------------------------*/
741
umrg_DeleteMRG(unif01_Gen * gen)742 void umrg_DeleteMRG (unif01_Gen * gen)
743 {
744 MRG_param *param;
745 if (NULL == gen)
746 return;
747 param = (MRG_param *) gen->param;
748 if (param->kind == MRG_ALL)
749 DeleteMRG_all (gen);
750 else
751 unif01_DeleteGen (gen);
752 }
753
754
755 /**************************************************************************/
756
MRGFloat_U01(void * vpar,void * vsta)757 static double MRGFloat_U01 (void *vpar, void *vsta)
758 /*
759 * Generator MRG of order k. Implementation uses floating-point.
760 */
761 {
762 MRGFloat_param *param = vpar;
763 MRGFloat_state *state = vsta;
764 double p;
765 long j;
766 int i;
767
768 p = 0.0;
769 for (i = state->k; i > 0; i--) {
770 if (param->A[i] != 0.0)
771 p += param->A[i] * state->S[i];
772 if (i > 1)
773 state->S[i] = state->S[i - 1];
774 }
775 j = p / param->M;
776 if (p >= 0.0)
777 p -= j * param->M; /* p = sum ... mod m */
778 else {
779 p += (1 - j) * param->M;
780 /* When truncating, if p % M = 0, we will get p = M */
781 if (p >= param->M)
782 p -= param->M;
783 }
784 state->S[1] = p; /* Must be in {0,..., m-1} */
785 return (p * param->Norm);
786 }
787
788 /*-------------------------------------------------------------------------*/
789
MRGFloat_Bits(void * vpar,void * vsta)790 static unsigned long MRGFloat_Bits (void *vpar, void *vsta)
791 {
792 return (unsigned long) (unif01_NORM32 * MRGFloat_U01 (vpar, vsta));
793 }
794
795 /*-------------------------------------------------------------------------*/
796
WrMRGFloat(void * vsta)797 static void WrMRGFloat (void *vsta)
798 {
799 MRGFloat_state *state = vsta;
800 int i;
801 if (unif01_WrLongStateFlag || (state->k < 8)) {
802 printf (" S = {\n ");
803 for (i = 1; i <= state->k; i++) {
804 printf (" %12.0f", state->S[i]);
805 if (i < state->k)
806 printf (",");
807 if ((i % 4) == 0)
808 printf ("\n ");
809 }
810 printf (" }\n");
811 } else
812 unif01_WrLongStateDef ();
813 }
814
815 /*-------------------------------------------------------------------------*/
816
umrg_CreateMRGFloat(long m,int k,long AA[],long SS[])817 unif01_Gen * umrg_CreateMRGFloat (long m, int k, long AA[], long SS[])
818 {
819 unif01_Gen *gen;
820 MRGFloat_param *param;
821 MRGFloat_state *state;
822 size_t leng;
823 char name[LEN + 1];
824 int i, n;
825 double *A, *S;
826 double pr1, pr2;
827
828 if ((k < 2) || (m < 2))
829 util_Error ("umrg_CreateMRGFloat: k < 2 or m < 2");
830
831 gen = util_Malloc (sizeof (unif01_Gen));
832 param = util_Malloc (sizeof (MRGFloat_param));
833 state = util_Malloc (sizeof (MRGFloat_state));
834
835 strncpy (name, "umrg_CreateMRGFloat:", (size_t) LEN);
836 addstr_Long (name, " m = ", m);
837 addstr_Long (name, ", k = ", k);
838 addstr_ArrayLong (name, ", A = ", k, AA);
839 addstr_ArrayLong (name, ", S = ", k, SS);
840 leng = strlen (name);
841 gen->name = util_Calloc (leng + 1, sizeof (char));
842 strncpy (gen->name, name, leng);
843
844 A = util_Calloc ((size_t) k + 1, sizeof (double));
845 S = util_Calloc ((size_t) k + 1, sizeof (double));
846
847 n = 0;
848 pr1 = 0.0;
849 pr2 = 0.0;
850 for (i = 1; i <= k; i++) {
851 if ((AA[i - 1] >= m) || (-AA[i - 1] >= m))
852 util_Error ("umrg_CreateMRGFloat: |A[i]| >= m");
853 if ((SS[i - 1] >= m) || (-SS[i - 1] >= m))
854 util_Error ("umrg_CreateMRGFloat: |S[i]| >= m");
855 A[i] = AA[i - 1];
856 S[i] = SS[i - 1];
857 if (SS[i - 1] < 0)
858 S[i] += m;
859 if (AA[i - 1] < 0)
860 pr2 -= A[i];
861 else
862 pr1 += A[i];
863 if (SS[i - 1] != 0)
864 n++;
865 }
866 if (n == 0)
867 util_Error (" umrg_CreateMRGFloat: all S[i] are 0");
868 if ((pr1 * m >= num_TwoExp[53]) || (pr2 * m >= num_TwoExp[53]))
869 util_Error ("umrg_CreateMRGFloat: Condition on a_i not valid");
870
871 param->M = m;
872 param->Norm = 1.0 / m;
873 param->A = A;
874 state->k = k;
875 state->S = S;
876
877 gen->param = param;
878 gen->state = state;
879 gen->GetBits = &MRGFloat_Bits;
880 gen->GetU01 = &MRGFloat_U01;
881 gen->Write = &WrMRGFloat;
882
883 return gen;
884 }
885
886 /*-------------------------------------------------------------------------*/
887
umrg_DeleteMRGFloat(unif01_Gen * gen)888 void umrg_DeleteMRGFloat (unif01_Gen * gen)
889 {
890 MRGFloat_param *param;
891 MRGFloat_state *state;
892
893 if (NULL == gen)
894 return;
895 param = gen->param;
896 state = gen->state;
897 util_Free (state->S);
898 util_Free (param->A);
899 gen->state = util_Free (gen->state);
900 gen->param = util_Free (gen->param);
901 gen->name = util_Free (gen->name);
902 util_Free (gen);
903 }
904
905
906 /**************************************************************************/
907
CombMRG3_U01(void * vpar,void * vsta)908 static double CombMRG3_U01 (void *vpar, void *vsta)
909 /*
910 * Implementation used for k = 3. It assumes that
911 * a11 = 0, a12 > 0, a13 < 0 and a21 > 0, a22 = 0, a23 < 0.
912 */
913 {
914 CombMRG3_param *param = vpar;
915 CombMRG3_state *state = vsta;
916 long h, p12, p13, p21, p23, Z;
917
918 /* Component 1 */
919 h = state->x10 / param->q13;
920 p13 = param->a13 * (state->x10 - h * param->q13) - h * param->r13;
921 if (p13 < 0)
922 p13 += param->M1;
923 util_Assert (p13 >= 0,
924 "umrg_CreateC2MRG: invalid parameters for a_{1,3}");
925 h = state->x11 / param->q12;
926 p12 = param->a12 * (state->x11 - h * param->q12) - h * param->r12;
927 if (p12 < 0)
928 p12 += param->M1;
929 util_Assert (p12 >= 0,
930 "umrg_CreateC2MRG: invalid parameters for a_{1,2}");
931 state->x10 = state->x11;
932 state->x11 = state->x12;
933 state->x12 = p12 - p13;
934 if (state->x12 < 0)
935 state->x12 += param->M1;
936
937 /* Component 2 */
938 h = state->x20 / param->q23;
939 p23 = param->a23 * (state->x20 - h * param->q23) - h * param->r23;
940 if (p23 < 0)
941 p23 += param->M2;
942 util_Assert (p23 >= 0,
943 "umrg_CreateC2MRG: invalid parameters for a_{2,3}");
944 h = state->x22 / param->q21;
945 p21 = param->a21 * (state->x22 - h * param->q21) - h * param->r21;
946 if (p21 < 0)
947 p21 += param->M2;
948 util_Assert (p21 >= 0,
949 "umrg_CreateC2MRG: invalid parameters for a_{2,1}");
950 state->x20 = state->x21;
951 state->x21 = state->x22;
952 state->x22 = p21 - p23;
953 if (state->x22 < 0)
954 state->x22 += param->M2;
955
956 /* Combinaison */
957 if (state->x12 < state->x22)
958 Z = state->x12 - state->x22 + param->M1;
959 else
960 Z = state->x12 - state->x22;
961 return (Z * param->Norm);
962 }
963
964 /*-------------------------------------------------------------------------*/
965
CombMRG3_Bits(void * vpar,void * vsta)966 static unsigned long CombMRG3_Bits (void *vpar, void *vsta)
967 {
968 return (unsigned long) (unif01_NORM32 * CombMRG3_U01 (vpar, vsta));
969 }
970
971 /*-------------------------------------------------------------------------*/
972
WrCombMRG3(void * vsta)973 static void WrCombMRG3 (void *vsta)
974 {
975 CombMRG3_state *state = vsta;
976 printf (" S1[0] = %1ld S1[1] = %1ld S1[2] = %1ld\n",
977 state->x10, state->x11, state->x12);
978 printf (" S2[0] = %1ld S2[1] = %1ld S2[2] = %1ld\n\n",
979 state->x20, state->x21, state->x22);
980 }
981
982 /*-------------------------------------------------------------------------*/
983
umrg_CreateC2MRG(long m1,long m2,int k,long A1[],long A2[],long S1[],long S2[])984 unif01_Gen * umrg_CreateC2MRG (long m1, long m2, int k, long A1[], long A2[],
985 long S1[], long S2[])
986 {
987 unif01_Gen *gen;
988 CombMRG3_param *param;
989 CombMRG3_state *state;
990 size_t leng;
991 char name[LEN + 1];
992 int i;
993 long cs1[4], cs2[4], cq1[4] = {0}, cq2[4] = {0}, cr1[4] = {0},
994 cr2[4] = {0}, ca1[4], ca2[4];
995
996 if (k != 3)
997 util_Error ("umrg_CreateC2MRG: k != 3");
998
999 gen = util_Malloc (sizeof (unif01_Gen));
1000 param = util_Malloc (sizeof (CombMRG3_param));
1001 state = util_Malloc (sizeof (CombMRG3_state));
1002
1003 strncpy (name, "umrg_CreateC2MRG:", (size_t) LEN);
1004 addstr_Long (name, " m1 = ", m1);
1005 addstr_Long (name, ", m2 = ", m2);
1006 addstr_Long (name, ", k = ", k);
1007 addstr_ArrayLong (name, ", A1 = ", k, A1);
1008 addstr_ArrayLong (name, ", S1 = ", k, S1);
1009 addstr_ArrayLong (name, ", A2 = ", k, A2);
1010 addstr_ArrayLong (name, ", S2 = ", k, S2);
1011 leng = strlen (name);
1012 gen->name = util_Calloc (leng + 1, sizeof (char));
1013 strncpy (gen->name, name, leng);
1014
1015 for (i = 1; i <= k; i++) {
1016 ca1[i] = A1[i - 1];
1017 ca2[i] = A2[i - 1];
1018 cs1[i] = S1[i - 1];
1019 cs2[i] = S2[i - 1];
1020 if (ca1[i] != 0) {
1021 cr1[i] = m1 % labs (ca1[i]);
1022 cq1[i] = m1 / labs (ca1[i]);
1023 util_Assert ((double) cr1[i] * labs(ca1[i]) < m1,
1024 "umrg_CreateC2MRG: |A1[i]| * (m1 mod |A1[i]|) >= m1");
1025 }
1026 if (ca2[i] != 0) {
1027 cr2[i] = m2 % labs (ca2[i]);
1028 cq2[i] = m2 / labs (ca2[i]);
1029 util_Assert ((double) cr2[i] * labs(ca2[i]) < m2,
1030 "umrg_CreateC2MRG:pp |A2[i]| * (m2 mod |A2[i]|) >= m2");
1031 }
1032 }
1033 param->M1 = m1;
1034 param->M2 = m2;
1035 param->Norm = 1.0 / m1;
1036 if (k == 3) {
1037 param->a12 = ca1[2];
1038 param->a13 = ca1[3];
1039 param->a21 = ca2[1];
1040 param->a23 = ca2[3];
1041 param->q12 = cq1[2];
1042 param->q13 = cq1[3];
1043 param->q21 = cq2[1];
1044 param->q23 = cq2[3];
1045 param->r12 = cr1[2];
1046 param->r13 = cr1[3];
1047 param->r21 = cr2[1];
1048 param->r23 = cr2[3];
1049 state->x10 = cs1[1];
1050 state->x11 = cs1[2];
1051 state->x12 = cs1[3];
1052 state->x20 = cs2[1];
1053 state->x21 = cs2[2];
1054 state->x22 = cs2[3];
1055
1056 gen->GetBits = &CombMRG3_Bits;
1057 gen->GetU01 = &CombMRG3_U01;
1058 gen->Write = &WrCombMRG3;
1059 }
1060 gen->param = param;
1061 gen->state = state;
1062 return gen;
1063 }
1064
1065 /*-------------------------------------------------------------------------*/
1066
umrg_DeleteC2MRG(unif01_Gen * gen)1067 void umrg_DeleteC2MRG (unif01_Gen *gen)
1068 {
1069 unif01_DeleteGen (gen);
1070 }
1071
1072
1073 /***************************************************************************/
1074
1075 #ifdef USE_GMP
1076
WrBigMRG(void * vsta)1077 static void WrBigMRG (void *vsta)
1078 {
1079 BigMRG_state *state = vsta;
1080 int i;
1081
1082 printf (" S = {\n ");
1083 for (i = 1; i <= state->k; i++) {
1084 mpz_out_str (NULL, 10, state->S[i]);
1085 if (i < state->k)
1086 printf (",\n ");
1087 }
1088 printf ("\n }\n");
1089 }
1090
1091 /*-------------------------------------------------------------------------*/
1092
BigMRG_U01(void * vpar,void * vsta)1093 static double BigMRG_U01 (void *vpar, void *vsta)
1094 {
1095 BigMRG_param *param = vpar;
1096 BigMRG_state *state = vsta;
1097 int i;
1098
1099 mpz_set_ui (param->T, 0); /* T = 0 */
1100
1101 for (i = state->k; i >= 1; i--) {
1102 if (param->AnonZero[i]) {
1103 mpz_mul (param->W, param->A[i], state->S[i]); /* W = A[i]*S[i] */
1104 mpz_mod (param->W, param->W, param->M); /* W = W % M */
1105 mpz_add (param->T, param->W, param->T); /* T = W + T */
1106 }
1107 if (i > 1) {
1108 mpz_set (state->S[i], state->S[i - 1]); /* S[i] = S[i - 1] */
1109 } else {
1110 mpz_mod (state->S[i], param->T, param->M); /* S[i] = T % M */
1111 }
1112 }
1113 mpf_set_z (param->F, state->S[1]); /* F = S[1] */
1114 mpf_mul (param->F, param->F, param->Norm); /* F = F * Norm */
1115 return mpf_get_d (param->F); /* U = F */
1116 }
1117
1118 /*-------------------------------------------------------------------------*/
1119
BigMRG_Bits(void * vpar,void * vsta)1120 static unsigned long BigMRG_Bits (void *vpar, void *vsta)
1121 {
1122 return (unsigned long) (unif01_NORM32 * BigMRG_U01 (vpar, vsta));
1123 }
1124
1125 /*-------------------------------------------------------------------------*/
1126
umrg_CreateBigMRG(char * m,int k,char * A[],char * S[])1127 unif01_Gen *umrg_CreateBigMRG (char *m, int k, char *A[], char *S[])
1128 {
1129 unif01_Gen *gen;
1130 BigMRG_param *param;
1131 BigMRG_state *state;
1132 size_t len1, len2;
1133 char name[LEN + 1];
1134 int i, flag;
1135
1136 gen = util_Malloc (sizeof (unif01_Gen));
1137 param = util_Malloc (sizeof (BigMRG_param));
1138 state = util_Malloc (sizeof (BigMRG_state));
1139
1140 /* These flags are set to 0 if the corresponding element is 0 */
1141 param->AnonZero = util_Calloc ((size_t) k + 1, sizeof (lebool));
1142
1143 strncpy (name, "umrg_CreateBigMRG:", (size_t) LEN);
1144 strcat (name, " m = ");
1145 strcat (name, m);
1146 addstr_Long (name, ", k = ", k);
1147 len1 = strlen (name);
1148 len2 = len1 + 4 * strlen (", A = ");
1149
1150 for (i = 0; i < k; i++) {
1151 if (A[i]) {
1152 len2 += strlen (A[i]);
1153 param->AnonZero[i + 1] = TRUE;
1154 } else
1155 /* If NULL pointer, set corresponding element to 0 */
1156 param->AnonZero[i + 1] = FALSE;
1157
1158 if (S[i])
1159 len2 += strlen (S[i]);
1160 }
1161
1162 gen->name = util_Calloc (len2 + LEN, sizeof (char));
1163 strncpy (gen->name, name, len2);
1164 AddArrayString (gen->name, ", A = ", k, A);
1165 AddArrayString (gen->name, ", S = ", k, S);
1166 len1 = strlen (gen->name);
1167 gen->name = util_Realloc (gen->name, (len1 + 1) * sizeof (char));
1168
1169 state->k = k;
1170 mpz_init (param->M);
1171 mpz_init (param->T);
1172 mpz_init (param->W);
1173 mpz_set_str (param->M, m, 0);
1174
1175 param->A = util_Calloc ((size_t) k + 1, sizeof (mpz_t));
1176 state->S = util_Calloc ((size_t) k + 1, sizeof (mpz_t));
1177
1178 for (i = 1; i <= k; i++) {
1179 if (param->AnonZero[i]) {
1180 mpz_init (param->A[i]);
1181 mpz_set_str (param->A[i], A[i - 1], 0);
1182 mpz_abs (param->T, param->A[i]);
1183 flag = mpz_cmp (param->T, param->M);
1184 util_Assert (flag < 0, "umrg_CreateBigMRG: one A >= m");
1185 /* Check if element is 0 */
1186 if (mpz_sgn (param->A[i]) == 0) {
1187 param->AnonZero[i] = FALSE;
1188 mpz_clear (param->A[i]);
1189 }
1190 }
1191 mpz_init (state->S[i]);
1192 if (S[i - 1]) {
1193 mpz_set_str (state->S[i], S[i - 1], 0);
1194 mpz_abs (param->T, state->S[i]);
1195 flag = mpz_cmp (param->T, param->M);
1196 util_Assert (flag < 0, "umrg_CreateBigMRG: one S >= m");
1197 }
1198 }
1199
1200 mpf_set_default_prec (DBL_MANT_DIG); /* Set precision to 53 bits */
1201 mpf_init (param->F);
1202 mpf_init (param->Norm);
1203 mpf_set_z (param->Norm, param->M);
1204 mpf_ui_div (param->Norm, 1, param->Norm); /* Norm = 1 / m */
1205
1206 gen->param = param;
1207 gen->state = state;
1208 gen->GetBits = &BigMRG_Bits;
1209 gen->GetU01 = &BigMRG_U01;
1210 gen->Write = &WrBigMRG;
1211
1212 return gen;
1213 }
1214
1215 /*-------------------------------------------------------------------------*/
1216
umrg_DeleteBigMRG(unif01_Gen * gen)1217 void umrg_DeleteBigMRG (unif01_Gen * gen)
1218 {
1219 BigMRG_param *param;
1220 BigMRG_state *state;
1221 int i;
1222
1223 if (NULL == gen)
1224 return;
1225 param = gen->param;
1226 state = gen->state;
1227 mpf_clear (param->F);
1228 mpf_clear (param->Norm);
1229
1230 for (i = 1; i <= state->k; i++) {
1231 if (param->AnonZero[i])
1232 mpz_clear (param->A[i]);
1233 mpz_clear (state->S[i]);
1234 }
1235
1236 mpz_clear (param->M);
1237 mpz_clear (param->T);
1238 mpz_clear (param->W);
1239
1240 util_Free (param->AnonZero);
1241 util_Free (state->S);
1242 util_Free (param->A);
1243
1244 gen->state = util_Free (gen->state);
1245 gen->param = util_Free (gen->param);
1246 gen->name = util_Free (gen->name);
1247 util_Free (gen);
1248 }
1249
1250
1251 /***************************************************************************/
1252
WrBigC2MRG(void * vsta)1253 static void WrBigC2MRG (void *vsta)
1254 {
1255 BigC2MRG_state *state = vsta;
1256 int i;
1257
1258 printf (" S1 = {\n ");
1259 for (i = 1; i <= state->k; i++) {
1260 mpz_out_str (NULL, 10, state->S1[i]);
1261 if (i < state->k)
1262 printf (",\n ");
1263 }
1264 printf ("\n }\n\n S2 = {\n ");
1265 for (i = 1; i <= state->k; i++) {
1266 mpz_out_str (NULL, 10, state->S2[i]);
1267 if (i < state->k)
1268 printf (",\n ");
1269 else
1270 printf ("\n }\n");
1271 }
1272 }
1273
1274 /*-------------------------------------------------------------------------*/
1275
BigC2MRG_U01(void * vpar,void * vsta)1276 static double BigC2MRG_U01 (void *vpar, void *vsta)
1277 {
1278 BigC2MRG_param *param = vpar;
1279 BigC2MRG_state *state = vsta;
1280 int i;
1281 double U, U2;
1282
1283 mpz_set_ui (param->T1, 0); /* T1 = 0 */
1284 mpz_set_ui (param->T2, 0); /* T2 = 0 */
1285
1286 for (i = state->k; i >= 1; i--) {
1287 if (param->A1nonZero[i]) {
1288 mpz_mul (param->W, param->A1[i], state->S1[i]); /* W = A1[i]*S1[i] */
1289 mpz_mod (param->W, param->W, param->M1); /* W = W % M1 */
1290 mpz_add (param->T1, param->W, param->T1); /* T1 = W + T1 */
1291 }
1292 if (param->A2nonZero[i]) {
1293 mpz_mul (param->W, param->A2[i], state->S2[i]); /* W = A2[i]*S2[i] */
1294 mpz_mod (param->W, param->W, param->M2); /* W = W % M2 */
1295 mpz_add (param->T2, param->W, param->T2); /* T2 = W + T2 */
1296 }
1297 if (i > 1) {
1298 mpz_set (state->S1[i], state->S1[i - 1]); /* S1[i] = S1[i - 1] */
1299 mpz_set (state->S2[i], state->S2[i - 1]); /* S2[i] = S2[i - 1] */
1300 } else {
1301 mpz_mod (state->S1[i], param->T1, param->M1); /* S1[i] = T1 % M1 */
1302 mpz_mod (state->S2[i], param->T2, param->M2); /* S2[i] = T2 % M2 */
1303 }
1304 }
1305 mpf_set_z (param->F1, state->S1[1]); /* F1 = S1[1] */
1306 mpf_set_z (param->F2, state->S2[1]); /* F2 = S2[1] */
1307 mpf_mul (param->F1, param->F1, param->Norm1); /* F1 = F1 * Norm1 */
1308 mpf_mul (param->F2, param->F2, param->Norm2); /* F2 = F2 * Norm2 */
1309 U = mpf_get_d (param->F1); /* U = F1 */
1310 U2 = mpf_get_d (param->F2); /* U2 = F2 */
1311
1312 U = U - U2;
1313 if (U < 0.0)
1314 return U + 1.0;
1315 if (U < 1.0)
1316 return U;
1317 return U - 1.0;
1318 }
1319
1320 /*-------------------------------------------------------------------------*/
1321
BigC2MRG_Bits(void * vpar,void * vsta)1322 static unsigned long BigC2MRG_Bits (void *vpar, void *vsta)
1323 {
1324 return (unsigned long) (unif01_NORM32 * BigC2MRG_U01 (vpar, vsta));
1325 }
1326
1327 /*-------------------------------------------------------------------------*/
1328
umrg_CreateBigC2MRG(char * m1,char * m2,int k,char * A1[],char * A2[],char * S1[],char * S2[])1329 unif01_Gen *umrg_CreateBigC2MRG (char *m1, char *m2, int k,
1330 char *A1[], char *A2[], char *S1[], char *S2[])
1331 {
1332 unif01_Gen *gen;
1333 BigC2MRG_param *param;
1334 BigC2MRG_state *state;
1335 size_t len1, len2;
1336 char name[LEN + 1];
1337 int i, flag;
1338
1339 gen = util_Malloc (sizeof (unif01_Gen));
1340 param = util_Malloc (sizeof (BigC2MRG_param));
1341 state = util_Malloc (sizeof (BigC2MRG_state));
1342
1343 /* These flags are set to 0 if the corresponding element is 0 */
1344 param->A1nonZero = util_Calloc ((size_t) k + 1, sizeof (lebool));
1345 param->A2nonZero = util_Calloc ((size_t) k + 1, sizeof (lebool));
1346
1347 strncpy (name, "umrg_CreateBigC2MRG:", (size_t) LEN);
1348 strcat (name, " m1 = ");
1349 strcat (name, m1);
1350 strcat (name, ", m2 = ");
1351 strcat (name, m2);
1352 addstr_Long (name, ", k = ", k);
1353 len1 = strlen (name);
1354 len2 = len1 + 4 * strlen (", A1 = ");
1355
1356 for (i = 0; i < k; i++) {
1357 if (A1[i]) {
1358 len2 += strlen (A1[i]);
1359 param->A1nonZero[i + 1] = TRUE;
1360 } else
1361 /* If NULL pointer, set corresponding element to 0 */
1362 param->A1nonZero[i + 1] = FALSE;
1363
1364 if (A2[i]) {
1365 len2 += strlen (A2[i]);
1366 param->A2nonZero[i + 1] = TRUE;
1367 } else
1368 param->A2nonZero[i + 1] = FALSE;
1369
1370 if (S1[i])
1371 len2 += strlen (S1[i]);
1372
1373 if (S2[i])
1374 len2 += strlen (S1[i]);
1375 }
1376
1377 gen->name = util_Calloc (len2 + LEN, sizeof (char));
1378 strncpy (gen->name, name, len2);
1379 AddArrayString (gen->name, ", A1 = ", k, A1);
1380 AddArrayString (gen->name, ", A2 = ", k, A2);
1381 AddArrayString (gen->name, ", S1 = ", k, S1);
1382 AddArrayString (gen->name, ", S2 = ", k, S2);
1383 len1 = strlen (gen->name);
1384 gen->name = util_Realloc (gen->name, (len1 + 1) * sizeof (char));
1385
1386 state->k = k;
1387 mpz_init (param->M1);
1388 mpz_init (param->M2);
1389 mpz_init (param->T1);
1390 mpz_init (param->T2);
1391 mpz_init (param->W);
1392 mpz_set_str (param->M1, m1, 0);
1393 mpz_set_str (param->M2, m2, 0);
1394
1395 param->A1 = util_Calloc ((size_t) k + 1, sizeof (mpz_t));
1396 param->A2 = util_Calloc ((size_t) k + 1, sizeof (mpz_t));
1397 state->S1 = util_Calloc ((size_t) k + 1, sizeof (mpz_t));
1398 state->S2 = util_Calloc ((size_t) k + 1, sizeof (mpz_t));
1399
1400 for (i = 1; i <= k; i++) {
1401 if (param->A1nonZero[i]) {
1402 mpz_init (param->A1[i]);
1403 mpz_set_str (param->A1[i], A1[i - 1], 0);
1404 mpz_abs (param->T1, param->A1[i]);
1405 flag = mpz_cmp (param->T1, param->M1);
1406 util_Assert (flag < 0, "umrg_CreateBigC2MRG: one A1 >= m1");
1407 /* Check if element is 0 */
1408 if (mpz_sgn (param->A1[i]) == 0) {
1409 param->A1nonZero[i] = FALSE;
1410 mpz_clear (param->A1[i]);
1411 }
1412 }
1413 if (param->A2nonZero[i]) {
1414 mpz_init (param->A2[i]);
1415 mpz_set_str (param->A2[i], A2[i - 1], 0);
1416 mpz_abs (param->T1, param->A2[i]);
1417 flag = mpz_cmp (param->T1, param->M2);
1418 util_Assert (flag < 0, "umrg_CreateBigC2MRG: one A2 >= m2");
1419 if (mpz_sgn (param->A2[i]) == 0) {
1420 param->A2nonZero[i] = FALSE;
1421 mpz_clear (param->A2[i]);
1422 }
1423 }
1424 mpz_init (state->S1[i]);
1425 if (S1[i - 1]) {
1426 mpz_set_str (state->S1[i], S1[i - 1], 0);
1427 mpz_abs (param->T1, state->S1[i]);
1428 flag = mpz_cmp (param->T1, param->M1);
1429 util_Assert (flag < 0, "umrg_CreateBigC2MRG: one S1 >= m1");
1430 }
1431
1432 mpz_init (state->S2[i]);
1433 if (S2[i - 1]) {
1434 mpz_set_str (state->S2[i], S2[i - 1], 0);
1435 mpz_abs (param->T1, state->S2[i]);
1436 flag = mpz_cmp (param->T1, param->M2);
1437 util_Assert (flag < 0, "umrg_CreateBigC2MRG: one S2 >= m2");
1438 }
1439 }
1440
1441 mpf_set_default_prec (DBL_MANT_DIG); /* Set precision to 53 bits */
1442 mpf_init (param->F1);
1443 mpf_init (param->F2);
1444 mpf_init (param->Norm1);
1445 mpf_init (param->Norm2);
1446 mpf_set_z (param->Norm1, param->M1);
1447 mpf_set_z (param->Norm2, param->M2);
1448 mpf_ui_div (param->Norm1, 1, param->Norm1); /* Norm1 = 1 / m1 */
1449 mpf_ui_div (param->Norm2, 1, param->Norm2); /* Norm2 = 1 / m2 */
1450
1451 gen->param = param;
1452 gen->state = state;
1453 gen->GetBits = &BigC2MRG_Bits;
1454 gen->GetU01 = &BigC2MRG_U01;
1455 gen->Write = &WrBigC2MRG;
1456
1457 return gen;
1458 }
1459
1460 /*-------------------------------------------------------------------------*/
1461
umrg_DeleteBigC2MRG(unif01_Gen * gen)1462 void umrg_DeleteBigC2MRG (unif01_Gen * gen)
1463 {
1464 BigC2MRG_param *param;
1465 BigC2MRG_state *state;
1466 int i;
1467
1468 if (NULL == gen)
1469 return;
1470 param = gen->param;
1471 state = gen->state;
1472 mpf_clear (param->F1);
1473 mpf_clear (param->F2);
1474 mpf_clear (param->Norm1);
1475 mpf_clear (param->Norm2);
1476
1477 for (i = 1; i <= state->k; i++) {
1478 if (param->A1nonZero[i])
1479 mpz_clear (param->A1[i]);
1480 if (param->A2nonZero[i])
1481 mpz_clear (param->A2[i]);
1482 mpz_clear (state->S1[i]);
1483 mpz_clear (state->S2[i]);
1484 }
1485
1486 mpz_clear (param->M1);
1487 mpz_clear (param->M2);
1488 mpz_clear (param->T1);
1489 mpz_clear (param->T2);
1490 mpz_clear (param->W);
1491
1492 util_Free (param->A1nonZero);
1493 util_Free (param->A2nonZero);
1494
1495 util_Free (state->S1);
1496 util_Free (state->S2);
1497 util_Free (param->A1);
1498 util_Free (param->A2);
1499
1500 gen->state = util_Free (gen->state);
1501 gen->param = util_Free (gen->param);
1502 gen->name = util_Free (gen->name);
1503 util_Free (gen);
1504 }
1505
1506 #endif
1507
1508 /***************************************************************************/
1509
LagFibAddFloat_U01(void * junk,void * vsta)1510 static double LagFibAddFloat_U01 (void *junk, void *vsta)
1511 {
1512 LagFibFloat_state *state = vsta;
1513 double temp;
1514
1515 temp = state->X[state->r] + state->X[state->s];
1516 if (temp >= 1.0)
1517 temp -= 1.0;
1518 state->X[state->r] = temp;
1519 if (--state->r == 0)
1520 state->r = state->Lag;
1521 if (--state->s == 0)
1522 state->s = state->Lag;
1523 return temp;
1524 }
1525
1526 /*-------------------------------------------------------------------------*/
1527
LagFibSousFloat_U01(void * vpar,void * vsta)1528 static double LagFibSousFloat_U01 (void *vpar, void *vsta)
1529 {
1530 LagFibFloat_param *param = vpar;
1531 LagFibFloat_state *state = vsta;
1532 double temp;
1533
1534 if (param->Flag)
1535 temp = state->X[state->r] - state->X[state->s];
1536 else
1537 temp = state->X[state->s] - state->X[state->r];
1538 if (temp < 0.0)
1539 temp += 1.0;
1540 state->X[state->r] = temp;
1541 if (--state->r == 0)
1542 state->r = state->Lag;
1543 if (--state->s == 0)
1544 state->s = state->Lag;
1545 return temp;
1546 }
1547
1548 /*-------------------------------------------------------------------------*/
1549
LagFibAddFloatLux_U01(void * vpar,void * vsta)1550 static double LagFibAddFloatLux_U01 (void *vpar, void *vsta)
1551 {
1552 LagFibFloat_param *param = vpar;
1553 LagFibFloat_state *state = vsta;
1554 double temp;
1555 int j;
1556
1557 if (--state->RR == 0) {
1558 state->RR = state->Lag;
1559 for (j = 0; j < param->Skip; j++) {
1560 temp = state->X[state->r] + state->X[state->s];
1561 state->X[state->r] = (temp >= 1.0) ? temp - 1.0 : temp;
1562 if (--state->r == 0)
1563 state->r = state->Lag;
1564 if (--state->s == 0)
1565 state->s = state->Lag;
1566 }
1567 }
1568 temp = state->X[state->r] + state->X[state->s];
1569 if (temp >= 1.0)
1570 temp -= 1.0;
1571 state->X[state->r] = temp;
1572 if (--state->r == 0)
1573 state->r = state->Lag;
1574 if (--state->s == 0)
1575 state->s = state->Lag;
1576 return temp;
1577 }
1578
1579 /*-------------------------------------------------------------------------*/
1580
LagFibSousFloatLux_U01(void * vpar,void * vsta)1581 static double LagFibSousFloatLux_U01 (void *vpar, void *vsta)
1582 {
1583 LagFibFloat_param *param = vpar;
1584 LagFibFloat_state *state = vsta;
1585 double temp;
1586 int j;
1587
1588 if (--state->RR == 0) {
1589 state->RR = state->Lag;
1590 for (j = 0; j < param->Skip; j++) {
1591 if (param->Flag)
1592 temp = state->X[state->r] - state->X[state->s];
1593 else
1594 temp = state->X[state->s] - state->X[state->r];
1595 state->X[state->r] = (temp < 0.0) ? temp + 1.0 : temp;
1596 if (--state->r == 0)
1597 state->r = state->Lag;
1598 if (--state->s == 0)
1599 state->s = state->Lag;
1600 }
1601 }
1602 if (param->Flag)
1603 temp = state->X[state->r] - state->X[state->s];
1604 else
1605 temp = state->X[state->s] - state->X[state->r];
1606 if (temp < 0.0)
1607 temp += 1.0;
1608 state->X[state->r] = temp;
1609 if (--state->r == 0)
1610 state->r = state->Lag;
1611 if (--state->s == 0)
1612 state->s = state->Lag;
1613 return temp;
1614 }
1615
1616 /*-------------------------------------------------------------------------*/
1617
LagFibAddFloat_Bits(void * vpar,void * vsta)1618 static unsigned long LagFibAddFloat_Bits (void *vpar, void *vsta)
1619 {
1620 return (unsigned long) (unif01_NORM32 * LagFibAddFloat_U01 (vpar, vsta));
1621 }
1622
1623 /*-------------------------------------------------------------------------*/
1624
LagFibSousFloat_Bits(void * vpar,void * vsta)1625 static unsigned long LagFibSousFloat_Bits (void *vpar, void *vsta)
1626 {
1627 return (unsigned long) (unif01_NORM32 * LagFibSousFloat_U01 (vpar, vsta));
1628 }
1629
1630 /*-------------------------------------------------------------------------*/
1631
LagFibAddFloatLux_Bits(void * vpar,void * vsta)1632 static unsigned long LagFibAddFloatLux_Bits (void *vpar, void *vsta)
1633 {
1634 return (unsigned long) (unif01_NORM32 * LagFibAddFloatLux_U01 (vpar, vsta));
1635 }
1636
1637 /*-------------------------------------------------------------------------*/
1638
LagFibSousFloatLux_Bits(void * vpar,void * vsta)1639 static unsigned long LagFibSousFloatLux_Bits (void *vpar, void *vsta)
1640 {
1641 return (unsigned long) (unif01_NORM32 * LagFibSousFloatLux_U01 (vpar, vsta));
1642 }
1643
1644 /*-------------------------------------------------------------------------*/
1645
WrLagFibFloat(void * vsta)1646 static void WrLagFibFloat (void *vsta)
1647 {
1648 LagFibFloat_state *state = vsta;
1649 int j;
1650 unsigned long Z;
1651
1652 if (unif01_WrLongStateFlag) {
1653 printf ("S = {\n");
1654 for (j = 0; j < state->Lag; j++) {
1655 Z = unif01_NORM32 * state->X[state->r];
1656 printf (" %12lu", Z);
1657 if (--state->r == 0)
1658 state->r = state->Lag;
1659 if (j < state->Lag - 1)
1660 printf (",");
1661 if ((j % 5) == 4)
1662 printf ("\n");
1663 }
1664 printf (" }\n");
1665 } else
1666 unif01_WrLongStateDef ();
1667 }
1668
1669 /*-------------------------------------------------------------------------*/
1670
umrg_CreateLagFibFloat(int k,int r,char Op,int Lux,unsigned long S[])1671 unif01_Gen * umrg_CreateLagFibFloat (int k, int r, char Op, int Lux,
1672 unsigned long S[])
1673 {
1674 unif01_Gen *gen;
1675 LagFibFloat_param *param;
1676 LagFibFloat_state *state;
1677 size_t leng;
1678 char name[LEN + 1];
1679 char chaine[2];
1680 int j;
1681 double *X;
1682
1683 /* util_Assert (k > r, "umrg_CreateLagFibFloat: k <= r"); */
1684 util_Assert ((Op == '-') || (Op == '+'),
1685 "umrg_CreateLagFibFloat: only + and - are implemented");
1686
1687 gen = util_Malloc (sizeof (unif01_Gen));
1688 param = util_Malloc (sizeof (LagFibFloat_param));
1689 state = util_Malloc (sizeof (LagFibFloat_state));
1690
1691 strncpy (name, "umrg_CreateLagFibFloat:", (size_t) LEN);
1692 addstr_Long (name, " k = ", k);
1693 addstr_Long (name, ", r = ", r);
1694 strcat (name, ", Op = ");
1695 sprintf (chaine, "%c", Op);
1696 strcat (name, chaine);
1697 addstr_Long (name, ", Lux = ", Lux);
1698 if (k >= r)
1699 addstr_ArrayUlong (name, ", S = ", k, S);
1700 else
1701 addstr_ArrayUlong (name, ", S = ", r, S);
1702 leng = strlen (name);
1703 gen->name = util_Calloc (leng + 1, sizeof (char));
1704 strncpy (gen->name, name, leng);
1705
1706 if (k >= r) {
1707 state->Lag = k;
1708 state->r = k;
1709 state->s = r;
1710 param->Flag = TRUE;
1711 } else {
1712 state->Lag = r;
1713 state->r = r;
1714 state->s = k;
1715 param->Flag = FALSE;
1716 }
1717 param->Skip = Lux - state->Lag;
1718
1719 if (param->Skip > 0) {
1720 X = util_Calloc ((size_t) Lux + 1, sizeof (double));
1721 state->RR = state->Lag;
1722 switch (Op) {
1723 case '+':
1724 gen->GetBits = &LagFibAddFloatLux_Bits;
1725 gen->GetU01 = &LagFibAddFloatLux_U01;
1726 break;
1727 case '-':
1728 gen->GetBits = &LagFibSousFloatLux_Bits;
1729 gen->GetU01 = &LagFibSousFloatLux_U01;
1730 break;
1731 }
1732 } else {
1733 X = util_Calloc ((size_t) state->Lag + 1, sizeof (double));
1734 switch (Op) {
1735 case '+':
1736 gen->GetBits = &LagFibAddFloat_Bits;
1737 gen->GetU01 = &LagFibAddFloat_U01;
1738 break;
1739 case '-':
1740 gen->GetBits = &LagFibSousFloat_Bits;
1741 gen->GetU01 = &LagFibSousFloat_U01;
1742 break;
1743 }
1744 }
1745
1746 for (j = 0; j < state->Lag; j++)
1747 X[state->Lag - j] = (S[j] & MASK32) / unif01_NORM32;
1748
1749 state->X = X;
1750 gen->param = param;
1751 gen->state = state;
1752 gen->Write = &WrLagFibFloat;
1753
1754 return gen;
1755 }
1756
1757 /*-------------------------------------------------------------------------*/
1758
umrg_DeleteLagFibFloat(unif01_Gen * gen)1759 void umrg_DeleteLagFibFloat (unif01_Gen * gen)
1760 {
1761 LagFibFloat_state *state;
1762
1763 if (NULL == gen)
1764 return;
1765 state = gen->state;
1766 util_Free (state->X);
1767 gen->state = util_Free (gen->state);
1768 gen->param = util_Free (gen->param);
1769 gen->name = util_Free (gen->name);
1770 util_Free (gen);
1771 }
1772
1773
1774 /**************************************************************************/
1775
LagFibAdd_Bits(void * vpar,void * vsta)1776 static unsigned long LagFibAdd_Bits (void *vpar, void *vsta)
1777 {
1778 LagFib_param *param = vpar;
1779 LagFib_state *state = vsta;
1780 unsigned long temp;
1781
1782 state->X[state->r] = (state->X[state->r] + state->X[state->s]) & param->Mask;
1783 if (param->LeftShift)
1784 temp = state->X[state->r] << param->b;
1785 else
1786 temp = state->X[state->r] >> param->b;
1787 if (--state->r == 0)
1788 state->r = state->Lag;
1789 if (--state->s == 0)
1790 state->s = state->Lag;
1791 return temp;
1792 }
1793
1794 /*-------------------------------------------------------------------------*/
1795
LagFibAddLux_Bits(void * vpar,void * vsta)1796 static unsigned long LagFibAddLux_Bits (void *vpar, void *vsta)
1797 {
1798 LagFib_param *param = vpar;
1799 LagFib_state *state = vsta;
1800 unsigned long temp;
1801
1802 if (--state->RR == 0) {
1803 int j;
1804 state->RR = state->Lag;
1805 for (j = 0; j < param->Skip; j++) {
1806 state->X[state->r] = (state->X[state->r] + state->X[state->s]) & param->Mask;
1807 if (--state->r == 0)
1808 state->r = state->Lag;
1809 if (--state->s == 0)
1810 state->s = state->Lag;
1811 }
1812 }
1813 state->X[state->r] = (state->X[state->r] + state->X[state->s]) & param->Mask;
1814 if (param->LeftShift)
1815 temp = state->X[state->r] << param->b;
1816 else
1817 temp = state->X[state->r] >> param->b;
1818 if (--state->r == 0)
1819 state->r = state->Lag;
1820 if (--state->s == 0)
1821 state->s = state->Lag;
1822 return temp;
1823 }
1824
1825 /*-------------------------------------------------------------------------*/
1826
LagFibAdd_U01(void * vpar,void * vsta)1827 static double LagFibAdd_U01 (void *vpar, void *vsta)
1828 {
1829 return LagFibAdd_Bits (vpar, vsta) * unif01_INV32;
1830 }
1831
1832 /*-------------------------------------------------------------------------*/
1833
LagFibAddLux_U01(void * vpar,void * vsta)1834 static double LagFibAddLux_U01 (void *vpar, void *vsta)
1835 {
1836 return LagFibAddLux_Bits (vpar, vsta) * unif01_INV32;
1837 }
1838
1839
1840 /**************************************************************************/
1841
LagFibSub_Bits(void * vpar,void * vsta)1842 static unsigned long LagFibSub_Bits (void *vpar, void *vsta)
1843 {
1844 LagFib_param *param = vpar;
1845 LagFib_state *state = vsta;
1846 unsigned long temp;
1847
1848 if (param->Flag)
1849 state->X[state->r] = (state->X[state->r] - state->X[state->s]) & param->Mask;
1850 else
1851 state->X[state->r] = (state->X[state->s] - state->X[state->r]) & param->Mask;
1852 if (param->LeftShift)
1853 temp = state->X[state->r] << param->b;
1854 else
1855 temp = state->X[state->r] >> param->b;
1856 if (--state->r == 0)
1857 state->r = state->Lag;
1858 if (--state->s == 0)
1859 state->s = state->Lag;
1860 return temp;
1861 }
1862
1863 /*-------------------------------------------------------------------------*/
1864
LagFibSubLux_Bits(void * vpar,void * vsta)1865 static unsigned long LagFibSubLux_Bits (void *vpar, void *vsta)
1866 {
1867 LagFib_param *param = vpar;
1868 LagFib_state *state = vsta;
1869 unsigned long temp;
1870
1871 if (--state->RR == 0) {
1872 int j;
1873 state->RR = state->Lag;
1874 for (j = 0; j < param->Skip; j++) {
1875 if (param->Flag)
1876 state->X[state->r] = (state->X[state->r] - state->X[state->s]) &
1877 param->Mask;
1878 else
1879 state->X[state->r] = (state->X[state->s] - state->X[state->r]) &
1880 param->Mask;
1881 if (--state->r == 0)
1882 state->r = state->Lag;
1883 if (--state->s == 0)
1884 state->s = state->Lag;
1885 }
1886 }
1887 if (param->Flag)
1888 state->X[state->r] = (state->X[state->r] - state->X[state->s]) & param->Mask;
1889 else
1890 state->X[state->r] = (state->X[state->s] - state->X[state->r]) & param->Mask;
1891 if (param->LeftShift)
1892 temp = state->X[state->r] << param->b;
1893 else
1894 temp = state->X[state->r] >> param->b;
1895 if (--state->r == 0)
1896 state->r = state->Lag;
1897 if (--state->s == 0)
1898 state->s = state->Lag;
1899 return temp;
1900 }
1901
1902 /*-------------------------------------------------------------------------*/
1903
LagFibSub_U01(void * vpar,void * vsta)1904 static double LagFibSub_U01 (void *vpar, void *vsta)
1905 {
1906 return LagFibSub_Bits (vpar, vsta) * unif01_INV32;
1907 }
1908
1909 /*-------------------------------------------------------------------------*/
1910
LagFibSubLux_U01(void * vpar,void * vsta)1911 static double LagFibSubLux_U01 (void *vpar, void *vsta)
1912 {
1913 return LagFibSubLux_Bits (vpar, vsta) * unif01_INV32;
1914 }
1915
1916
1917 /**************************************************************************/
1918
LagFibXor_Bits(void * vpar,void * vsta)1919 static unsigned long LagFibXor_Bits (void *vpar, void *vsta)
1920 {
1921 LagFib_param *param = vpar;
1922 LagFib_state *state = vsta;
1923 unsigned long temp;
1924
1925 state->X[state->r] = state->X[state->s] ^ state->X[state->r];
1926 if (param->LeftShift)
1927 temp = state->X[state->r] << param->b;
1928 else
1929 temp = state->X[state->r] >> param->b;
1930 if (--state->r == 0)
1931 state->r = state->Lag;
1932 if (--state->s == 0)
1933 state->s = state->Lag;
1934 return temp;
1935 }
1936
1937 /*-------------------------------------------------------------------------*/
1938
LagFibXorLux_Bits(void * vpar,void * vsta)1939 static unsigned long LagFibXorLux_Bits (void *vpar, void *vsta)
1940 {
1941 LagFib_param *param = vpar;
1942 LagFib_state *state = vsta;
1943 unsigned long temp;
1944
1945 if (--state->RR == 0) {
1946 int j;
1947 state->RR = state->Lag;
1948 for (j = 0; j < param->Skip; j++) {
1949 state->X[state->r] = (state->X[state->r] ^ state->X[state->s]);
1950 if (--state->r == 0)
1951 state->r = state->Lag;
1952 if (--state->s == 0)
1953 state->s = state->Lag;
1954 }
1955 }
1956 state->X[state->r] = (state->X[state->r] ^ state->X[state->s]);
1957 if (param->LeftShift)
1958 temp = state->X[state->r] << param->b;
1959 else
1960 temp = state->X[state->r] >> param->b;
1961 if (--state->r == 0)
1962 state->r = state->Lag;
1963 if (--state->s == 0)
1964 state->s = state->Lag;
1965 return temp;
1966 }
1967
1968 /*-------------------------------------------------------------------------*/
1969
LagFibXor_U01(void * vpar,void * vsta)1970 static double LagFibXor_U01 (void *vpar, void *vsta)
1971 {
1972 return LagFibXor_Bits (vpar, vsta) * unif01_INV32;
1973 }
1974
1975 /*-------------------------------------------------------------------------*/
1976
LagFibXorLux_U01(void * vpar,void * vsta)1977 static double LagFibXorLux_U01 (void *vpar, void *vsta)
1978 {
1979 return LagFibXorLux_Bits (vpar, vsta) * unif01_INV32;
1980 }
1981
1982
1983 /**************************************************************************/
1984
LagFibMult_Bits(void * vpar,void * vsta)1985 static unsigned long LagFibMult_Bits (void *vpar, void *vsta)
1986 {
1987 LagFib_param *param = vpar;
1988 LagFib_state *state = vsta;
1989 unsigned long temp;
1990
1991 state->X[state->r] = (state->X[state->r] * state->X[state->s]) & param->Mask;
1992 if (param->LeftShift)
1993 temp = state->X[state->r] << param->b;
1994 else
1995 temp = state->X[state->r] >> param->b;
1996 if (--state->r == 0)
1997 state->r = state->Lag;
1998 if (--state->s == 0)
1999 state->s = state->Lag;
2000 return temp;
2001 }
2002
2003 /*-------------------------------------------------------------------------*/
2004
LagFibMultLux_Bits(void * vpar,void * vsta)2005 static unsigned long LagFibMultLux_Bits (void *vpar, void *vsta)
2006 {
2007 LagFib_param *param = vpar;
2008 LagFib_state *state = vsta;
2009 unsigned long temp;
2010
2011 if (--state->RR == 0) {
2012 int j;
2013 state->RR = state->Lag;
2014 for (j = 0; j < param->Skip; j++) {
2015 state->X[state->r] = (state->X[state->r] * state->X[state->s]) & param->Mask;
2016 if (--state->r == 0)
2017 state->r = state->Lag;
2018 if (--state->s == 0)
2019 state->s = state->Lag;
2020 }
2021 }
2022 state->X[state->r] = (state->X[state->r] * state->X[state->s]) & param->Mask;
2023 if (param->LeftShift)
2024 temp = state->X[state->r] << param->b;
2025 else
2026 temp = state->X[state->r] >> param->b;
2027 if (--state->r == 0)
2028 state->r = state->Lag;
2029 if (--state->s == 0)
2030 state->s = state->Lag;
2031 return temp;
2032 }
2033
2034 /*-------------------------------------------------------------------------*/
2035
LagFibMult_U01(void * vpar,void * vsta)2036 static double LagFibMult_U01 (void *vpar, void *vsta)
2037 {
2038 return LagFibMult_Bits (vpar, vsta) * unif01_INV32;
2039 }
2040
2041 /*-------------------------------------------------------------------------*/
2042
LagFibMultLux_U01(void * vpar,void * vsta)2043 static double LagFibMultLux_U01 (void *vpar, void *vsta)
2044 {
2045 return LagFibMultLux_Bits (vpar, vsta) * unif01_INV32;
2046 }
2047
2048 /*-------------------------------------------------------------------------*/
2049
WrLagFib(void * vsta)2050 static void WrLagFib (void *vsta)
2051 {
2052 LagFib_state *state = vsta;
2053 int j;
2054
2055 if (unif01_WrLongStateFlag) {
2056 printf ("S = {\n");
2057 for (j = 0; j < state->Lag; j++) {
2058 printf (" %12lu", state->X[state->r]);
2059 if (--state->r == 0)
2060 state->r = state->Lag;
2061 if (j < state->Lag - 1)
2062 printf (",");
2063 if ((j % 5) == 4)
2064 printf ("\n");
2065 }
2066 printf (" }\n");
2067 } else
2068 unif01_WrLongStateDef ();
2069 }
2070
2071 /*-------------------------------------------------------------------------*/
2072
umrg_CreateLagFib(int t,int k,int r,char Op,int Lux,unsigned long S[])2073 unif01_Gen * umrg_CreateLagFib (int t, int k, int r, char Op, int Lux,
2074 unsigned long S[])
2075 {
2076 unif01_Gen *gen;
2077 LagFib_param *param;
2078 LagFib_state *state;
2079 size_t leng;
2080 char name[LEN + 1];
2081 char chaine[2];
2082 int j;
2083 unsigned long *X;
2084
2085 util_Assert (t > 0, "umrg_CreateLagFib: t <= 0");
2086 #ifdef IS_ULONG32
2087 util_Assert (t <= 32, "umrg_CreateLagFib: t > 32");
2088 #else
2089 util_Assert (t <= 64, "umrg_CreateLagFib: t > 64");
2090 #endif
2091 /* util_Assert (k > r, "umrg_CreateLagFib: k <= r"); */
2092 util_Assert ((Op == '*') || (Op == '+') || (Op == '-') || (Op == 'x'),
2093 "umrg_CreateLagFib: Op must be one of { +, -, *, x }");
2094
2095 gen = util_Malloc (sizeof (unif01_Gen));
2096 param = util_Malloc (sizeof (LagFib_param));
2097 state = util_Malloc (sizeof (LagFib_state));
2098
2099 strcpy (name, "umrg_CreateLagFib:");
2100 addstr_Long (name, " t = ", t);
2101 addstr_Long (name, ", k = ", k);
2102 addstr_Long (name, ", r = ", r);
2103 strcat (name, ", Op = ");
2104 sprintf (chaine, "%c", Op);
2105 strcat (name, chaine);
2106 addstr_Long (name, ", Lux = ", Lux);
2107
2108 if (k >= r) {
2109 state->Lag = k;
2110 state->r = k;
2111 state->s = r;
2112 param->Flag = TRUE;
2113 } else {
2114 state->Lag = r;
2115 state->r = r;
2116 state->s = k;
2117 param->Flag = FALSE;
2118 }
2119 param->Skip = Lux - state->Lag;
2120
2121 addstr_ArrayUlong (name, ", S = ", state->Lag, S);
2122 leng = strlen (name);
2123 gen->name = util_Calloc (leng + 1, sizeof (char));
2124 strncpy (gen->name, name, leng);
2125
2126 param->b = t - 32;
2127 if (param->b <= 0) {
2128 param->b = -param->b;
2129 param->LeftShift = TRUE;
2130 } else {
2131 param->LeftShift = FALSE;
2132 }
2133 param->Mask = (1UL << t) - 1;
2134 #ifndef IS_ULONG32
2135 if (64 == t)
2136 param->Mask = 0xffffffffffffffffUL;
2137 #endif
2138 if (32 == t)
2139 param->Mask = 0xffffffffUL;
2140
2141 if (param->Skip > 0) {
2142 X = util_Calloc ((size_t) Lux + 1, sizeof (unsigned long));
2143 state->RR = state->Lag;
2144
2145 switch (Op) {
2146 case '*':
2147 gen->GetBits = &LagFibMultLux_Bits;
2148 gen->GetU01 = &LagFibMultLux_U01;
2149 break;
2150 case '-':
2151 gen->GetBits = &LagFibSubLux_Bits;
2152 gen->GetU01 = &LagFibSubLux_U01;
2153 break;
2154 case '+':
2155 gen->GetBits = &LagFibAddLux_Bits;
2156 gen->GetU01 = &LagFibAddLux_U01;
2157 break;
2158 case 'x':
2159 gen->GetBits = &LagFibXorLux_Bits;
2160 gen->GetU01 = &LagFibXorLux_U01;
2161 break;
2162 }
2163
2164 } else {
2165 X = util_Calloc ((size_t) state->Lag + 1, sizeof (unsigned long));
2166 switch (Op) {
2167 case '*':
2168 gen->GetBits = &LagFibMult_Bits;
2169 gen->GetU01 = &LagFibMult_U01;
2170 break;
2171 case '-':
2172 gen->GetBits = &LagFibSub_Bits;
2173 gen->GetU01 = &LagFibSub_U01;
2174 break;
2175 case '+':
2176 gen->GetBits = &LagFibAdd_Bits;
2177 gen->GetU01 = &LagFibAdd_U01;
2178 break;
2179 case 'x':
2180 gen->GetBits = &LagFibXor_Bits;
2181 gen->GetU01 = &LagFibXor_U01;
2182 break;
2183 }
2184 }
2185
2186 if (Op == '*') {
2187 int flag = 0;
2188 for (j = 0; j < state->Lag; j++) {
2189 X[state->Lag - j] = (S[j] & param->Mask) | 1;
2190 if (X[state->Lag - j] % 4 != 1)
2191 flag = 1;
2192 }
2193 /* Make sure that not all seeds are = 1 mod 4 */
2194 if (!flag)
2195 X[1] = (X[1] + 2) & param->Mask;
2196 } else {
2197 for (j = 0; j < state->Lag; j++)
2198 X[state->Lag - j] = S[j] & param->Mask;
2199 }
2200
2201 state->X = X;
2202 gen->param = param;
2203 gen->state = state;
2204 gen->Write = &WrLagFib;
2205 return gen;
2206 }
2207
2208
2209 /*-------------------------------------------------------------------------*/
2210
umrg_DeleteLagFib(unif01_Gen * gen)2211 void umrg_DeleteLagFib (unif01_Gen * gen)
2212 {
2213 LagFib_state *state;
2214
2215 if (NULL == gen)
2216 return;
2217 state = gen->state;
2218 util_Free (state->X);
2219 gen->state = util_Free (gen->state);
2220 gen->param = util_Free (gen->param);
2221 gen->name = util_Free (gen->name);
2222 util_Free (gen);
2223 }
2224
2225 /**************************************************************************/
2226