1 /*************************************************************************\
2 *
3 * Package: TestU01
4 * File: ucubic.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 "ucubic.h"
36 #include "unif01.h"
37
38 #include <stdio.h>
39 #include <stdlib.h>
40 #include <string.h>
41
42
43
44 #define LEN 150 /* Max length of strings */
45
46
47
48 /*================================ Types ================================*/
49
50 typedef struct {
51 unsigned long m, a, b, c, d;
52 double Norm;
53 } Cubic_param;
54
55 typedef struct {
56 unsigned long X;
57 } Cubic_state;
58
59
60 /*---------------*/
61 #ifdef USE_LONGLONG
62
63 typedef struct {
64 ulonglong m, a, b, c, d;
65 double Norm;
66 } CubicL_param;
67
68 typedef struct {
69 ulonglong X;
70 } CubicL_state;
71
72 /*---------------*/
73 #else
74
75 typedef struct {
76 long m, a, b, c, d;
77 double Norm;
78 } CubicL_param;
79
80 typedef struct {
81 long X;
82 } CubicL_state;
83
84 #endif
85 /*---------------*/
86
87
88 /*-------------------------------------------------------------------------*/
89
90 typedef struct {
91 double a, b, c, d, m;
92 double Norm;
93 } CubicFloat_param;
94
95 typedef struct {
96 double X;
97 } CubicFloat_state;
98
99
100 /*-------------------------------------------------------------------------*/
101
102 typedef struct {
103 unsigned long A, M;
104 double Norm;
105 } Cubic1_param;
106
107 typedef struct {
108 unsigned long X;
109 } Cubic1_state;
110
111
112 typedef struct {
113 long M, A;
114 double Norm;
115 } Cubic1L_param;
116
117 typedef struct {
118 long X;
119 } Cubic1L_state;
120
121
122 /*-------------------------------------------------------------------------*/
123
124 typedef struct {
125 double M, A;
126 double Norm;
127 } Cubic1Float_param;
128
129 typedef struct {
130 double X;
131 } Cubic1Float_state;
132
133
134 /*-------------------------------------------------------------------------*/
135
136 typedef struct {
137 unsigned long M1, M2, A1, A2;
138 double Norm1, Norm2;
139 } CombCubic2_param;
140
141 typedef struct {
142 unsigned long S1, S2;
143 } CombCubic2_state;
144
145 typedef struct {
146 long M1, M2, A1, A2;
147 double Norm1, Norm2;
148 } CombCubic2L_param;
149
150 typedef struct {
151 long S1, S2;
152 } CombCubic2L_state;
153
154
155 /*-------------------------------------------------------------------------*/
156
157 typedef struct {
158 long M, A, C;
159 double Norm;
160 } CubicOut_param;
161
162 typedef struct {
163 long X;
164 } CubicOut_state;
165
166
167 /*-------------------------------------------------------------------------*/
168
169 typedef struct {
170 double M, A, C;
171 double Norm;
172 } CubicOutFloat_param;
173
174 typedef struct {
175 double X;
176 } CubicOutFloat_state;
177
178
179
180
181 /**************************************************************************/
182
Cubic_U01(void * vpar,void * vsta)183 static double Cubic_U01 (void *vpar, void *vsta)
184 /*
185 * Implementation used when all intermediary results hold in an int of 32
186 * bits. This is the case when m < 2^{16} since a, b, c, d are all < m.
187 */
188 {
189 Cubic_param *param = vpar;
190 Cubic_state *state = vsta;
191
192 state->X = (param->d + state->X * ((param->c + state->X * ((param->b +
193 state->X * param->a) % param->m)) % param->m)) % param->m;
194 return (state->X * param->Norm);
195 }
196
197
198 /*---------------*/
199 #ifdef USE_LONGLONG
200
CubicL_U01(void * vpar,void * vsta)201 static double CubicL_U01 (void *vpar, void *vsta)
202 /*
203 * Implementation of Cubic used in the general case.
204 */
205 {
206 CubicL_param *param = vpar;
207 CubicL_state *state = vsta;
208
209 state->X = (param->d + state->X * ((param->c + state->X * ((param->b +
210 state->X * param->a) % param->m)) % param->m)) % param->m;
211 return state->X * param->Norm;
212 }
213
214 /*---------------*/
215 #else
216
CubicL_U01(void * vpar,void * vsta)217 static double CubicL_U01 (void *vpar, void *vsta)
218 /*
219 * Implementation of Cubic used in the general case. Very slow.
220 */
221 {
222 CubicL_param *param = vpar;
223 CubicL_state *state = vsta;
224 long k;
225
226 k = num_MultModL (param->a, state->X, param->b, param->m);
227 k = num_MultModL (k, state->X, param->c, param->m);
228 state->X = num_MultModL (k, state->X, param->d, param->m);
229 return state->X * param->Norm;
230 }
231
232 #endif
233 /*---------------*/
234
235
Cubic_Bits(void * vpar,void * vsta)236 static unsigned long Cubic_Bits (void *vpar, void *vsta)
237 {
238 return (unsigned long) (unif01_NORM32 * Cubic_U01 (vpar, vsta));
239 }
240
CubicL_Bits(void * vpar,void * vsta)241 static unsigned long CubicL_Bits (void *vpar, void *vsta)
242 {
243 return (unsigned long) (unif01_NORM32 * CubicL_U01 (vpar, vsta));
244 }
245
WrCubic(void * vsta)246 static void WrCubic (void *vsta)
247 {
248 Cubic_state *state = vsta;
249 printf (" X = %lu\n", state->X);
250 }
251
WrCubicL(void * vsta)252 static void WrCubicL (void *vsta)
253 {
254 CubicL_state *state = vsta;
255 printf (" X = %1lu\n", (unsigned long) state->X);
256 }
257
258
SetCubic_0(unif01_Gen * gen,long m,long a,long b,long c,long d,long s)259 static void SetCubic_0 (unif01_Gen * gen, long m, long a, long b,
260 long c, long d, long s)
261 {
262 Cubic_param *param;
263 Cubic_state *state;
264
265 param = util_Malloc (sizeof (Cubic_param));
266 state = util_Malloc (sizeof (Cubic_state));
267
268 param->Norm = 1.0 / m;
269 param->m = m;
270 param->a = a;
271 param->b = b;
272 param->c = c;
273 param->d = d;
274 state->X = s;
275
276 gen->GetBits = &Cubic_Bits;
277 gen->GetU01 = &Cubic_U01;
278 gen->Write = &WrCubic;
279 gen->param = param;
280 gen->state = state;
281 }
282
SetCubicL_0(unif01_Gen * gen,long m,long a,long b,long c,long d,long s)283 static void SetCubicL_0 (unif01_Gen * gen, long m, long a, long b,
284 long c, long d, long s)
285 {
286 CubicL_param *param;
287 CubicL_state *state;
288
289 param = util_Malloc (sizeof (CubicL_param));
290 state = util_Malloc (sizeof (CubicL_state));
291
292 param->Norm = 1.0 / m;
293 param->m = m;
294 param->a = a;
295 param->b = b;
296 param->c = c;
297 param->d = d;
298 state->X = s;
299
300 gen->GetBits = &CubicL_Bits;
301 gen->GetU01 = &CubicL_U01;
302 gen->Write = &WrCubicL;
303 gen->param = param;
304 gen->state = state;
305 }
306
ucubic_CreateCubic(long m,long a,long b,long c,long d,long s)307 unif01_Gen *ucubic_CreateCubic (long m, long a, long b, long c, long d,
308 long s)
309 {
310 unif01_Gen *gen;
311 size_t leng;
312 char name[LEN + 1];
313
314 util_Assert (m > 0, "ucubic_CreateCubic: m <= 0");
315 util_Assert ((a > 0) && (a < m),
316 "ucubic_CreateCubic: a must be in (0, m)");
317 util_Assert ((b >= 0) && (b < m),
318 "ucubic_CreateCubic: b must be in [0, m)");
319 util_Assert ((c >= 0) && (c < m),
320 "ucubic_CreateCubic: c must be in [0, m)");
321 util_Assert ((d >= 0) && (d < m),
322 "ucubic_CreateCubic: d must be in [0, m)");
323 util_Assert ((s >= 0) && (s < m),
324 "ucubic_CreateCubic: s must be in [0, m)");
325
326 gen = util_Malloc (sizeof (unif01_Gen));
327
328 strncpy (name, "ucubic_CreateCubic:", (size_t) LEN);
329 addstr_Long (name, " m = ", m);
330 addstr_Long (name, ", a = ", a);
331 addstr_Long (name, ", b = ", b);
332 addstr_Long (name, ", c = ", c);
333 addstr_Long (name, ", d = ", d);
334 addstr_Long (name, ", s = ", s);
335 leng = strlen (name);
336 gen->name = util_Calloc (leng + 1, sizeof (char));
337 strncpy (gen->name, name, leng);
338
339 if (m < num_TwoExp[16])
340 SetCubic_0 (gen, m, a, b, c, d, s);
341 else
342 SetCubicL_0 (gen, m, a, b, c, d, s);
343
344 return gen;
345 }
346
347
348 /**************************************************************************/
349
WrCubicFloat(void * vsta)350 static void WrCubicFloat (void *vsta)
351 {
352 CubicFloat_state *state = vsta;
353 printf (" X = %1ld\n", (long) state->X);
354 }
355
CubicFloatA_U01(void * vpar,void * vsta)356 static double CubicFloatA_U01 (void *vpar, void *vsta)
357 /*
358 * Implementation used when
359 * [ a*(m-1)^3 + b*(m-1)^2 + c*(m-1) + d ] / m < 2^31
360 * else k may not hold in a long int.
361 */
362 {
363 CubicFloat_param *param = vpar;
364 CubicFloat_state *state = vsta;
365 long k;
366
367 state->X = param->d + state->X * (param->c + state->X * (param->b +
368 state->X * param->a));
369 k = state->X * param->Norm;
370 state->X -= k * param->m;
371 return (state->X * param->Norm);
372 }
373
CubicFloatB_U01(void * vpar,void * vsta)374 static double CubicFloatB_U01 (void *vpar, void *vsta)
375 /*
376 * Implementation used when (m-1)*m < 2^53
377 * else, R*R + R may not hold in a double.
378 */
379 {
380 CubicFloat_param *param = vpar;
381 CubicFloat_state *state = vsta;
382 long k;
383 double y;
384
385 y = param->a * state->X + param->b;
386 k = y * param->Norm;
387 y -= k * param->m;
388 y = y * state->X + param->c;
389 k = y * param->Norm;
390 y -= k * param->m;
391 state->X = y * state->X + param->d;
392 k = state->X * param->Norm;
393 state->X -= k * param->m;
394 return state->X * param->Norm;
395 }
396
CubicFloatC_U01(void * vpar,void * vsta)397 static double CubicFloatC_U01 (void *vpar, void *vsta)
398 /*
399 * Implementation used in the general case.
400 */
401 {
402 CubicFloat_param *param = vpar;
403 CubicFloat_state *state = vsta;
404 double y;
405
406 y = num_MultModD (param->a, state->X, param->b, param->m);
407 y = num_MultModD (y, state->X, param->c, param->m);
408 state->X = num_MultModD (y, state->X, param->d, param->m);
409 return (state->X * param->Norm);
410 }
411
CubicFloatA_Bits(void * vpar,void * vsta)412 static unsigned long CubicFloatA_Bits (void *vpar, void *vsta)
413 {
414 return (unsigned long) (unif01_NORM32 * CubicFloatA_U01 (vpar, vsta));
415 }
416
CubicFloatB_Bits(void * vpar,void * vsta)417 static unsigned long CubicFloatB_Bits (void *vpar, void *vsta)
418 {
419 return (unsigned long) (unif01_NORM32 * CubicFloatB_U01 (vpar, vsta));
420 }
421
CubicFloatC_Bits(void * vpar,void * vsta)422 static unsigned long CubicFloatC_Bits (void *vpar, void *vsta)
423 {
424 return (unsigned long) (unif01_NORM32 * CubicFloatC_U01 (vpar, vsta));
425 }
426
427
ucubic_CreateCubicFloat(long m,long a,long b,long c,long d,long s)428 unif01_Gen *ucubic_CreateCubicFloat (long m, long a, long b, long c, long d,
429 long s)
430 {
431 unif01_Gen *gen;
432 CubicFloat_param *param;
433 CubicFloat_state *state;
434 size_t leng;
435 char name[LEN + 1];
436 double y, v;
437
438 util_Assert (m > 0, "ucubic_CreateCubicFloat: m <= 0");
439 util_Assert ((a > 0) && (a < m),
440 "ucubic_CreateCubicFloat: a must be in (0, m)");
441 util_Assert ((b >= 0) && (b < m),
442 "ucubic_CreateCubicFloat: b must be in [0, m)");
443 util_Assert ((c >= 0) && (c < m),
444 "ucubic_CreateCubicFloat: c must be in [0, m)");
445 util_Assert ((d >= 0) && (d < m),
446 "ucubic_CreateCubicFloat: d must be in [0, m)");
447 util_Assert ((s >= 0) && (s < m),
448 "ucubic_CreateCubicFloat: s must be in [0, m)");
449
450 gen = util_Malloc (sizeof (unif01_Gen));
451 param = util_Malloc (sizeof (CubicFloat_param));
452 state = util_Malloc (sizeof (CubicFloat_state));
453
454 strncpy (name, "ucubic_CreateCubicFloat:", (size_t) LEN);
455 addstr_Long (name, " m = ", m);
456 addstr_Long (name, ", a = ", a);
457 addstr_Long (name, ", b = ", b);
458 addstr_Long (name, ", c = ", c);
459 addstr_Long (name, ", d = ", d);
460 addstr_Long (name, ", s = ", s);
461 leng = strlen (name);
462 gen->name = util_Calloc (leng + 1, sizeof (char));
463 strncpy (gen->name, name, leng);
464
465 param->Norm = 1.0 / m;
466 param->m = m;
467 param->a = a;
468 param->b = b;
469 param->c = c;
470 param->d = d;
471 state->X = s;
472
473 gen->Write = &WrCubicFloat;
474 gen->param = param;
475 gen->state = state;
476
477 y = m - 1;
478 v = d + y * (c + y * (b + y * a));
479
480 if (v / m < num_TwoExp[31]) {
481 gen->GetU01 = CubicFloatA_U01;
482 gen->GetBits = CubicFloatA_Bits;
483
484 } else if (y * y < num_TwoExp[53]) {
485 gen->GetU01 = CubicFloatB_U01;
486 gen->GetBits = CubicFloatB_Bits;
487
488 } else {
489 gen->GetU01 = CubicFloatC_U01;
490 gen->GetBits = CubicFloatC_Bits;
491 }
492
493 return gen;
494 }
495
496
497 /**************************************************************************/
498
WrCubic1(void * vsta)499 static void WrCubic1 (void *vsta)
500 {
501 Cubic1_state *state = vsta;
502 printf (" X = %1lu\n", state->X);
503 }
504
Cubic1_U01(void * vpar,void * vsta)505 static double Cubic1_U01 (void *vpar, void *vsta)
506 /*
507 * Implementation used when all intermediary results hold in an int of 32
508 * bits. This is the case when M < 2^{16}.
509 */
510 {
511 Cubic1_param *param = vpar;
512 Cubic1_state *state = vsta;
513
514 state->X = (state->X * ((state->X * state->X) % param->M)) % param->M;
515 state->X = (param->A * state->X + 1) % param->M;
516 return (state->X * param->Norm);
517 }
518
Cubic1L_U01(void * vpar,void * vsta)519 static double Cubic1L_U01 (void *vpar, void *vsta)
520 /*
521 * Implementation of Cubic1 used in the general case. Very slow.
522 */
523 {
524 Cubic1L_param *param = vpar;
525 Cubic1L_state *state = vsta;
526 long S;
527
528 S = num_MultModL (state->X, state->X, 0L, param->M);
529 S = num_MultModL (state->X, S, 0L, param->M);
530 state->X = num_MultModL (param->A, S, 1L, param->M);
531 return (state->X * param->Norm);
532 }
533
Cubic1_Bits(void * vpar,void * vsta)534 static unsigned long Cubic1_Bits (void *vpar, void *vsta)
535 {
536 return (unsigned long) (unif01_NORM32 * Cubic1_U01 (vpar, vsta));
537 }
538
Cubic1L_Bits(void * vpar,void * vsta)539 static unsigned long Cubic1L_Bits (void *vpar, void *vsta)
540 {
541 return (unsigned long) (unif01_NORM32 * Cubic1L_U01 (vpar, vsta));
542 }
543
544
SetCubic1_0(unif01_Gen * gen,long m,long a,long s)545 static void SetCubic1_0 (unif01_Gen * gen, long m, long a, long s)
546 {
547 Cubic1_param *param;
548 Cubic1_state *state;
549
550 param = util_Malloc (sizeof (Cubic1_param));
551 state = util_Malloc (sizeof (Cubic1_state));
552
553 param->Norm = 1.0 / m;
554 param->M = m;
555 param->A = a;
556 state->X = s;
557
558 gen->GetBits = &Cubic1_Bits;
559 gen->GetU01 = &Cubic1_U01;
560 gen->Write = &WrCubic1;
561 gen->param = param;
562 gen->state = state;
563 }
564
SetCubic1L_0(unif01_Gen * gen,long m,long a,long s)565 static void SetCubic1L_0 (unif01_Gen * gen, long m, long a, long s)
566 {
567 Cubic1L_param *param;
568 Cubic1L_state *state;
569
570 param = util_Malloc (sizeof (Cubic1L_param));
571 state = util_Malloc (sizeof (Cubic1L_state));
572
573 param->Norm = 1.0 / m;
574 param->M = m;
575 param->A = a;
576 state->X = s;
577
578 gen->GetBits = &Cubic1L_Bits;
579 gen->GetU01 = &Cubic1L_U01;
580 gen->Write = &WrCubic1;
581 gen->param = param;
582 gen->state = state;
583 }
584
585
ucubic_CreateCubic1(long m,long a,long s)586 unif01_Gen *ucubic_CreateCubic1 (long m, long a, long s)
587 {
588 unif01_Gen *gen;
589 size_t leng;
590 char name[LEN + 1];
591
592 util_Assert ((m > 0), "ucubic_CreateCubic1: m <= 0");
593 util_Assert ((a > 0) && (a < m),
594 "ucubic_CreateCubic1: a must be in (0, m)");
595 util_Assert ((s >= 0) && (s < m),
596 "ucubic_CreateCubic1: s must be in [0, m)");
597
598 gen = util_Malloc (sizeof (unif01_Gen));
599
600 strncpy (name, "ucubic_CreateCubic1:", (size_t) LEN);
601 addstr_Long (name, " m = ", m);
602 addstr_Long (name, ", a = ", a);
603 addstr_Long (name, ", s = ", s);
604 leng = strlen (name);
605 gen->name = util_Calloc (leng + 1, sizeof (char));
606 strncpy (gen->name, name, leng);
607
608 if (m < num_TwoExp[16])
609 SetCubic1_0 (gen, m, a, s);
610 else
611 SetCubic1L_0 (gen, m, a, s);
612
613 return gen;
614 }
615
616
617 /**************************************************************************/
618
WrCubic1Float(void * vsta)619 static void WrCubic1Float (void *vsta)
620 {
621 Cubic1Float_state *state = vsta;
622 printf (" X = %1ld\n", (long) state->X);
623 }
624
Cubic1FloatA_U01(void * vpar,void * vsta)625 static double Cubic1FloatA_U01 (void *vpar, void *vsta)
626 /*
627 * Used when condition 1 + a*(m-1)^3/m < 2^31;
628 * else k may not hold in a long.
629 */
630 {
631 Cubic1Float_param *param = vpar;
632 Cubic1Float_state *state = vsta;
633 long k;
634
635 state->X = param->A * state->X * state->X * state->X + 1.0;
636 k = state->X * param->Norm;
637 state->X -= k * param->M;
638 return (state->X * param->Norm);
639 }
640
Cubic1FloatB_U01(void * vpar,void * vsta)641 static double Cubic1FloatB_U01 (void *vpar, void *vsta)
642 /*
643 * Used when condition (m-1)^2 < 2^53;
644 * else R*R may not hold in a double.
645 */
646 {
647 Cubic1Float_param *param = vpar;
648 Cubic1Float_state *state = vsta;
649 long k;
650 double y;
651
652 y = param->A * state->X;
653 k = y * param->Norm;
654 y -= k * param->M;
655 y *= state->X;
656 k = y * param->Norm;
657 y -= k * param->M;
658 state->X = y * state->X + 1.0;
659 k = state->X * param->Norm;
660 state->X -= k * param->M;
661 return (state->X * param->Norm);
662 }
663
Cubic1FloatC_U01(void * vpar,void * vsta)664 static double Cubic1FloatC_U01 (void *vpar, void *vsta)
665 /*
666 * Inmplementation used in the general case
667 */
668 {
669 Cubic1Float_param *param = vpar;
670 Cubic1Float_state *state = vsta;
671 double y;
672
673 y = num_MultModD (state->X, state->X, 0.0, param->M);
674 y = num_MultModD (state->X, y, 0.0, param->M);
675 state->X = num_MultModD (param->A, y, 1.0, param->M);
676 return (state->X * param->Norm);
677 }
678
Cubic1FloatA_Bits(void * vpar,void * vsta)679 static unsigned long Cubic1FloatA_Bits (void *vpar, void *vsta)
680 {
681 return (unsigned long) (unif01_NORM32 * Cubic1FloatA_U01 (vpar, vsta));
682 }
683
Cubic1FloatB_Bits(void * vpar,void * vsta)684 static unsigned long Cubic1FloatB_Bits (void *vpar, void *vsta)
685 {
686 return (unsigned long) (unif01_NORM32 * Cubic1FloatB_U01 (vpar, vsta));
687 }
688
Cubic1FloatC_Bits(void * vpar,void * vsta)689 static unsigned long Cubic1FloatC_Bits (void *vpar, void *vsta)
690 {
691 return (unsigned long) (unif01_NORM32 * Cubic1FloatC_U01 (vpar, vsta));
692 }
693
ucubic_CreateCubic1Float(long m,long a,long s)694 unif01_Gen *ucubic_CreateCubic1Float (long m, long a, long s)
695 {
696 unif01_Gen *gen;
697 Cubic1Float_param *param;
698 Cubic1Float_state *state;
699 size_t leng;
700 char name[LEN + 1];
701 double y;
702
703 util_Assert ((m > 0), "ucubic_CreateCubic1Float: m <= 0");
704 util_Assert ((a > 0) && (a < m),
705 "ucubic_CreateCubic1Float: a must be in (0, m)");
706 util_Assert ((s >= 0) && (s < m),
707 "ucubic_CreateCubic1Float: s must be in [0, m)");
708
709 gen = util_Malloc (sizeof (unif01_Gen));
710 param = util_Malloc (sizeof (Cubic1Float_param));
711 state = util_Malloc (sizeof (Cubic1Float_state));
712
713 strncpy (name, "ucubic_CreateCubic1Float:", (size_t) LEN);
714 addstr_Long (name, " m = ", m);
715 addstr_Long (name, ", a = ", a);
716 addstr_Long (name, ", s = ", s);
717 leng = strlen (name);
718 gen->name = util_Calloc (leng + 1, sizeof (char));
719 strncpy (gen->name, name, leng);
720
721 param->Norm = 1.0 / m;
722 param->M = m;
723 param->A = a;
724 state->X = s;
725
726 y = m - 1;
727 if (1.0 + a * y * y * y / m < num_TwoExp[31]) {
728 gen->GetU01 = Cubic1FloatA_U01;
729 gen->GetBits = Cubic1FloatA_Bits;
730
731 } else if (y * y < num_TwoExp[53]) {
732 gen->GetU01 = Cubic1FloatB_U01;
733 gen->GetBits = Cubic1FloatB_Bits;
734
735 } else {
736 gen->GetU01 = Cubic1FloatC_U01;
737 gen->GetBits = Cubic1FloatC_Bits;
738 }
739
740 gen->Write = &WrCubic1Float;
741 gen->param = param;
742 gen->state = state;
743
744 return gen;
745 }
746
747
748 /**************************************************************************/
749
WrCombCubic2(void * vsta)750 static void WrCombCubic2 (void *vsta)
751 {
752 CombCubic2_state *state = vsta;
753 printf (" X1 = %1lu, X2 = %1lu\n", state->S1, state->S2);
754 }
755
CombCubic2_U01(void * vpar,void * vsta)756 static double CombCubic2_U01 (void *vpar, void *vsta)
757 /*
758 * Implementation used when all intermediary results hold in an int of 32
759 * bits. This is the case when M < 2^{16}.
760 */
761 {
762 CombCubic2_param *param = vpar;
763 CombCubic2_state *state = vsta;
764 unsigned long Y;
765 double U;
766
767 Y = (state->S1 * ((state->S1 * state->S1) % param->M1)) % param->M1;
768 state->S1 = (param->A1 * Y + 1) % param->M1;
769 Y = (state->S2 * ((state->S2 * state->S2) % param->M2)) % param->M2;
770 state->S2 = (param->A2 * Y + 1) % param->M2;
771 U = state->S1 * param->Norm1 + state->S2 * param->Norm2;
772 if (U >= 1.0)
773 return (U - 1.0);
774 else
775 return U;
776 }
777
CombCubic2L_U01(void * vpar,void * vsta)778 static double CombCubic2L_U01 (void *vpar, void *vsta)
779 /*
780 * Implementation of CombCubic2 used in the general case. Very slow.
781 */
782 {
783 CombCubic2L_param *param = vpar;
784 CombCubic2L_state *state = vsta;
785 long Y;
786 double U;
787
788 Y = num_MultModL (state->S1, state->S1, 0L, param->M1);
789 Y = num_MultModL (state->S1, Y, 0L, param->M1);
790 state->S1 = num_MultModL (param->A1, Y, 1L, param->M1);
791 Y = num_MultModL (state->S2, state->S2, 0L, param->M2);
792 Y = num_MultModL (state->S2, Y, 0L, param->M2);
793 state->S2 = num_MultModL (param->A2, Y, 1L, param->M2);
794 U = state->S1 * param->Norm1 + state->S2 * param->Norm2;
795 if (U >= 1.0)
796 return (U - 1.0);
797 else
798 return U;
799 }
800
CombCubic2_Bits(void * vpar,void * vsta)801 static unsigned long CombCubic2_Bits (void *vpar, void *vsta)
802 {
803 return (unsigned long) (unif01_NORM32 * CombCubic2_U01 (vpar, vsta));
804 }
805
CombCubic2L_Bits(void * vpar,void * vsta)806 static unsigned long CombCubic2L_Bits (void *vpar, void *vsta)
807 {
808 return (unsigned long) (unif01_NORM32 * CombCubic2L_U01 (vpar, vsta));
809 }
810
811
SetCombCubic2_0(unif01_Gen * gen,long m1,long m2,long a1,long a2,long s1,long s2)812 static void SetCombCubic2_0 (unif01_Gen * gen, long m1, long m2, long a1,
813 long a2, long s1, long s2)
814 {
815 CombCubic2_param *param;
816 CombCubic2_state *state;
817
818 param = util_Malloc (sizeof (CombCubic2_param));
819 state = util_Malloc (sizeof (CombCubic2_state));
820
821 param->Norm1 = 1.0 / m1;
822 param->Norm2 = 1.0 / m2;
823 param->M1 = m1;
824 param->A1 = a1;
825 param->M2 = m2;
826 param->A2 = a2;
827 state->S1 = s1;
828 state->S2 = s2;
829
830 gen->GetBits = &CombCubic2_Bits;
831 gen->GetU01 = &CombCubic2_U01;
832 gen->Write = &WrCombCubic2;
833 gen->param = param;
834 gen->state = state;
835 }
836
SetCombCubic2L_0(unif01_Gen * gen,long m1,long m2,long a1,long a2,long s1,long s2)837 static void SetCombCubic2L_0 (unif01_Gen * gen, long m1, long m2, long a1,
838 long a2, long s1, long s2)
839 {
840 CombCubic2L_param *param;
841 CombCubic2L_state *state;
842
843 param = util_Malloc (sizeof (CombCubic2L_param));
844 state = util_Malloc (sizeof (CombCubic2L_state));
845
846 param->Norm1 = 1.0 / m1;
847 param->Norm2 = 1.0 / m2;
848 param->M1 = m1;
849 param->A1 = a1;
850 param->M2 = m2;
851 param->A2 = a2;
852 state->S1 = s1;
853 state->S2 = s2;
854
855 gen->GetBits = &CombCubic2L_Bits;
856 gen->GetU01 = &CombCubic2L_U01;
857 gen->Write = &WrCombCubic2;
858 gen->param = param;
859 gen->state = state;
860 }
861
862
ucubic_CreateCombCubic2(long m1,long m2,long a1,long a2,long s1,long s2)863 unif01_Gen *ucubic_CreateCombCubic2 (long m1, long m2, long a1, long a2,
864 long s1, long s2)
865 {
866 unif01_Gen *gen;
867 size_t leng;
868 char name[LEN + 1];
869
870 if ((a1 <= 0) || (a1 >= m1) || (s1 < 0) || (s1 >= m1) || (m1 <= 0) ||
871 (a2 <= 0) || (a2 >= m2) || (s2 < 0) || (s2 >= m2) || (m2 <= 0)) {
872 printf ("m1 = %1ld, m2 = %1ld, a1 = %1ld, a2 = %1ld,"
873 " s1 = %1ld, s2 = %1ld\n", m1, m2, a1, a2, s1, s2);
874 util_Error ("ucubic_CreateCombCubic2: Invalid parameter");
875 }
876
877 gen = util_Malloc (sizeof (unif01_Gen));
878
879 strncpy (name, "ucubic_CreateCombCubic2:", (size_t) LEN);
880 addstr_Long (name, " m1 = ", m1);
881 addstr_Long (name, ", m2 = ", m2);
882 addstr_Long (name, ", a1 = ", a1);
883 addstr_Long (name, ", a2 = ", a2);
884 addstr_Long (name, ", s1 = ", s1);
885 addstr_Long (name, ", s2 = ", s2);
886 leng = strlen (name);
887 gen->name = util_Calloc (leng + 1, sizeof (char));
888 strncpy (gen->name, name, leng);
889
890 if ((m1 < num_TwoExp[16]) && (m2 < num_TwoExp[16]))
891 SetCombCubic2_0 (gen, m1, m2, a1, a2, s1, s2);
892 else
893 SetCombCubic2L_0 (gen, m1, m2, a1, a2, s1, s2);
894 return gen;
895 }
896
897
898 /**************************************************************************/
899
WrCubicOut(void * vsta)900 static void WrCubicOut (void *vsta)
901 {
902 CubicOut_state *state = vsta;
903 printf (" X = %1ld\n", state->X);
904 }
905
CubicOut_U01(void * vpar,void * vsta)906 static double CubicOut_U01 (void *vpar, void *vsta)
907 {
908 CubicOut_param *param = vpar;
909 CubicOut_state *state = vsta;
910 long Y;
911
912 state->X = num_MultModL (param->A, state->X, param->C, param->M);
913 Y = num_MultModL (state->X, state->X, 0L, param->M);
914 Y = num_MultModL (state->X, Y, 0L, param->M);
915 return (Y * param->Norm);
916 }
917
CubicOut_Bits(void * vpar,void * vsta)918 static unsigned long CubicOut_Bits (void *vpar, void *vsta)
919 {
920 return (unsigned long) (unif01_NORM32 * CubicOut_U01 (vpar, vsta));
921 }
922
923
ucubic_CreateCubicOut(long m,long a,long c,long s)924 unif01_Gen *ucubic_CreateCubicOut (long m, long a, long c, long s)
925 {
926 unif01_Gen *gen;
927 CubicOut_param *param;
928 CubicOut_state *state;
929 size_t leng;
930 char name[LEN + 1];
931
932 util_Assert ((m > 0), "ucubic_CreateCubicOut: m <= 0");
933 util_Assert ((a > 0) && (a < m),
934 "ucubic_CreateCubicOut: a must be in (0, m)");
935 util_Assert ((c >= 0) && (c < m),
936 "ucubic_CreateCubicOut: c must be in [0, m)");
937 util_Assert ((s >= 0) && (s < m),
938 "ucubic_CreateCubicOut: s must be in [0, m)");
939
940 gen = util_Malloc (sizeof (unif01_Gen));
941 param = util_Malloc (sizeof (CubicOut_param));
942 state = util_Malloc (sizeof (CubicOut_state));
943
944 strncpy (name, "ucubic_CreateCubicOut:", (size_t) LEN);
945 addstr_Long (name, " m = ", m);
946 addstr_Long (name, ", a = ", a);
947 addstr_Long (name, ", c = ", c);
948 addstr_Long (name, ", s = ", s);
949 leng = strlen (name);
950 gen->name = util_Calloc (leng + 1, sizeof (char));
951 strncpy (gen->name, name, leng);
952
953 param->Norm = 1.0 / m;
954 param->M = m;
955 param->A = a;
956 param->C = c;
957 state->X = s;
958
959 gen->GetU01 = CubicOut_U01;
960 gen->GetBits = CubicOut_Bits;
961 gen->Write = &WrCubicOut;
962 gen->param = param;
963 gen->state = state;
964 return gen;
965 }
966
967
968 /**************************************************************************/
969
WrCubicOutFloat(void * vsta)970 static void WrCubicOutFloat (void *vsta)
971 {
972 CubicOutFloat_state *state = vsta;
973 printf (" X = %1ld\n", (long) state->X);
974 }
975
CubicOutFloatA_U01(void * vpar,void * vsta)976 static double CubicOutFloatA_U01 (void *vpar, void *vsta)
977 /*
978 * Used when condition (m-1)^3/m < 2^31;
979 * else k may not hold in a long.
980 */
981 {
982 CubicOutFloat_param *param = vpar;
983 CubicOutFloat_state *state = vsta;
984 long k;
985 double y;
986
987 state->X = param->A * state->X + param->C;
988 k = state->X * param->Norm;
989 state->X -= k * param->M;
990 y = state->X * state->X * state->X;
991 k = y * param->Norm;
992 y -= k * param->M;
993 return (y * param->Norm);
994 }
995
CubicOutFloatB_U01(void * vpar,void * vsta)996 static double CubicOutFloatB_U01 (void *vpar, void *vsta)
997 /*
998 * Floating-point implementation for the case when m <= 94906266 ~ 2^26.5;
999 * follows the condition (m-1)^2 < 2^53; else, R*R may not hold in a double.
1000 */
1001 {
1002 CubicOutFloat_param *param = vpar;
1003 CubicOutFloat_state *state = vsta;
1004 long k;
1005 double y;
1006
1007 state->X = param->A * state->X + param->C;
1008 k = state->X * param->Norm;
1009 state->X -= k * param->M;
1010 y = state->X * state->X;
1011 k = y * param->Norm;
1012 y -= k * param->M;
1013 y = state->X * y;
1014 k = y * param->Norm;
1015 y -= k * param->M;
1016 return (y * param->Norm);
1017 }
1018
CubicOutFloatC_U01(void * vpar,void * vsta)1019 static double CubicOutFloatC_U01 (void *vpar, void *vsta)
1020 {
1021 CubicOutFloat_param *param = vpar;
1022 CubicOutFloat_state *state = vsta;
1023 double y;
1024
1025 state->X = num_MultModD (param->A, state->X, param->C, param->M);
1026 y = num_MultModD (state->X, state->X, 0.0, param->M);
1027 y = num_MultModD (state->X, y, 0.0, param->M);
1028 return (y * param->Norm);
1029 }
1030
CubicOutFloatA_Bits(void * vpar,void * vsta)1031 static unsigned long CubicOutFloatA_Bits (void *vpar, void *vsta)
1032 {
1033 return (unsigned long) (unif01_NORM32 * CubicOutFloatA_U01 (vpar, vsta));
1034 }
1035
CubicOutFloatB_Bits(void * vpar,void * vsta)1036 static unsigned long CubicOutFloatB_Bits (void *vpar, void *vsta)
1037 {
1038 return (unsigned long) (unif01_NORM32 * CubicOutFloatB_U01 (vpar, vsta));
1039 }
1040
CubicOutFloatC_Bits(void * vpar,void * vsta)1041 static unsigned long CubicOutFloatC_Bits (void *vpar, void *vsta)
1042 {
1043 return (unsigned long) (unif01_NORM32 * CubicOutFloatC_U01 (vpar, vsta));
1044 }
1045
ucubic_CreateCubicOutFloat(long m,long a,long c,long s)1046 unif01_Gen *ucubic_CreateCubicOutFloat (long m, long a, long c, long s)
1047 {
1048 unif01_Gen *gen;
1049 CubicOutFloat_param *param;
1050 CubicOutFloat_state *state;
1051 size_t leng;
1052 char name[LEN + 1];
1053 double y;
1054
1055 util_Assert ((m > 0), "ucubic_CreateCubicOutFloat: m <= 0");
1056 util_Assert ((a > 0) && (a < m),
1057 "ucubic_CreateCubicOutFloat: a must be in (0, m)");
1058 util_Assert ((c >= 0) && (c < m),
1059 "ucubic_CreateCubicOutFloat: c must be in [0, m)");
1060 util_Assert ((s >= 0) && (s < m),
1061 "ucubic_CreateCubicOutFloat: s must be in [0, m)");
1062
1063 gen = util_Malloc (sizeof (unif01_Gen));
1064 param = util_Malloc (sizeof (CubicOutFloat_param));
1065 state = util_Malloc (sizeof (CubicOutFloat_state));
1066
1067 strncpy (name, "ucubic_CreateCubicOutFloat:", (size_t) LEN);
1068 addstr_Long (name, " m = ", m);
1069 addstr_Long (name, ", a = ", a);
1070 addstr_Long (name, ", c = ", c);
1071 addstr_Long (name, ", s = ", s);
1072 leng = strlen (name);
1073 gen->name = util_Calloc (leng + 1, sizeof (char));
1074 strncpy (gen->name, name, leng);
1075
1076 param->M = m;
1077 param->A = a;
1078 state->X = s;
1079 param->C = c;
1080 param->Norm = 1.0 / m;
1081
1082 y = m - 1;
1083 if (y * y * y / m < num_TwoExp[31]) {
1084 gen->GetU01 = CubicOutFloatA_U01;
1085 gen->GetBits = CubicOutFloatA_Bits;
1086
1087 } else if (y * y < num_TwoExp[53]) {
1088 gen->GetU01 = CubicOutFloatB_U01;
1089 gen->GetBits = CubicOutFloatB_Bits;
1090
1091 } else {
1092 gen->GetU01 = CubicOutFloatC_U01;
1093 gen->GetBits = CubicOutFloatC_Bits;
1094 }
1095
1096 gen->Write = &WrCubicOutFloat;
1097 gen->param = param;
1098 gen->state = state;
1099
1100 return gen;
1101 }
1102
1103
1104 /**************************************************************************/
1105
ucubic_DeleteGen(unif01_Gen * gen)1106 void ucubic_DeleteGen (unif01_Gen * gen)
1107 {
1108 unif01_DeleteGen (gen);
1109 }
1110