1 /*************************************************************************\
2 *
3 * Package: TestU01
4 * File: unif01.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 "chrono.h"
35 #include "swrite.h"
36 #include "unif01.h"
37
38 #include <math.h>
39 #include <stdio.h>
40 #include <stdlib.h>
41 #include <string.h>
42
43
44 #define LEN0 500 /* Length of strings */
45 #define LEN1 100 /* Length of strings */
46
47 #define MASK32 0xffffffffUL /* 2^32 - 1 */
48
49
50
51
52 /*------------------------- extern variables ------------------------------*/
53
54 lebool unif01_WrLongStateFlag = FALSE;
55
56
57
58
59
60 /* ========================== functions ================================== */
61
unif01_WriteNameGen(unif01_Gen * gen)62 void unif01_WriteNameGen (unif01_Gen *gen)
63 {
64 if (gen->name)
65 printf ("%s\n\n", gen->name);
66 }
67
unif01_WriteState(unif01_Gen * gen)68 void unif01_WriteState (unif01_Gen *gen)
69 {
70 printf ("\nGenerator state:\n");
71 gen->Write (gen->state);
72 printf ("\n");
73 }
74
unif01_WrLongStateDef(void)75 void unif01_WrLongStateDef (void)
76 {
77 printf (" Not shown here ... takes too much space\n");
78 }
79
80
81 /**************************************************************************/
82
unif01_StripD(unif01_Gen * gen,int r)83 double unif01_StripD (unif01_Gen *gen, int r)
84 {
85 if (r == 0) {
86 return (gen->GetU01) (gen->param, gen->state);
87 } else {
88 double u = num_TwoExp[r] * (gen->GetU01) (gen->param, gen->state);
89 return (u - (long) u);
90 }
91 }
92
unif01_StripL(unif01_Gen * gen,int r,long d)93 long unif01_StripL (unif01_Gen *gen, int r, long d)
94 {
95 if (r == 0)
96 return (long) (d * gen->GetU01 (gen->param, gen->state));
97 else {
98 double u = num_TwoExp[r] * (gen->GetU01) (gen->param, gen->state);
99 return (long) (d * (u - (long) u));
100 }
101 }
102
unif01_StripB(unif01_Gen * gen,int r,int s)103 unsigned long unif01_StripB (unif01_Gen *gen, int r, int s)
104 {
105 if (r == 0) {
106 return gen->GetBits (gen->param, gen->state) >> (32 - s);
107 } else {
108 unsigned long u = gen->GetBits (gen->param, gen->state);
109 return ((u << r) & MASK32) >> (32 - s);
110 }
111 }
112
113
114 /*************************************************************************/
115
116 /* Dummy generator, always return 0. */
117
DummyGen_U01(void * junk1,void * junk2)118 static double DummyGen_U01 (void *junk1, void *junk2)
119 {
120 return 0.0;
121 }
122
DummyGen_Bits(void * junk1,void * junk2)123 static unsigned long DummyGen_Bits (void *junk1, void *junk2)
124 {
125 return 0;
126 }
127
WrDummyGen(void * junk)128 static void WrDummyGen (void *junk)
129 {
130 printf (" Empty Generator (no state)\n");
131 }
132
unif01_CreateDummyGen(void)133 unif01_Gen * unif01_CreateDummyGen (void)
134 {
135 unif01_Gen *gen;
136 size_t len;
137
138 gen = util_Malloc (sizeof (unif01_Gen));
139 len = strlen ("Dummy generator that always returns 0");
140 gen->name = util_Calloc (len + 1, sizeof (char));
141 strncpy (gen->name, "Dummy generator that always returns 0", len);
142 gen->param = NULL;
143 gen->state = NULL;
144 gen->Write = &WrDummyGen;
145 gen->GetBits = &DummyGen_Bits;
146 gen->GetU01 = &DummyGen_U01;
147 return gen;
148 }
149
unif01_DeleteDummyGen(unif01_Gen * gen)150 void unif01_DeleteDummyGen (unif01_Gen *gen)
151 {
152 if (NULL == gen) return;
153 gen->name = util_Free (gen->name);
154 util_Free (gen);
155 }
156
unif01_DeleteGen(unif01_Gen * gen)157 void unif01_DeleteGen (unif01_Gen *gen)
158 {
159 if (NULL == gen) return;
160 gen->state = util_Free (gen->state);
161 gen->param = util_Free (gen->param);
162 gen->name = util_Free (gen->name);
163 util_Free (gen);
164 }
165
166
167 /**************************************************************************/
168 /*
169 * The original generator is gen0. The position of the bit at which the
170 * increased precision is applied is s, counting from the most significant
171 * bit; v = 1 / 2^s.
172 */
173 typedef struct {
174 unif01_Gen *gen0;
175 double v;
176 int s;
177 } DoubleGen_param;
178
179
DoubleGen_U01(void * vpar,void * junk)180 static double DoubleGen_U01 (void *vpar, void *junk)
181 {
182 double U;
183 DoubleGen_param *paramD = vpar;
184 unif01_Gen *gen = paramD->gen0;
185
186 U = gen->GetU01 (gen->param, gen->state);
187 U += paramD->v * gen->GetU01 (gen->param, gen->state);
188 if (U < 1.0)
189 return U;
190 else
191 return U - 1.0;
192 }
193
DoubleGen_Bits(void * vpar,void * junk)194 static unsigned long DoubleGen_Bits (void *vpar, void *junk)
195 {
196 return (unsigned long) (unif01_NORM32 * DoubleGen_U01 (vpar, junk));
197 }
198
199
unif01_CreateDoubleGen2(unif01_Gen * gen,double v)200 unif01_Gen * unif01_CreateDoubleGen2 (unif01_Gen *gen, double v)
201 {
202 unif01_Gen *genD;
203 DoubleGen_param *paramD;
204 char *name;
205 char str[20];
206 size_t len, len2, len3;
207
208 util_Assert (v > 0.0, "unif01_CreateDoubleGen2: h <= 0");
209 util_Assert (v < 1.0, "unif01_CreateDoubleGen2: h >= 1");
210 genD = util_Malloc (sizeof (unif01_Gen));
211 paramD = util_Malloc (sizeof (DoubleGen_param));
212 paramD->s = -num_Log2(v);
213 paramD->v = v;
214 paramD->gen0 = gen;
215
216 len = strlen (gen->name);
217 len2 = strlen ("\nunif01_CreateDoubleGen2 with h = ");
218 len += len2;
219 sprintf (str, "%-g", v);
220 len3 = strlen (str);
221 len += len3;
222 name = util_Calloc (len + 1, sizeof (char));
223 strncpy (name, gen->name, len);
224 strncat (name, "\nunif01_CreateDoubleGen2 with h = ", len2);
225 strncat (name, str, len3);
226
227 /* The state of the double generator is simply the state of the original
228 generator */
229 genD->name = name;
230 genD->param = paramD;
231 genD->state = gen->state;
232 genD->Write = gen->Write;
233 genD->GetBits = &DoubleGen_Bits;
234 genD->GetU01 = &DoubleGen_U01;
235 return genD;
236 }
237
unif01_CreateDoubleGen(unif01_Gen * gen,int s)238 unif01_Gen * unif01_CreateDoubleGen (unif01_Gen *gen, int s)
239 {
240 unif01_Gen *genD;
241 DoubleGen_param *paramD;
242 char *name;
243 char str[8];
244 size_t len, len2, len3;
245
246 util_Assert (s > 0, "unif01_CreateDoubleGen: s <= 0");
247 genD = unif01_CreateDoubleGen2 (gen, 1.0 / num_TwoExp[s]);
248 paramD = genD->param;
249 paramD->s = s;
250
251 len = strlen (gen->name);
252 len2 = strlen ("\nunif01_CreateDoubleGen with s = ");
253 len += len2;
254 sprintf (str, "%-d", paramD->s);
255 len3 = strlen (str);
256 len += len3;
257 name = util_Calloc (len + 1, sizeof (char));
258 strncpy (name, gen->name, len);
259 strncat (name, "\nunif01_CreateDoubleGen with s = ", len2);
260 strncat (name, str, len3);
261 genD->name = name;
262 return genD;
263 }
264
unif01_DeleteDoubleGen(unif01_Gen * gen)265 void unif01_DeleteDoubleGen (unif01_Gen *gen)
266 {
267 if (NULL == gen) return;
268 gen->param = util_Free (gen->param);
269 gen->name = util_Free (gen->name);
270 util_Free (gen);
271 }
272
273
274 /**************************************************************************/
275
276 typedef struct {
277 unif01_Gen *gen0; /* Original generator */
278 long *ILac; /* Table of lacunary indices */
279 int k; /* Size of ILac */
280 int cur; /* Current index in the table */
281 long n;
282 } LacGen_param;
283
LacGen_U01(void * vpar,void * junk)284 static double LacGen_U01 (void *vpar, void *junk)
285 {
286 LacGen_param *paramL = vpar;
287 unif01_Gen *gen = paramL->gen0;
288 int cur = paramL->cur;
289 long *ILac = paramL->ILac;
290 long j;
291 #if 1
292 if (cur > 0) {
293 for (j = 2; j <= ILac[cur] - ILac[cur - 1]; j++)
294 gen->GetU01 (gen->param, gen->state);
295 } else {
296 for (j = 0; j < ILac[0]; j++)
297 gen->GetU01 (gen->param, gen->state);
298 }
299 cur++;
300 if (cur >= paramL->k)
301 cur = 0;
302 paramL->cur = cur;
303
304 #else
305 /* For debugging: write the lacunary indices of the random numbers outputted */
306 if (cur > 0) {
307 for (j = 2; j <= ILac[cur] - ILac[cur - 1]; j++) {
308 gen->GetU01 (gen->param, gen->state);
309 paramL->n++;
310 }
311 } else {
312 for (j = 0; j < ILac[0]; j++) {
313 gen->GetU01 (gen->param, gen->state);
314 paramL->n++;
315 }
316 }
317 cur++;
318 if (cur >= paramL->k)
319 cur = 0;
320 paramL->cur = cur;
321 printf ("Lac = %ld\n", paramL->n);
322 paramL->n++;
323 #endif
324
325 return gen->GetU01 (gen->param, gen->state);
326 }
327
LacGen_Bits(void * vpar,void * junk)328 static unsigned long LacGen_Bits (void *vpar, void *junk)
329 {
330 LacGen_param *paramL = vpar;
331 unif01_Gen *gen = paramL->gen0;
332 int cur = paramL->cur;
333 long *ILac = paramL->ILac;
334 long j;
335
336 if (cur > 0) {
337 for (j = 2; j <= ILac[cur] - ILac[cur - 1]; j++)
338 gen->GetBits (gen->param, gen->state);
339 } else {
340 for (j = 0; j < ILac[0]; j++)
341 gen->GetBits (gen->param, gen->state);
342 }
343 cur++;
344 if (cur >= paramL->k)
345 cur = 0;
346 paramL->cur = cur;
347 return gen->GetBits (gen->param, gen->state);
348 }
349
unif01_CreateLacGen(unif01_Gen * gen,int k,long I[])350 unif01_Gen * unif01_CreateLacGen (unif01_Gen *gen, int k, long I[])
351 {
352 unif01_Gen *genL;
353 LacGen_param *paramL;
354 char name[LEN0 + 1] = "";
355 char str[16];
356 size_t len, len2;
357 int j;
358
359 genL = util_Malloc (sizeof (unif01_Gen));
360 paramL = util_Malloc (sizeof (LacGen_param));
361 paramL->gen0 = gen;
362 paramL->k = k;
363 paramL->cur = 0;
364 paramL->n = 0;
365 paramL->ILac = util_Calloc ((size_t) k, sizeof (long));
366 for (j = 0; j < k; j++)
367 paramL->ILac[j] = I[j];
368
369 len = strlen (gen->name);
370 strncpy (name, gen->name, len);
371 len2 = strlen ("\nunif01_CreateLacGen with k = ");
372 len += len2;
373 strncat (name, "\nunif01_CreateLacGen with k = ", len2);
374 sprintf (str, "%-d", k);
375 strncat (name, str, 16);
376 strncat (name, ", I = (", 8);
377
378 for (j = 0; j < k; j++) {
379 sprintf (str, "%-ld", I[j]);
380 strncat (name, str, 16);
381 if (j < k - 1)
382 strncat (name, ", ", 2);
383 else
384 strncat (name, ")", 1);
385 }
386
387 len = strlen (name);
388 genL->name = util_Calloc (1 + len, sizeof (char));
389 strncpy (genL->name, name, len);
390
391 /* The state of the lacunary generator is simply the state of the original
392 generator */
393 genL->param = paramL;
394 genL->state = gen->state;
395 genL->Write = gen->Write;
396 genL->GetBits = &LacGen_Bits;
397 genL->GetU01 = &LacGen_U01;
398 return genL;
399 }
400
unif01_DeleteLacGen(unif01_Gen * gen)401 void unif01_DeleteLacGen (unif01_Gen *gen)
402 {
403 LacGen_param *param;
404 if (NULL == gen) return;
405 param = gen->param;
406 param->ILac = util_Free (param->ILac);
407 gen->param = util_Free (gen->param);
408 gen->name = util_Free (gen->name);
409 util_Free (gen);
410 }
411
412
413 /**************************************************************************/
414
415 typedef struct {
416 unif01_Gen *gen0; /* The original generator */
417 double R; /* Total probability over [0, a) */
418 double S; /* (R - a) / (1 - a) */
419 double invp; /* Inverse of probability density over [0, a) */
420 double invq; /* Inverse of probability density over [a, 1) */
421 } BiasGen_param;
422
423
BiasGen_U01(void * vpar,void * junk)424 static double BiasGen_U01 (void *vpar, void *junk)
425 {
426 double U;
427 BiasGen_param *paramB = vpar;
428 unif01_Gen *gen = paramB->gen0;
429
430 U = gen->GetU01 (gen->param, gen->state);
431 if (U < paramB->R)
432 return (U * paramB->invp);
433 else
434 return (U - paramB->S) * paramB->invq;
435 }
436
437
BiasGen_Bits(void * vpar,void * junk)438 static unsigned long BiasGen_Bits (void *vpar, void *junk)
439 {
440 return (unsigned long) (unif01_NORM32 * BiasGen_U01 (vpar, junk));
441 }
442
443
unif01_CreateBiasGen(unif01_Gen * gen,double a,double R)444 unif01_Gen * unif01_CreateBiasGen (unif01_Gen *gen, double a, double R)
445 {
446 const double Epsilon = 2.0E-16;
447 unif01_Gen *genB;
448 BiasGen_param *paramB;
449 double p; /* probability density over [0, a) */
450 double q; /* probability density over [a, 1) */
451 char name[LEN0 + 1] = "";
452 char str[16];
453 size_t len;
454
455 util_Assert (R >= 0.0 && R <= 1.0,
456 "unif01_CreateBiasGen: P must be in [0, 1]");
457 util_Assert (a > 0.0 && a < 1.0,
458 "unif01_CreateBiasGen: a must be in (0, 1)");
459
460 genB = util_Malloc (sizeof (unif01_Gen));
461 paramB = util_Malloc (sizeof (BiasGen_param));
462 paramB->gen0 = gen;
463
464 p = R / a;
465 q = (1.0 - R) / (1.0 - a);
466 if (p < Epsilon)
467 paramB->invp = 0.0;
468 else
469 paramB->invp = 1.0 / p;
470 if (q < Epsilon)
471 paramB->invq = 0.0;
472 else
473 paramB->invq = 1.0 / q;
474 paramB->R = R;
475 paramB->S = (p - q) * a;
476
477 strncpy (name, gen->name, LEN0);
478 len = strlen ("\nunif01_CreateBiasGen with P = ");
479 strncat (name, "\nunif01_CreateBiasGen with P = ", len);
480 sprintf (str, "%.4f", R);
481 len = strlen (str);
482 strncat (name, str, len);
483 strncat (name, ", a = ", 8);
484 sprintf (str, "%.4f", a);
485 len = strlen (str);
486 strncat (name, str, len);
487
488 len = strlen (name);
489 genB->name = util_Calloc (1 + len, sizeof (char));
490 strncpy (genB->name, name, len);
491
492 /* The state of the bias generator is simply the state of the original
493 generator */
494 genB->param = paramB;
495 genB->state = gen->state;
496 genB->Write = gen->Write;
497 genB->GetBits = &BiasGen_Bits;
498 genB->GetU01 = &BiasGen_U01;
499 return genB;
500 }
501
unif01_DeleteBiasGen(unif01_Gen * gen)502 void unif01_DeleteBiasGen (unif01_Gen *gen)
503 {
504 if (NULL == gen) return;
505 gen->param = util_Free (gen->param);
506 gen->name = util_Free (gen->name);
507 util_Free (gen);
508 }
509
510
511 /*************************************************************************/
512 typedef struct {
513 unif01_Gen *gen0; /* The original generator */
514 int k; /* keep k numbers */
515 int s; /* skip s numbers */
516 int n; /* state */
517 } LuxGen_param;
518
519
LuxGen_Bits(void * vpar,void * junk)520 static unsigned long LuxGen_Bits (void *vpar, void *junk)
521 {
522 LuxGen_param *paramL = vpar;
523 unif01_Gen *gen = paramL->gen0;
524 if (0 == paramL->n) {
525 int i;
526 for (i = paramL->s; i > 0; --i)
527 gen->GetBits (gen->param, gen->state);
528 paramL->n = paramL->k;
529 }
530 --paramL->n;
531 return gen->GetBits (gen->param, gen->state);
532 }
533
534
LuxGen_U01(void * vpar,void * junk)535 static double LuxGen_U01 (void *vpar, void *junk)
536 {
537 LuxGen_param *paramL = vpar;
538 unif01_Gen *gen = paramL->gen0;
539 if (0 == paramL->n) {
540 int i;
541 for (i = paramL->s; i > 0; --i)
542 gen->GetU01 (gen->param, gen->state);
543 paramL->n = paramL->k;
544 }
545 --paramL->n;
546 return gen->GetU01 (gen->param, gen->state);
547 }
548
549
unif01_CreateLuxGen(unif01_Gen * gen,int k,int L)550 unif01_Gen * unif01_CreateLuxGen (unif01_Gen *gen, int k, int L)
551 {
552 unif01_Gen *genL;
553 LuxGen_param *paramL;
554 char name[LEN0 + 1] = "";
555 char str[26];
556 size_t len;
557 const int s = L - k;
558
559 util_Assert (k > 0, "unif01_CreateLuxGen: k <= 0");
560 util_Assert (k <= L, "unif01_CreateLuxGen: L < k");
561
562 genL = util_Malloc (sizeof (unif01_Gen));
563 paramL = util_Malloc (sizeof (LuxGen_param));
564 paramL->gen0 = gen;
565 paramL->s = s;
566 paramL->k = k;
567 paramL->n = k;
568
569 strncpy (name, gen->name, LEN0);
570 len = strlen ("\nunif01_CreateLuxGen: k = ");
571 strncat (name, "\nunif01_CreateLuxGen: k = ", len);
572 sprintf (str, "%-d, L = %-d", k, L);
573 len = strlen (str);
574 strncat (name, str, len);
575 len = strlen (name);
576 genL->name = util_Calloc (1 + len, sizeof (char));
577 strncpy (genL->name, name, len);
578
579 genL->param = paramL;
580 genL->state = gen->state;
581 genL->Write = gen->Write;
582 genL->GetBits = &LuxGen_Bits;
583 genL->GetU01 = &LuxGen_U01;
584 return genL;
585 }
586
587
unif01_DeleteLuxGen(unif01_Gen * gen)588 void unif01_DeleteLuxGen (unif01_Gen *gen)
589 {
590 if (NULL == gen) return;
591 gen->param = util_Free (gen->param);
592 gen->name = util_Free (gen->name);
593 util_Free (gen);
594 }
595
596
597 /*************************************************************************/
598 typedef struct {
599 unif01_Gen *gen0; /* The original generator */
600 unsigned long mask; /* s most significant bits */
601 } TruncGen_param;
602
603
TruncGen_Bits(void * vpar,void * junk)604 static unsigned long TruncGen_Bits (void *vpar, void *junk)
605 {
606 TruncGen_param *paramT = vpar;
607 unif01_Gen *gen = paramT->gen0;
608
609 return paramT->mask & gen->GetBits (gen->param, gen->state);
610 }
611
TruncGen_U01(void * vpar,void * vsta)612 static double TruncGen_U01 (void *vpar, void *vsta)
613 {
614 return TruncGen_Bits (vpar, vsta) * unif01_INV32;
615 }
616
unif01_CreateTruncGen(unif01_Gen * gen,int b)617 unif01_Gen * unif01_CreateTruncGen (unif01_Gen *gen, int b)
618 {
619 unif01_Gen *genT;
620 TruncGen_param *paramT;
621 char name[LEN0 + 1] = "";
622 char str[16];
623 size_t len;
624
625 if (b < 0)
626 util_Error ("unif01_CreateTruncGen: s < 0");
627 if (b > 32)
628 util_Error ("unif01_CreateTruncGen: s > 32");
629
630 genT = util_Malloc (sizeof (unif01_Gen));
631 paramT = util_Malloc (sizeof (TruncGen_param));
632 paramT->gen0 = gen;
633 if (b >= 32)
634 paramT->mask = 0xffffffffU;
635 else
636 paramT->mask = (0xffffffffU >> (32 - b)) << (32 - b);
637
638 strncpy (name, gen->name, LEN0);
639 len = strlen ("\nunif01_CreateTruncGen with b = ");
640 strncat (name, "\nunif01_CreateTruncGen with b = ", len);
641 sprintf (str, "%-d", b);
642 len = strlen (str);
643 strncat (name, str, len);
644 strncat (name, " bits:", 8);
645
646 len = strlen (name);
647 genT->name = util_Calloc (1 + len, sizeof (char));
648 strncpy (genT->name, name, len);
649
650 /* The state of the trunc generator is simply the state of the original
651 generator */
652 genT->param = paramT;
653 genT->state = gen->state;
654 genT->Write = gen->Write;
655 genT->GetBits = &TruncGen_Bits;
656 genT->GetU01 = &TruncGen_U01;
657 return genT;
658 }
659
unif01_DeleteTruncGen(unif01_Gen * gen)660 void unif01_DeleteTruncGen (unif01_Gen *gen)
661 {
662 if (NULL == gen) return;
663 gen->param = util_Free (gen->param);
664 gen->name = util_Free (gen->name);
665 util_Free (gen);
666 }
667
668
669 /************************************************************************/
670 /*
671 * The original generator is gen0. The r most significant bits of each random
672 * number are dropped, and the s following bits are kept.
673 */
674 typedef struct {
675 unif01_Gen *gen0;
676 int nrows; /* Number of integers used in making a new random number */
677 int B; /* Number of blocks in 1 s-bits group */
678 int w; /* Number of bits in a block */
679 unsigned long maskw; /* Mask of w bits = 2^w - 1 */
680 int r;
681 int s;
682 } BitBlock_param;
683
684
685 typedef struct {
686 unsigned long *Z;
687 int n; /* Build n random numbers at a time, n <= 32 */
688 BitBlock_param *param;
689 } BitBlock_state;
690
691
BitBlock_Bits(void * vpar,void * vsta)692 static unsigned long BitBlock_Bits (void *vpar, void *vsta)
693 {
694 BitBlock_state *state = vsta;
695
696 if (state->n <= 0) {
697 BitBlock_param *param = vpar;
698 unsigned long X;
699 int i, j;
700
701 /* Generate B random integers Z from the bits of nrows random integers
702 X from the original gen0 */
703 for (j = 0; j < param->B; j++)
704 state->Z[j] = 0;
705 for (i = 0; i < param->nrows; i++) {
706 /* Get a random integer X of s bits */
707 X = unif01_StripB (param->gen0, param->r, param->s);
708 /* Take w of the s bits of X to make bits of Z[j] */
709 for (j = 0; j < param->B; j++) {
710 state->Z[j] <<= param->w;
711 state->Z[j] |= X & param->maskw;
712 X >>= param->w;
713 }
714 }
715 state->n = param->B;
716 }
717 return state->Z[--state->n];
718 }
719
720
BitBlock_U01(void * vpar,void * vsta)721 static double BitBlock_U01 (void *vpar, void *vsta)
722 {
723 return BitBlock_Bits (vpar, vsta) * unif01_INV32;
724 }
725
726
WrBitBlock(void * vsta)727 static void WrBitBlock (void *vsta)
728 {
729 BitBlock_state *state = vsta;
730 state->param->gen0->Write (state->param->gen0->state);
731 }
732
733
unif01_CreateBitBlockGen(unif01_Gen * gen,int r,int s,int w)734 unif01_Gen * unif01_CreateBitBlockGen (unif01_Gen *gen, int r, int s, int w)
735 {
736 unif01_Gen *genV;
737 BitBlock_param *paramV;
738 BitBlock_state *stateV;
739 char *name;
740 char str[64];
741 size_t len, len2, len3;
742
743 util_Assert (s > 0, "unif01_CreateBitBlockGen: s <= 0");
744 util_Assert (r >= 0, "unif01_CreateBitBlockGen: r < 0");
745 util_Assert (r + s <= 32, "unif01_CreateBitBlockGen: r + s > 32");
746 util_Assert (w > 0, "unif01_CreateBitBlockGen: w < 1");
747 util_Assert (32 % w == 0, "unif01_CreateBitBlockGen: w must divide 32");
748
749 genV = util_Malloc (sizeof (unif01_Gen));
750 paramV = util_Malloc (sizeof (BitBlock_param));
751 stateV = util_Malloc (sizeof (BitBlock_state));
752 paramV->gen0 = gen;
753 paramV->s = s;
754 paramV->r = r;
755 paramV->w = w;
756 paramV->B = s / w;
757 paramV->maskw = num_TwoExp[paramV->w] - 1.0;
758 paramV->nrows = 32 / w;
759 stateV->param = paramV;
760 stateV->n = 0;
761 stateV->Z = util_Calloc ((size_t) paramV->B, sizeof (unsigned long));
762
763 len = strlen (gen->name);
764 len2 = strlen ("\nunif01_CreateBitBlockGen: ");
765 len += len2;
766 sprintf (str, "r = %1d, s = %1d, w = %1d", r, s, w);
767 len3 = strlen (str);
768 len += len3;
769 name = util_Calloc (len + 1, sizeof (char));
770 strncpy (name, gen->name, len);
771 strncat (name, "\nunif01_CreateBitBlockGen: ", len2);
772 strncat (name, str, len3);
773 genV->name = name;
774 genV->param = paramV;
775 genV->state = stateV;
776 genV->Write = &WrBitBlock;
777 genV->GetBits = &BitBlock_Bits;
778 genV->GetU01 = &BitBlock_U01;
779 return genV;
780 }
781
782
unif01_DeleteBitBlockGen(unif01_Gen * gen)783 void unif01_DeleteBitBlockGen (unif01_Gen *gen)
784 {
785 BitBlock_state *state;
786 if (NULL == gen) return;
787 state = gen->state;
788 state->Z = util_Free (state->Z);
789 gen->param = util_Free (gen->param);
790 gen->state = util_Free (gen->state);
791 gen->name = util_Free (gen->name);
792 util_Free (gen);
793 }
794
795
796 /************************************************************************/
797
CombGen2_U01_Add(void * vpar,void * junk)798 static double CombGen2_U01_Add (void *vpar, void *junk)
799 {
800 unif01_Comb2_Param *g = vpar;
801 unif01_Gen *gen1 = g->gen1;
802 unif01_Gen *gen2 = g->gen2;
803 double U;
804
805 U = gen1->GetU01 (gen1->param, gen1->state) +
806 gen2->GetU01 (gen2->param, gen2->state);
807 if (U >= 1.0)
808 return (U - 1.0);
809 else
810 return U;
811 }
812
813
CombGen2_Bits_Add(void * vpar,void * junk)814 static unsigned long CombGen2_Bits_Add (void *vpar, void *junk)
815 {
816 return (unsigned long) (unif01_NORM32 * CombGen2_U01_Add (vpar, junk));
817 }
818
819
CombGen2_Bits_Xor(void * vpar,void * junk)820 static unsigned long CombGen2_Bits_Xor (void *vpar, void *junk)
821 {
822 unif01_Comb2_Param *g = vpar;
823 unif01_Gen *gen1 = g->gen1;
824 unif01_Gen *gen2 = g->gen2;
825
826 return gen1->GetBits (gen1->param, gen1->state) ^
827 gen2->GetBits (gen2->param, gen2->state);
828 }
829
830
CombGen2_U01_Xor(void * vpar,void * junk)831 static double CombGen2_U01_Xor (void *vpar, void *junk)
832 {
833 return CombGen2_Bits_Xor (vpar, junk) * unif01_INV32;
834 }
835
836
WrCombGen2(void * vsta)837 static void WrCombGen2 (void *vsta)
838 {
839 unif01_Comb2_Param *g = vsta;
840 printf ("2 Combined Generators:\n");
841 g->gen1->Write (g->gen1->state);
842 g->gen2->Write (g->gen2->state);
843 }
844
845
CreateCombGen2(unif01_Gen * g1,unif01_Gen * g2,char * mess,char * name)846 static unif01_Gen * CreateCombGen2 (unif01_Gen *g1, unif01_Gen *g2,
847 char *mess, char *name)
848 {
849 unif01_Gen *gen;
850 unif01_Comb2_Param *paramC;
851 size_t len, L;
852
853 gen = util_Malloc (sizeof (unif01_Gen));
854 paramC = util_Malloc (sizeof (unif01_Comb2_Param));
855 paramC->gen1 = g1;
856 paramC->gen2 = g2;
857
858 len = strlen (g1->name) + strlen (g2->name) + strlen (name) + strlen (mess);
859 len += 5;
860 gen->name = util_Calloc (len + 1, sizeof (char));
861 L = strlen (mess);
862 if (L > 0) {
863 strncpy (gen->name, mess, len);
864 if (mess[L - 1] != ':')
865 strncat (gen->name, ":", 3);
866 strncat (gen->name, "\n", 3);
867 }
868 strncat (gen->name, g1->name, len);
869 strncat (gen->name, "\n", 3);
870 strncat (gen->name, g2->name, len);
871 strncat (gen->name, name, len);
872
873 gen->param = paramC;
874 gen->state = paramC;
875 gen->Write = &WrCombGen2;
876 return gen;
877 }
878
879
unif01_CreateCombAdd2(unif01_Gen * g1,unif01_Gen * g2,char * Mess)880 unif01_Gen * unif01_CreateCombAdd2 (unif01_Gen *g1, unif01_Gen *g2, char *Mess)
881 {
882 unif01_Gen *gen;
883 gen = CreateCombGen2 (g1, g2, Mess, "\nunif01_CreateCombAdd2");
884 gen->GetU01 = &CombGen2_U01_Add;
885 gen->GetBits = &CombGen2_Bits_Add;
886 return gen;
887 }
888
889
unif01_CreateCombXor2(unif01_Gen * g1,unif01_Gen * g2,char * Mess)890 unif01_Gen * unif01_CreateCombXor2 (unif01_Gen *g1, unif01_Gen *g2,
891 char *Mess)
892 {
893 unif01_Gen *gen;
894 gen = CreateCombGen2 (g1, g2, Mess, "\nunif01_CreateCombXor2");
895 gen->GetU01 = &CombGen2_U01_Xor;
896 gen->GetBits = &CombGen2_Bits_Xor;
897 return gen;
898 }
899
unif01_DeleteCombGen(unif01_Gen * gen)900 void unif01_DeleteCombGen (unif01_Gen *gen)
901 {
902 if (NULL == gen) return;
903 gen->param = util_Free (gen->param);
904 gen->name = util_Free (gen->name);
905 util_Free (gen);
906 }
907
908
909 /************************************************************************/
910
911 typedef struct {
912 unif01_Gen *gen1;
913 unif01_Gen *gen2;
914 unif01_Gen *gen3;
915 } Comb3_Param;
916
917
CombGen3_U01_Add(void * vpar,void * junk)918 static double CombGen3_U01_Add (void *vpar, void *junk)
919 {
920 Comb3_Param *g = vpar;
921 unif01_Gen *gen1 = g->gen1;
922 unif01_Gen *gen2 = g->gen2;
923 unif01_Gen *gen3 = g->gen3;
924 double U;
925
926 /*
927 When the combined generator is used to generate random integers, in rare
928 cases, an integer may differ by 1 unit depending on the order of
929 addition of the 3 terms (one from each component). This is due
930 to the last bit (bit 53) of the value returned which may be affected by
931 floating-point numerical errors. Furthermore, the result
932 may be different if the addition is done without function calls
933 (as in the pre-programmed version of Wichmann-Hill for example in
934 {\tt ulcg\_CreateCombWH3}), in which case, the 2 extra guard bits
935 required by the IEEE-754 standard in floating-point arithmetic
936 operations may give a more precise result.
937 */
938 U = gen1->GetU01 (gen1->param, gen1->state) +
939 gen2->GetU01 (gen2->param, gen2->state) +
940 gen3->GetU01 (gen3->param, gen3->state);
941
942 if (U < 1.0)
943 return U;
944 if (U < 2.0)
945 return (U - 1.0);
946 return U - 2.0;
947 }
948
949
CombGen3_Bits_Add(void * vpar,void * junk)950 static unsigned long CombGen3_Bits_Add (void *vpar, void *junk)
951 {
952 return (unsigned long) (CombGen3_U01_Add (vpar, junk) * unif01_NORM32);
953 }
954
955
CombGen3_Bits_Xor(void * vpar,void * junk)956 static unsigned long CombGen3_Bits_Xor (void *vpar, void *junk)
957 {
958 Comb3_Param *g = vpar;
959 unif01_Gen *gen1 = g->gen1;
960 unif01_Gen *gen2 = g->gen2;
961 unif01_Gen *gen3 = g->gen3;
962
963 return gen1->GetBits (gen1->param, gen1->state) ^
964 gen2->GetBits (gen2->param, gen2->state) ^
965 gen3->GetBits (gen3->param, gen3->state);
966 }
967
968
CombGen3_U01_Xor(void * vpar,void * junk)969 static double CombGen3_U01_Xor (void *vpar, void *junk)
970 {
971 return CombGen3_Bits_Xor (vpar, junk) * unif01_INV32;
972 }
973
974
WrCombGen3(void * vsta)975 static void WrCombGen3 (void *vsta )
976 {
977 Comb3_Param *g = vsta;
978 printf ("3 Combined Generators:\n");
979 g->gen1->Write (g->gen1->state);
980 g->gen2->Write (g->gen2->state);
981 g->gen3->Write (g->gen3->state);
982 }
983
984
CreateCombGen3(unif01_Gen * g1,unif01_Gen * g2,unif01_Gen * g3,const char * mess,const char * name)985 static unif01_Gen * CreateCombGen3 (unif01_Gen *g1, unif01_Gen *g2,
986 unif01_Gen *g3, const char *mess, const char *name)
987 {
988 unif01_Gen *gen;
989 Comb3_Param *paramC;
990 size_t len, L;
991
992 gen = util_Malloc (sizeof (unif01_Gen));
993 paramC = util_Malloc (sizeof (Comb3_Param));
994 paramC->gen1 = g1;
995 paramC->gen2 = g2;
996 paramC->gen3 = g3;
997
998 len = strlen (g1->name) + strlen (g2->name) + strlen (g3->name) +
999 strlen (name) + strlen (mess);
1000 len += 5;
1001 gen->name = util_Calloc (len + 1, sizeof (char));
1002 L = strlen (mess);
1003 if (L > 0) {
1004 strncpy (gen->name, mess, len);
1005 if (mess[L - 1] != ':')
1006 strncat (gen->name, ":", 3);
1007 strncat (gen->name, "\n", 3);
1008 }
1009 strncat (gen->name, g1->name, len);
1010 strncat (gen->name, "\n", 3);
1011 strncat (gen->name, g2->name, len);
1012 strncat (gen->name, "\n", 3);
1013 strncat (gen->name, g3->name, len);
1014 strncat (gen->name, name, len);
1015
1016 gen->param = paramC;
1017 gen->state = paramC;
1018 gen->Write = &WrCombGen3;
1019 return gen;
1020 }
1021
1022
unif01_CreateCombAdd3(unif01_Gen * g1,unif01_Gen * g2,unif01_Gen * g3,char * mess)1023 unif01_Gen * unif01_CreateCombAdd3 (unif01_Gen *g1, unif01_Gen *g2,
1024 unif01_Gen *g3, char *mess)
1025 {
1026 unif01_Gen *gen;
1027 gen = CreateCombGen3 (g1, g2, g3, mess, "\nunif01_CreateCombAdd3");
1028 gen->GetU01 = &CombGen3_U01_Add;
1029 gen->GetBits = &CombGen3_Bits_Add;
1030 return gen;
1031 }
1032
1033
unif01_CreateCombXor3(unif01_Gen * g1,unif01_Gen * g2,unif01_Gen * g3,char * mess)1034 unif01_Gen * unif01_CreateCombXor3 (unif01_Gen *g1, unif01_Gen *g2,
1035 unif01_Gen *g3, char *mess)
1036 {
1037 unif01_Gen *gen;
1038 gen = CreateCombGen3 (g1, g2, g3, mess, "\nunif01_CreateCombXor3");
1039 gen->GetU01 = &CombGen3_U01_Xor;
1040 gen->GetBits = &CombGen3_Bits_Xor;
1041 return gen;
1042 }
1043
1044
1045 /*=========================================================================*/
1046
1047 typedef struct {
1048 int j; /* Which random number */
1049 int i; /* Which generator */
1050 int L;
1051 int k; /* Number of parallel generators */
1052 unif01_Gen **agen; /* Parallel generators */
1053 } ParallelGen_state;
1054
1055
ParallelGen_U01(void * junk,void * vsta)1056 static double ParallelGen_U01 (void *junk, void *vsta)
1057 {
1058 ParallelGen_state *stateP = vsta;
1059 unif01_Gen *g;
1060
1061 if (++stateP->j >= stateP->L) {
1062 stateP->j = 0;
1063 if (++stateP->i >= stateP->k)
1064 stateP->i = 0;
1065 }
1066 g = stateP->agen[stateP->i];
1067 return g->GetU01 (g->param, g->state);
1068 }
1069
1070
ParallelGen_Bits(void * junk,void * vsta)1071 static unsigned long ParallelGen_Bits (void *junk, void *vsta)
1072 {
1073 ParallelGen_state *stateP = vsta;
1074 unif01_Gen *g;
1075
1076 if (++stateP->j >= stateP->L) {
1077 stateP->j = 0;
1078 if (++stateP->i >= stateP->k)
1079 stateP->i = 0;
1080 }
1081 g = stateP->agen[stateP->i];
1082 return g->GetBits (g->param, g->state);
1083 }
1084
1085
WrParallelGen(void * vsta)1086 static void WrParallelGen (void *vsta)
1087 {
1088 int i;
1089 ParallelGen_state *state = vsta;
1090 printf (" i = %d, j = %d\n\nParallel Generators:\n", state->i, state->j);
1091 for (i = 0; i < state->k; ++i)
1092 unif01_WriteNameGen(state->agen[i]);
1093 }
1094
1095
unif01_CreateParallelGen(int k,unif01_Gen * gen[],int L)1096 unif01_Gen * unif01_CreateParallelGen (int k, unif01_Gen *gen[], int L)
1097 {
1098 #define NCAT 16
1099 unif01_Gen *genP;
1100 ParallelGen_state *stateP;
1101 char name[LEN0 + 1] = {0};
1102 char str[NCAT + 1];
1103 size_t len;
1104 int j;
1105
1106 genP = util_Malloc (sizeof (unif01_Gen));
1107 stateP = util_Malloc (sizeof (ParallelGen_state));
1108 stateP->k = k;
1109 stateP->L = L;
1110 stateP->i = k;
1111 stateP->j = L;
1112 stateP->agen = util_Calloc ((size_t) k, sizeof (unif01_Gen *));
1113 for (j = 0; j < k; j++)
1114 stateP->agen[j] = gen[j];
1115
1116 len = strlen ("unif01_CreateParallelGen: k = ");
1117 strncpy (name, "unif01_CreateParallelGen: k = ", len);
1118 sprintf (str, "%-d", k);
1119 strncat (name, str, NCAT);
1120 strncat (name, ", L = ", NCAT);
1121 sprintf (str, "%-d", L);
1122 strncat (name, str, NCAT);
1123 len = strlen (name);
1124 genP->name = util_Calloc (1 + len, sizeof (char));
1125 strncpy (genP->name, name, len);
1126
1127 genP->state = stateP;
1128 genP->Write = &WrParallelGen;
1129 genP->GetBits = &ParallelGen_Bits;
1130 genP->GetU01 = &ParallelGen_U01;
1131 return genP;
1132 #undef NCAT
1133 }
1134
1135
unif01_DeleteParallelGen(unif01_Gen * gen)1136 void unif01_DeleteParallelGen (unif01_Gen *gen)
1137 {
1138 ParallelGen_state *state;
1139 if (NULL == gen) return;
1140 state = gen->state;
1141 state->agen = util_Free (state->agen);
1142 gen->state = util_Free (gen->state);
1143 gen->name = util_Free (gen->name);
1144 util_Free (gen);
1145 }
1146
1147
1148 /*=========================================================================*/
1149
1150 static double (*externGen_U01)(void); /* The external generator U01 */
1151 static int coGU = 0; /* Counter for GU_U01 */
1152
1153
GU_U01(void * junk1,void * junk2)1154 static double GU_U01 (void *junk1, void *junk2)
1155 {
1156 return externGen_U01 ();
1157 }
1158
1159
GU_Bits(void * junk1,void * junk2)1160 static unsigned long GU_Bits (void *junk1, void *junk2)
1161 {
1162 return (unsigned long) (externGen_U01 () * unif01_NORM32);
1163 }
1164
1165
WrExternGen(void * junk2)1166 static void WrExternGen (void *junk2)
1167 {
1168 }
1169
1170
unif01_CreateExternGen01(char * name,double (* f_U01)(void))1171 unif01_Gen *unif01_CreateExternGen01 (char *name, double (*f_U01)(void))
1172 {
1173 unif01_Gen *gen;
1174 size_t leng;
1175
1176 util_Assert (coGU == 0,
1177 "unif01_CreateExternGen01: only 1 such generator can be in use");
1178 coGU++;
1179 gen = util_Malloc (sizeof (unif01_Gen));
1180 gen->state = NULL;
1181 gen->param = NULL;
1182 gen->Write = WrExternGen;
1183 externGen_U01 = f_U01;
1184 gen->GetU01 = GU_U01;
1185 gen->GetBits = GU_Bits;
1186
1187 if (name) {
1188 leng = strlen (name);
1189 gen->name = util_Calloc (leng + 2, sizeof (char));
1190 strncpy (gen->name, name, leng);
1191 } else {
1192 gen->name = util_Calloc (1, sizeof (char));
1193 gen->name[0] = '\0';
1194 }
1195 return gen;
1196 }
1197
1198
unif01_DeleteExternGen01(unif01_Gen * gen)1199 void unif01_DeleteExternGen01 (unif01_Gen * gen)
1200 {
1201 if (NULL == gen)
1202 return;
1203 gen->name = util_Free (gen->name);
1204 util_Free (gen);
1205 coGU--;
1206 }
1207
1208
1209 /*=========================================================================*/
1210
1211 static unsigned int (*externGen_Bits)(void);
1212 static int coGB = 0; /* Counter for GB_U01 */
1213
GB_U01(void * junk1,void * junk2)1214 static double GB_U01 (void *junk1, void *junk2)
1215 {
1216 return externGen_Bits () / unif01_NORM32;
1217 }
1218
1219
GB_Bits(void * junk1,void * junk2)1220 static unsigned long GB_Bits (void *junk1, void *junk2)
1221 {
1222 return externGen_Bits ();
1223 }
1224
1225
unif01_CreateExternGenBits(char * name,unsigned int (* f_Bits)(void))1226 unif01_Gen *unif01_CreateExternGenBits (char *name,
1227 unsigned int (*f_Bits)(void))
1228 {
1229 unif01_Gen *gen;
1230 size_t leng;
1231
1232 util_Assert (coGB == 0,
1233 "unif01_CreateExternGenBits: only 1 such generator can be in use");
1234 coGB++;
1235 gen = util_Malloc (sizeof (unif01_Gen));
1236 gen->state = NULL;
1237 gen->param = NULL;
1238 gen->Write = WrExternGen;
1239 externGen_Bits = f_Bits;
1240 gen->GetU01 = GB_U01;
1241 gen->GetBits = GB_Bits;
1242
1243 if (name) {
1244 leng = strlen (name);
1245 gen->name = util_Calloc (leng + 2, sizeof (char));
1246 strncpy (gen->name, name, leng);
1247 } else {
1248 gen->name = util_Calloc (1, sizeof (char));
1249 gen->name[0] = '\0';
1250 }
1251 return gen;
1252 }
1253
1254
unif01_DeleteExternGenBits(unif01_Gen * gen)1255 void unif01_DeleteExternGenBits (unif01_Gen * gen)
1256 {
1257 if (NULL == gen)
1258 return;
1259 gen->name = util_Free (gen->name);
1260 util_Free (gen);
1261 coGB--;
1262 }
1263
1264
1265 /*=========================================================================*/
1266
1267 static unsigned long (*externGenLong_Bits)(void);
1268 static int coGBL = 0; /* Counter for GBLong_U01 */
1269
GBLong_U01(void * junk1,void * junk2)1270 static double GBLong_U01 (void *junk1, void *junk2)
1271 {
1272 return externGenLong_Bits () / unif01_NORM32;
1273 }
1274
1275
GBLong_Bits(void * junk1,void * junk2)1276 static unsigned long GBLong_Bits (void *junk1, void *junk2)
1277 {
1278 return externGenLong_Bits ();
1279 }
1280
1281
unif01_CreateExternGenBitsL(char * name,unsigned long (* f_Bits)(void))1282 unif01_Gen *unif01_CreateExternGenBitsL (char *name,
1283 unsigned long (*f_Bits)(void))
1284 {
1285 unif01_Gen *gen;
1286 size_t leng;
1287
1288 util_Assert (coGBL == 0,
1289 "unif01_CreateExternGenBitsL: only 1 such generator can be in use");
1290 coGBL++;
1291 gen = util_Malloc (sizeof (unif01_Gen));
1292 gen->state = NULL;
1293 gen->param = NULL;
1294 gen->Write = WrExternGen;
1295 externGenLong_Bits = f_Bits;
1296 gen->GetU01 = GBLong_U01;
1297 gen->GetBits = GBLong_Bits;
1298
1299 if (name) {
1300 leng = strlen (name);
1301 gen->name = util_Calloc (leng + 2, sizeof (char));
1302 strncpy (gen->name, name, leng);
1303 } else {
1304 gen->name = util_Calloc (1, sizeof (char));
1305 gen->name[0] = '\0';
1306 }
1307 return gen;
1308 }
1309
1310
unif01_DeleteExternGenBitsL(unif01_Gen * gen)1311 void unif01_DeleteExternGenBitsL (unif01_Gen * gen)
1312 {
1313 if (NULL == gen)
1314 return;
1315 gen->name = util_Free (gen->name);
1316 util_Free (gen);
1317 coGBL--;
1318 }
1319
1320
1321 /**************************************************************************/
1322
unif01_TimerGen(unif01_Gen * gen,unif01_TimerRec * pt,long n,lebool fU01)1323 void unif01_TimerGen (unif01_Gen *gen, unif01_TimerRec * pt, long n,
1324 lebool fU01)
1325 {
1326 chrono_Chrono *C1;
1327 double U;
1328 unsigned long V;
1329 long i;
1330
1331 C1 = chrono_Create ();
1332 if (fU01)
1333 for (i = 0; i < n; i++)
1334 U = gen->GetU01 (gen->param, gen->state);
1335 else
1336 for (i = 0; i < n; i++)
1337 V = gen->GetBits (gen->param, gen->state);
1338 pt->time = chrono_Val (C1, chrono_sec);
1339 pt->mean = 0.0;
1340 pt->n = n;
1341 pt->fU01 = fU01;
1342 pt->gen = gen;
1343 chrono_Delete (C1);
1344 }
1345
unif01_TimerSumGen(unif01_Gen * gen,unif01_TimerRec * pt,long n,lebool fU01)1346 void unif01_TimerSumGen (unif01_Gen *gen, unif01_TimerRec * pt, long n,
1347 lebool fU01)
1348 {
1349 chrono_Chrono *C1;
1350 double Sum = 0.0;
1351 unsigned long Y = 0;
1352 long i;
1353
1354 C1 = chrono_Create ();
1355 if (fU01)
1356 for (i = 0; i < n; i++)
1357 Sum += gen->GetU01 (gen->param, gen->state);
1358 else
1359 for (i = 0; i < n; i++)
1360 Y += gen->GetBits (gen->param, gen->state);
1361 pt->time = chrono_Val (C1, chrono_sec);
1362 if (fU01)
1363 pt->mean = Sum / n;
1364 else
1365 pt->mean = (double) Y / n;
1366 pt->n = n;
1367 pt->gen = gen;
1368 pt->fU01 = fU01;
1369 chrono_Delete (C1);
1370 }
1371
unif01_WriteTimerRec(unif01_TimerRec * R)1372 void unif01_WriteTimerRec (unif01_TimerRec *R)
1373 {
1374 unif01_Gen *gen = R->gen;
1375 char stri [LEN1 + 1] = "";
1376 char *p;
1377 size_t len;
1378
1379 printf ("\n------------- Results of speed test ---------------");
1380 printf ("\n\n Host: ");
1381 if (swrite_Host)
1382 gdef_WriteHostName ();
1383 else
1384 printf ("\n");
1385
1386 /* Print only the generator name, without the parameters or seeds. */
1387 /* The parameters start after the first blank; name ends with ':' */
1388 printf (" Generator: ");
1389 len = strcspn (gen->name, ":");
1390 strncpy (stri, gen->name, len);
1391 stri [len] = '\0';
1392 printf ("%s", stri);
1393 p = strstr (gen->name, "unif01");
1394 while (p != NULL) {
1395 /* For Filters or Combined generators */
1396 len = strcspn (p, " \0");
1397 strncpy (stri, p, len);
1398 stri [len] = '\0';
1399 printf (", %s", stri);
1400 p += len;
1401 p = strstr (p, "unif01");
1402 }
1403 if (R->fU01) {
1404 printf ("\n Method: GetU01");
1405 if (R->mean > 0.0)
1406 printf ("\n Mean = %.15f", R->mean);
1407 } else {
1408 printf ("\n Method: GetBits");
1409 if (R->mean > 0.0)
1410 printf ("\n Mean = %.16g", R->mean);
1411 }
1412 printf ("\n Number of calls: %ld", R->n);
1413 printf ("\n Total CPU time: ");
1414 printf ("%6.2f sec\n\n", R->time);
1415 }
1416
unif01_TimerGenWr(unif01_Gen * gen,long n,lebool fU01)1417 void unif01_TimerGenWr (unif01_Gen *gen, long n, lebool fU01)
1418 {
1419 unif01_TimerRec timer;
1420 unif01_TimerGen (gen, &timer, n, fU01);
1421 unif01_WriteTimerRec (&timer);
1422 }
1423
unif01_TimerSumGenWr(unif01_Gen * gen,long n,lebool fU01)1424 void unif01_TimerSumGenWr (unif01_Gen *gen, long n, lebool fU01)
1425 {
1426 unif01_TimerRec timer;
1427 unif01_TimerSumGen (gen, &timer, n, fU01);
1428 unif01_WriteTimerRec (&timer);
1429 }
1430