1 # include <cstdlib>
2 # include <iostream>
3 # include <iomanip>
4 # include <ctime>
5
6 using namespace std;
7
8 # include "rnglib.hpp"
9
10 namespace dream {
11
12 //****************************************************************************80
13
advance_state(int k)14 void advance_state ( int k )
15
16 //****************************************************************************80
17 //
18 // Purpose:
19 //
20 // ADVANCE_STATE advances the state of the current generator.
21 //
22 // Discussion:
23 //
24 // This procedure advances the state of the current generator by 2^K
25 // values and resets the initial seed to that value.
26 //
27 // Licensing:
28 //
29 // This code is distributed under the GNU LGPL license.
30 //
31 // Modified:
32 //
33 // 30 March 2013
34 //
35 // Author:
36 //
37 // Original Pascal version by Pierre L'Ecuyer, Serge Cote.
38 // C++ version by John Burkardt.
39 //
40 // Reference:
41 //
42 // Pierre LEcuyer, Serge Cote,
43 // Implementing a Random Number Package with Splitting Facilities,
44 // ACM Transactions on Mathematical Software,
45 // Volume 17, Number 1, March 1991, pages 98-111.
46 //
47 // Parameters:
48 //
49 // Input, int K, indicates that the generator is to be
50 // advanced by 2^K values.
51 // 0 <= K.
52 //
53 {
54 const int a1 = 40014;
55 const int a2 = 40692;
56 int b1;
57 int b2;
58 int cg1;
59 int cg2;
60 int g;
61 int i;
62 const int m1 = 2147483563;
63 const int m2 = 2147483399;
64
65 if ( k < 0 )
66 {
67 cerr << "\n";
68 cerr << "ADVANCE_STATE - Fatal error!\n";
69 cerr << " Input exponent K is out of bounds.\n";
70 exit ( 1 );
71 }
72 //
73 // Check whether the package must be initialized.
74 //
75 if ( ! initialized_get ( ) )
76 {
77 cout << "\n";
78 cout << "ADVANCE_STATE - Note:\n";
79 cout << " Initializing RNGLIB package.\n";
80 initialize ( );
81 }
82 //
83 // Get the current generator index.
84 //
85 g = cgn_get ( );
86
87 b1 = a1;
88 b2 = a2;
89
90 for ( i = 1; i <= k; k++ )
91 {
92 b1 = multmod ( b1, b1, m1 );
93 b2 = multmod ( b2, b2, m2 );
94 }
95
96 cg_get ( g, cg1, cg2 );
97 cg1 = multmod ( b1, cg1, m1 );
98 cg2 = multmod ( b2, cg2, m2 );
99 cg_set ( g, cg1, cg2 );
100
101 return;
102 }
103 //****************************************************************************80
104
antithetic_get()105 bool antithetic_get ( )
106
107 //***************************************************************************80
108 //
109 // Purpose:
110 //
111 // ANTITHETIC_GET queries the antithetic value for a given generator.
112 //
113 // Licensing:
114 //
115 // This code is distributed under the GNU LGPL license.
116 //
117 // Modified:
118 //
119 // 02 April 2013
120 //
121 // Author:
122 //
123 // John Burkardt
124 //
125 // Parameters:
126 //
127 // Output, bool ANTITHETIC_GET, is TRUE if generator G is antithetic.
128 //
129 {
130 int i;
131 bool value;
132
133 i = -1;
134 antithetic_memory ( i, value );
135
136 return value;
137 }
138 //****************************************************************************80
139
antithetic_memory(int i,bool & value)140 void antithetic_memory ( int i, bool &value )
141
142 //****************************************************************************80
143 //
144 // Purpose:
145 //
146 // ANTITHETIC_MEMORY stores the antithetic value for all generators.
147 //
148 // Licensing:
149 //
150 // This code is distributed under the GNU LGPL license.
151 //
152 // Modified:
153 //
154 // 02 April 2013
155 //
156 // Author:
157 //
158 // John Burkardt
159 //
160 // Parameters:
161 //
162 // Input, int I, the desired action.
163 // -1, get a value.
164 // 0, initialize all values.
165 // 1, set a value.
166 //
167 // Input/output, bool &VALUE. For I = -1, VALUE is an output
168 // quantity, if I = +1, then VALUE is an input quantity.
169 //
170 {
171 # define G_MAX 32
172
173 static bool a_save[G_MAX];
174 int g;
175 const int g_max = 32;
176 int j;
177
178 if ( i < 0 )
179 {
180 g = cgn_get ( );
181 value = a_save[g];
182 }
183 else if ( i == 0 )
184 {
185 for ( j = 0; j < g_max; j++ )
186 {
187 a_save[j] = false;
188 }
189 }
190 else if ( 0 < i )
191 {
192 g = cgn_get ( );
193 a_save[g] = value;
194 }
195
196 return;
197 # undef G_MAX
198 }
199 //****************************************************************************80
200
antithetic_set(bool value)201 void antithetic_set ( bool value )
202
203 //****************************************************************************80
204 //
205 // Purpose:
206 //
207 // ANTITHETIC_SET sets the antithetic value for a given generator.
208 //
209 // Licensing:
210 //
211 // This code is distributed under the GNU LGPL license.
212 //
213 // Modified:
214 //
215 // 02 April 2013
216 //
217 // Author:
218 //
219 // John Burkardt
220 //
221 // Parameters:
222 //
223 // Input, bool VALUE, is TRUE if generator G is to be antithetic.
224 //
225 {
226 int i;
227
228 i = +1;
229 antithetic_memory ( i, value );
230
231 return;
232 }
233 //****************************************************************************80
234
cg_get(int g,int & cg1,int & cg2)235 void cg_get ( int g, int &cg1, int &cg2 )
236
237 //****************************************************************************80
238 //
239 // Purpose:
240 //
241 // CG_GET queries the CG values for a given generator.
242 //
243 // Licensing:
244 //
245 // This code is distributed under the GNU LGPL license.
246 //
247 // Modified:
248 //
249 // 27 March 2013
250 //
251 // Author:
252 //
253 // John Burkardt
254 //
255 // Parameters:
256 //
257 // Input, int G, the index of the generator.
258 // 0 <= G <= 31.
259 //
260 // Output, int &CG1, &CG2, the CG values for generator G.
261 //
262 {
263 int i;
264
265 i = -1;
266 cg_memory ( i, g, cg1, cg2 );
267
268 return;
269 }
270 //****************************************************************************80
271
cg_memory(int i,int g,int & cg1,int & cg2)272 void cg_memory ( int i, int g, int &cg1, int &cg2 )
273
274 //****************************************************************************80
275 //
276 // Purpose:
277 //
278 // CG_MEMORY stores the CG values for all generators.
279 //
280 // Licensing:
281 //
282 // This code is distributed under the GNU LGPL license.
283 //
284 // Modified:
285 //
286 // 30 March 2013
287 //
288 // Author:
289 //
290 // John Burkardt
291 //
292 // Parameters:
293 //
294 // Input, int I, the desired action.
295 // -1, get a value.
296 // 0, initialize all values.
297 // 1, set a value.
298 //
299 // Input, int G, for I = -1 or +1, the index of
300 // the generator, with 0 <= G <= 31.
301 //
302 // Input/output, int &CG1, &CG2. For I = -1,
303 // these are output, for I = +1, these are input, for I = 0,
304 // these arguments are ignored. When used, the arguments are
305 // old or new values of the CG parameter for generator G.
306 //
307 {
308 # define G_MAX 32
309
310 static int cg1_save[G_MAX];
311 static int cg2_save[G_MAX];
312 const int g_max = 32;
313 int j;
314
315 if ( g < 0 || g_max <= g )
316 {
317 cerr << "\n";
318 cerr << "CG_MEMORY - Fatal error!\n";
319 cerr << " Input generator index G is out of bounds.\n";
320 exit ( 1 );
321 }
322
323 if ( i < 0 )
324 {
325 cg1 = cg1_save[g];
326 cg2 = cg2_save[g];
327 }
328 else if ( i == 0 )
329 {
330 for ( j = 0; j < g_max; j++ )
331 {
332 cg1_save[j] = 0;
333 cg2_save[j] = 0;
334 }
335 }
336 else if ( 0 < i )
337 {
338 cg1_save[g] = cg1;
339 cg2_save[g] = cg2;
340 }
341
342 return;
343 # undef G_MAX
344 }
345 //****************************************************************************80
346
cg_set(int g,int cg1,int cg2)347 void cg_set ( int g, int cg1, int cg2 )
348
349 //****************************************************************************80
350 //
351 // Purpose:
352 //
353 // CG_SET sets the CG values for a given generator.
354 //
355 // Licensing:
356 //
357 // This code is distributed under the GNU LGPL license.
358 //
359 // Modified:
360 //
361 // 27 March 2013
362 //
363 // Author:
364 //
365 // John Burkardt
366 //
367 // Parameters:
368 //
369 // Input, int G, the index of the generator.
370 // 0 <= G <= 31.
371 //
372 // Input, int CG1, CG2, the CG values for generator G.
373 //
374 {
375 int i;
376
377 i = +1;
378 cg_memory ( i, g, cg1, cg2 );
379
380 return;
381 }
382 //****************************************************************************80
383
cgn_get()384 int cgn_get ( )
385
386 //****************************************************************************80
387 //
388 // Purpose:
389 //
390 // CGN_GET gets the current generator index.
391 //
392 // Licensing:
393 //
394 // This code is distributed under the GNU LGPL license.
395 //
396 // Modified:
397 //
398 // 30 March 2013
399 //
400 // Author:
401 //
402 // John Burkardt
403 //
404 // Parameters:
405 //
406 // Output, int CGN_GET, the current generator index.
407 //
408 {
409 int g;
410 int i;
411
412 i = -1;
413 cgn_memory ( i, g );
414
415 return g;
416 }
417 //****************************************************************************80
418
cgn_memory(int i,int & g)419 void cgn_memory ( int i, int &g )
420
421 //****************************************************************************80
422 //
423 // Purpose:
424 //
425 // CGN_MEMORY stores the current generator index.
426 //
427 // Licensing:
428 //
429 // This code is distributed under the GNU LGPL license.
430 //
431 // Modified:
432 //
433 // 02 April 2013
434 //
435 // Author:
436 //
437 // John Burkardt
438 //
439 // Parameters:
440 //
441 // Input, int I, the desired action.
442 // -1, get the value.
443 // 0, initialize the value.
444 // 1, set the value.
445 //
446 // Input/output, int &G. For I = -1 or 0, this is output.
447 // For I = +1, this is input.
448 //
449 {
450 # define G_MAX 32
451
452 static int g_save = 0;
453 const int g_max = 32;
454
455 if ( i < 0 )
456 {
457 g = g_save;
458 }
459 else if ( i == 0 )
460 {
461 g_save = 0;
462 g = g_save;
463 }
464 else if ( 0 < i )
465 {
466
467 if ( g < 0 || g_max <= g )
468 {
469 cerr << "\n";
470 cerr << "CGN_MEMORY - Fatal error!\n";
471 cerr << " Input generator index G is out of bounds.\n";
472 exit ( 1 );
473 }
474
475 g_save = g;
476 }
477
478 return;
479 # undef G_MAX
480 }
481 //****************************************************************************80
482
cgn_set(int g)483 void cgn_set ( int g )
484
485 //****************************************************************************80
486 //
487 // Purpose:
488 //
489 // CGN_SET sets the current generator index.
490 //
491 // Licensing:
492 //
493 // This code is distributed under the GNU LGPL license.
494 //
495 // Modified:
496 //
497 // 30 March 2013
498 //
499 // Author:
500 //
501 // John Burkardt
502 //
503 // Parameters:
504 //
505 // Input, int G, the current generator index.
506 // 0 <= G <= 31.
507 //
508 {
509 int i;
510
511 i = +1;
512 cgn_memory ( i, g );
513
514 return;
515 }
516 //****************************************************************************80
517
get_state(int & cg1,int & cg2)518 void get_state ( int &cg1, int &cg2 )
519
520 //****************************************************************************80
521 //
522 // Purpose:
523 //
524 // GET_STATE returns the state of the current generator.
525 //
526 // Licensing:
527 //
528 // This code is distributed under the GNU LGPL license.
529 //
530 // Modified:
531 //
532 // 30 March 2013
533 //
534 // Author:
535 //
536 // Original Pascal version by Pierre L'Ecuyer, Serge Cote.
537 // C++ version by John Burkardt.
538 //
539 // Reference:
540 //
541 // Pierre LEcuyer, Serge Cote,
542 // Implementing a Random Number Package with Splitting Facilities,
543 // ACM Transactions on Mathematical Software,
544 // Volume 17, Number 1, March 1991, pages 98-111.
545 //
546 // Parameters:
547 //
548 // Output, int *CG1, *CG2, the CG values for the current generator.
549 //
550 {
551 int g;
552 //
553 // Check whether the package must be initialized.
554 //
555 if ( ! initialized_get ( ) )
556 {
557 cout << "\n";
558 cout << "GET_STATE - Note:\n";
559 cout << " Initializing RNGLIB package.\n";
560 initialize ( );
561 }
562 //
563 // Get the current generator index.
564 //
565 g = cgn_get ( );
566 //
567 // Retrieve the seed values for this generator.
568 //
569 cg_get ( g, cg1, cg2 );
570
571 return;
572 }
573 //****************************************************************************80
574
i4_uniform()575 int i4_uniform ( )
576
577 //****************************************************************************80
578 //
579 // Purpose:
580 //
581 // I4_UNIFORM generates a random positive integer.
582 //
583 // Discussion:
584 //
585 // This procedure returns a random integer following a uniform distribution
586 // over (1, 2147483562) using the current generator.
587 //
588 // The original name of this function was "random()", but this conflicts
589 // with a standard library function name in C.
590 //
591 // Licensing:
592 //
593 // This code is distributed under the GNU LGPL license.
594 //
595 // Modified:
596 //
597 // 02 April 2013
598 //
599 // Author:
600 //
601 // Original Pascal version by Pierre L'Ecuyer, Serge Cote.
602 // C++ version by John Burkardt.
603 //
604 // Reference:
605 //
606 // Pierre LEcuyer, Serge Cote,
607 // Implementing a Random Number Package with Splitting Facilities,
608 // ACM Transactions on Mathematical Software,
609 // Volume 17, Number 1, March 1991, pages 98-111.
610 //
611 // Parameters:
612 //
613 // Output, int I4_UNIFORM, the random integer.
614 //
615 {
616 const int a1 = 40014;
617 const int a2 = 40692;
618 int cg1;
619 int cg2;
620 int g;
621 int k;
622 const int m1 = 2147483563;
623 const int m2 = 2147483399;
624 bool value;
625 int z;
626 //
627 // Check whether the package must be initialized.
628 //
629 if ( ! initialized_get ( ) )
630 {
631 cout << "\n";
632 cout << "I4_UNIFORM - Note:\n";
633 cout << " Initializing RNGLIB package.\n";
634 initialize ( );
635 }
636 //
637 // Get the current generator index.
638 //
639 g = cgn_get ( );
640 //
641 // Retrieve the current seeds.
642 //
643 cg_get ( g, cg1, cg2 );
644 //
645 // Update the seeds.
646 //
647 k = cg1 / 53668;
648 cg1 = a1 * ( cg1 - k * 53668 ) - k * 12211;
649
650 if ( cg1 < 0 )
651 {
652 cg1 = cg1 + m1;
653 }
654
655 k = cg2 / 52774;
656 cg2 = a2 * ( cg2 - k * 52774 ) - k * 3791;
657
658 if ( cg2 < 0 )
659 {
660 cg2 = cg2 + m2;
661 }
662 //
663 // Store the updated seeds.
664 //
665 cg_set ( g, cg1, cg2 );
666 //
667 // Form the random integer.
668 //
669 z = cg1 - cg2;
670
671 if ( z < 1 )
672 {
673 z = z + m1 - 1;
674 }
675 //
676 // If the generator is antithetic, reflect the value.
677 //
678 value = antithetic_get ( );
679
680 if ( value )
681 {
682 z = m1 - z;
683 }
684 return z;
685 }
686 //****************************************************************************80
687
ig_get(int g,int & ig1,int & ig2)688 void ig_get ( int g, int &ig1, int &ig2 )
689
690 //****************************************************************************80
691 //
692 // Purpose:
693 //
694 // IG_GET queries the IG values for a given generator.
695 //
696 // Licensing:
697 //
698 // This code is distributed under the GNU LGPL license.
699 //
700 // Modified:
701 //
702 // 27 March 2013
703 //
704 // Author:
705 //
706 // John Burkardt
707 //
708 // Parameters:
709 //
710 // Input, int G, the index of the generator.
711 // 0 <= G <= 31.
712 //
713 // Output, int &IG1, &IG2, the IG values for generator G.
714 //
715 {
716 int i;
717
718 i = -1;
719 ig_memory ( i, g, ig1, ig2 );
720
721 return;
722 }
723 //****************************************************************************80
724
ig_memory(int i,int g,int & ig1,int & ig2)725 void ig_memory ( int i, int g, int &ig1, int &ig2 )
726
727 //****************************************************************************80
728 //
729 // Purpose:
730 //
731 // IG_MEMORY stores the IG values for all generators.
732 //
733 // Licensing:
734 //
735 // This code is distributed under the GNU LGPL license.
736 //
737 // Modified:
738 //
739 // 30 March 2013
740 //
741 // Author:
742 //
743 // John Burkardt
744 //
745 // Parameters:
746 //
747 // Input, int I, the desired action.
748 // -1, get a value.
749 // 0, initialize all values.
750 // 1, set a value.
751 //
752 // Input, int G, for I = -1 or +1, the index of
753 // the generator, with 0 <= G <= 31.
754 //
755 // Input/output, int &IG1, &IG2. For I = -1,
756 // these are output, for I = +1, these are input, for I = 0,
757 // these arguments are ignored. When used, the arguments are
758 // old or new values of the IG parameter for generator G.
759 //
760 {
761 # define G_MAX 32
762
763 const int g_max = 32;
764 static int ig1_save[G_MAX];
765 static int ig2_save[G_MAX];
766 int j;
767
768 if ( g < 0 || g_max <= g )
769 {
770 cerr << "\n";
771 cerr << "IG_MEMORY - Fatal error!\n";
772 cerr << " Input generator index G is out of bounds.\n";
773 exit ( 1 );
774 }
775
776 if ( i < 0 )
777 {
778 ig1 = ig1_save[g];
779 ig2 = ig2_save[g];
780 }
781 else if ( i == 0 )
782 {
783 for ( j = 0; j < g_max; j++ )
784 {
785 ig1_save[j] = 0;
786 ig2_save[j] = 0;
787 }
788 }
789 else if ( 0 < i )
790 {
791 ig1_save[g] = ig1;
792 ig2_save[g] = ig2;
793 }
794
795 return;
796 # undef G_MAX
797 }
798 //****************************************************************************80
799
ig_set(int g,int ig1,int ig2)800 void ig_set ( int g, int ig1, int ig2 )
801
802 //****************************************************************************80
803 //
804 // Purpose:
805 //
806 // IG_SET sets the IG values for a given generator.
807 //
808 // Licensing:
809 //
810 // This code is distributed under the GNU LGPL license.
811 //
812 // Modified:
813 //
814 // 27 March 2013
815 //
816 // Author:
817 //
818 // John Burkardt
819 //
820 // Parameters:
821 //
822 // Input, int G, the index of the generator.
823 // 0 <= G <= 31.
824 //
825 // Input, int IG1, IG2, the IG values for generator G.
826 //
827 {
828 int i;
829
830 i = +1;
831 ig_memory ( i, g, ig1, ig2 );
832
833 return;
834 }
835 //****************************************************************************80
836
init_generator(int t)837 void init_generator ( int t )
838
839 //****************************************************************************80
840 //
841 // Purpose:
842 //
843 // INIT_GENERATOR sets the state of generator G to initial, last or new seed.
844 //
845 // Licensing:
846 //
847 // This code is distributed under the GNU LGPL license.
848 //
849 // Modified:
850 //
851 // 02 April 2013
852 //
853 // Author:
854 //
855 // Original Pascal version by Pierre L'Ecuyer, Serge Cote.
856 // C++ version by John Burkardt.
857 //
858 // Reference:
859 //
860 // Pierre LEcuyer, Serge Cote,
861 // Implementing a Random Number Package with Splitting Facilities,
862 // ACM Transactions on Mathematical Software,
863 // Volume 17, Number 1, March 1991, pages 98-111.
864 //
865 // Parameters:
866 //
867 // Input, int T, the seed type:
868 // 0, use the seed chosen at initialization time.
869 // 1, use the last seed.
870 // 2, use a new seed set 2^30 values away.
871 //
872 {
873 const int a1_w = 1033780774;
874 const int a2_w = 1494757890;
875 int cg1;
876 int cg2;
877 int g;
878 int ig1;
879 int ig2;
880 int lg1;
881 int lg2;
882 const int m1 = 2147483563;
883 const int m2 = 2147483399;
884 //
885 // Check whether the package must be initialized.
886 //
887 if ( ! initialized_get ( ) )
888 {
889 cout << "\n";
890 cout << "INIT_GENERATOR - Note:\n";
891 cout << " Initializing RNGLIB package.\n";
892 initialize ( );
893 }
894 //
895 // Get the current generator index.
896 //
897 g = cgn_get ( );
898 //
899 // 0: restore the initial seed.
900 //
901 if ( t == 0 )
902 {
903 ig_get ( g, ig1, ig2 );
904 lg1 = ig1;
905 lg2 = ig2;
906 lg_set ( g, lg1, lg2 );
907 }
908 //
909 // 1: restore the last seed.
910 //
911 else if ( t == 1 )
912 {
913 lg_get ( g, lg1, lg2 );
914 }
915 //
916 // 2: Advance to a new seed.
917 //
918 else if ( t == 2 )
919 {
920 lg_get ( g, lg1, lg2 );
921 lg1 = multmod ( a1_w, lg1, m1 );
922 lg2 = multmod ( a2_w, lg2, m2 );
923 lg_set ( g, lg1, lg2 );
924 }
925 else
926 {
927 cerr << "\n";
928 cerr << "INIT_GENERATOR - Fatal error!\n";
929 cerr << " Input parameter T out of bounds.\n";
930 exit ( 1 );
931 }
932 //
933 // Store the new seed.
934 //
935 cg1 = lg1;
936 cg2 = lg2;
937 cg_set ( g, cg1, cg2 );
938
939 return;
940 }
941 //****************************************************************************80
942
initialize()943 void initialize ( )
944
945 //****************************************************************************80
946 //
947 // Purpose:
948 //
949 // INITIALIZE initializes the random number generator library.
950 //
951 // Licensing:
952 //
953 // This code is distributed under the GNU LGPL license.
954 //
955 // Modified:
956 //
957 // 30 March 2013
958 //
959 // Author:
960 //
961 // Original Pascal version by Pierre L'Ecuyer, Serge Cote.
962 // C++ version by John Burkardt.
963 //
964 // Reference:
965 //
966 // Pierre LEcuyer, Serge Cote,
967 // Implementing a Random Number Package with Splitting Facilities,
968 // ACM Transactions on Mathematical Software,
969 // Volume 17, Number 1, March 1991, pages 98-111.
970 //
971 // Parameters:
972 //
973 // None
974 //
975 {
976 int g;
977 const int g_max = 32;
978 int ig1;
979 int ig2;
980 int value;
981 //
982 // Remember that we have called INITIALIZE().
983 //
984 initialized_set ( );
985 //
986 // Initialize all generators to have FALSE antithetic value.
987 //
988 value = 0;
989 for ( g = 0; g < g_max; g++ )
990 {
991 cgn_set ( g );
992 antithetic_set ( value );
993 }
994 //
995 // Set the initial seeds.
996 //
997 ig1 = 1234567890;
998 ig2 = 123456789;
999 set_initial_seed ( ig1, ig2 );
1000 //
1001 // Initialize the current generator index to 0.
1002 //
1003 g = 0;
1004 cgn_set ( g );
1005
1006 cout << "\n";
1007 cout << "INITIALIZE - Note:\n";
1008 cout << " The RNGLIB package has been initialized.\n";
1009
1010 return;
1011 }
1012 //****************************************************************************80
1013
initialized_get()1014 bool initialized_get ( )
1015
1016 //****************************************************************************80
1017 //
1018 // Purpose:
1019 //
1020 // INITIALIZED_GET queries the INITIALIZED value.
1021 //
1022 // Licensing:
1023 //
1024 // This code is distributed under the GNU LGPL license.
1025 //
1026 // Modified:
1027 //
1028 // 28 March 2013
1029 //
1030 // Author:
1031 //
1032 // John Burkardt
1033 //
1034 // Parameters:
1035 //
1036 // Output, bool INITIALIZED_GET, is TRUE (1) if the package has been initialized.
1037 //
1038 {
1039 int i;
1040 bool value;
1041
1042 i = -1;
1043 initialized_memory ( i, value );
1044
1045 return value;
1046 }
1047 //****************************************************************************80
1048
initialized_memory(int i,bool & initialized)1049 void initialized_memory ( int i, bool &initialized )
1050
1051 //****************************************************************************80
1052 //
1053 // Purpose:
1054 //
1055 // INITIALIZED_MEMORY stores the INITIALIZED value for the package.
1056 //
1057 // Licensing:
1058 //
1059 // This code is distributed under the GNU LGPL license.
1060 //
1061 // Modified:
1062 //
1063 // 28 March 2013
1064 //
1065 // Author:
1066 //
1067 // John Burkardt
1068 //
1069 // Parameters:
1070 //
1071 // Input, int I, the desired action.
1072 // -1, get the value.
1073 // 0, initialize the value.
1074 // 1, set the value.
1075 //
1076 // Input/output, bool &INITIALIZED. For I = -1, this is an output
1077 // quantity. If I = +1, this is an input quantity. If I = 0,
1078 // this is ignored.
1079 //
1080 {
1081 static bool initialized_save = false;
1082
1083 if ( i < 0 )
1084 {
1085 initialized = initialized_save;
1086 }
1087 else if ( i == 0 )
1088 {
1089 initialized_save = false;
1090 }
1091 else if ( 0 < i )
1092 {
1093 initialized_save = initialized;
1094 }
1095
1096 return;
1097 }
1098 //****************************************************************************80
1099
initialized_set()1100 void initialized_set ( )
1101
1102 //****************************************************************************80
1103 //
1104 // Purpose:
1105 //
1106 // INITIALIZED_SET sets the INITIALIZED value true.
1107 //
1108 // Licensing:
1109 //
1110 // This code is distributed under the GNU LGPL license.
1111 //
1112 // Modified:
1113 //
1114 // 28 March 2013
1115 //
1116 // Author:
1117 //
1118 // John Burkardt
1119 //
1120 // Parameters:
1121 //
1122 // None
1123 //
1124 {
1125 int i;
1126 bool initialized;
1127
1128 i = +1;
1129 initialized = true;
1130 initialized_memory ( i, initialized );
1131
1132 return;
1133 }
1134 //****************************************************************************80
1135
lg_get(int g,int & lg1,int & lg2)1136 void lg_get ( int g, int &lg1, int &lg2 )
1137
1138 //****************************************************************************80
1139 //
1140 // Purpose:
1141 //
1142 // LG_GET queries the LG values for a given generator.
1143 //
1144 // Licensing:
1145 //
1146 // This code is distributed under the GNU LGPL license.
1147 //
1148 // Modified:
1149 //
1150 // 27 March 2013
1151 //
1152 // Author:
1153 //
1154 // John Burkardt
1155 //
1156 // Parameters:
1157 //
1158 // Input, int G, the index of the generator.
1159 // 0 <= G <= 31.
1160 //
1161 // Output, int &LG1, &LG2, the LG values for generator G.
1162 //
1163 {
1164 int i;
1165
1166 i = -1;
1167 lg_memory ( i, g, lg1, lg2 );
1168
1169 return;
1170 }
1171 //****************************************************************************80
1172
lg_memory(int i,int g,int & lg1,int & lg2)1173 void lg_memory ( int i, int g, int &lg1, int &lg2 )
1174
1175 //****************************************************************************80
1176 //
1177 // Purpose:
1178 //
1179 // LG_MEMORY stores the LG values for all generators.
1180 //
1181 // Licensing:
1182 //
1183 // This code is distributed under the GNU LGPL license.
1184 //
1185 // Modified:
1186 //
1187 // 30 March 2013
1188 //
1189 // Author:
1190 //
1191 // John Burkardt
1192 //
1193 // Parameters:
1194 //
1195 // Input, int I, the desired action.
1196 // -1, get a value.
1197 // 0, initialize all values.
1198 // 1, set a value.
1199 //
1200 // Input, int G, for I = -1 or +1, the index of
1201 // the generator, with 0 <= G <= 31.
1202 //
1203 // Input/output, int &LG1, &LG2. For I = -1,
1204 // these are output, for I = +1, these are input, for I = 0,
1205 // these arguments are ignored. When used, the arguments are
1206 // old or new values of the LG parameter for generator G.
1207 //
1208 {
1209 # define G_MAX 32
1210
1211 const int g_max = 32;
1212
1213 int j;
1214 static int lg1_save[G_MAX];
1215 static int lg2_save[G_MAX];
1216
1217 if ( g < 0 || g_max <= g )
1218 {
1219 cerr << "\n";
1220 cerr << "LG_MEMORY - Fatal error!\n";
1221 cerr << " Input generator index G is out of bounds.\n";
1222 exit ( 1 );
1223 }
1224
1225 if ( i < 0 )
1226 {
1227 lg1 = lg1_save[g];
1228 lg2 = lg2_save[g];
1229 }
1230 else if ( i == 0 )
1231 {
1232 for ( j = 0; j < g_max; j++ )
1233 {
1234 lg1_save[j] = 0;
1235 lg2_save[j] = 0;
1236 }
1237 }
1238 else if ( 0 < i )
1239 {
1240 lg1_save[g] = lg1;
1241 lg2_save[g] = lg2;
1242 }
1243 return;
1244 # undef G_MAX
1245 }
1246 //****************************************************************************80
1247
lg_set(int g,int lg1,int lg2)1248 void lg_set ( int g, int lg1, int lg2 )
1249
1250 //****************************************************************************80
1251 //
1252 // Purpose:
1253 //
1254 // LG_SET sets the LG values for a given generator.
1255 //
1256 // Licensing:
1257 //
1258 // This code is distributed under the GNU LGPL license.
1259 //
1260 // Modified:
1261 //
1262 // 27 March 2013
1263 //
1264 // Author:
1265 //
1266 // John Burkardt
1267 //
1268 // Parameters:
1269 //
1270 // Input, int G, the index of the generator.
1271 // 0 <= G <= 31.
1272 //
1273 // Input, int LG1, LG2, the LG values for generator G.
1274 //
1275 {
1276 int i;
1277
1278 i = +1;
1279 lg_memory ( i, g, lg1, lg2 );
1280
1281 return;
1282 }
1283 //****************************************************************************80
1284
multmod(int a,int s,int m)1285 int multmod ( int a, int s, int m )
1286
1287 //****************************************************************************80
1288 //
1289 // Purpose:
1290 //
1291 // MULTMOD carries out modular multiplication.
1292 //
1293 // Discussion:
1294 //
1295 // This procedure returns
1296 //
1297 // ( A * S ) mod M
1298 //
1299 // Licensing:
1300 //
1301 // This code is distributed under the GNU LGPL license.
1302 //
1303 // Modified:
1304 //
1305 // 27 March 2013
1306 //
1307 // Author:
1308 //
1309 // Original Pascal version by Pierre L'Ecuyer, Serge Cote.
1310 // C++ version by John Burkardt.
1311 //
1312 // Reference:
1313 //
1314 // Pierre LEcuyer, Serge Cote,
1315 // Implementing a Random Number Package with Splitting Facilities,
1316 // ACM Transactions on Mathematical Software,
1317 // Volume 17, Number 1, March 1991, pages 98-111.
1318 //
1319 // Parameters:
1320 //
1321 // Input, int A, S, M, the arguments.
1322 //
1323 // Output, int MULTMOD, the value of the product of A and S,
1324 // modulo M.
1325 //
1326 {
1327 int a0;
1328 int a1;
1329 const int h = 32768;
1330 int k;
1331 int p;
1332 int q;
1333 int qh;
1334 int rh;
1335
1336 if ( a <= 0 )
1337 {
1338 cerr << "\n";
1339 cerr << "MULTMOD - Fatal error!\n";
1340 cerr << " A <= 0.\n";
1341 exit ( 1 );
1342 }
1343
1344 if ( m <= a )
1345 {
1346 cerr << "\n";
1347 cerr << "MULTMOD - Fatal error!\n";
1348 cerr << " M <= A.\n";
1349 exit ( 1 );
1350 }
1351
1352 if ( s <= 0 )
1353 {
1354 cerr << "\n";
1355 cerr << "MULTMOD - Fatal error!\n";
1356 cerr << " S <= 0.\n";
1357 exit ( 1 );
1358 }
1359
1360 if ( m <= s )
1361 {
1362 cerr << "\n";
1363 cerr << "MULTMOD - Fatal error!\n";
1364 cerr << " M <= S.\n";
1365 exit ( 1 );
1366 }
1367
1368 if ( a < h )
1369 {
1370 a0 = a;
1371 p = 0;
1372 }
1373 else
1374 {
1375 a1 = a / h;
1376 a0 = a - h * a1;
1377 qh = m / h;
1378 rh = m - h * qh;
1379
1380 if ( h <= a1 )
1381 {
1382 a1 = a1 - h;
1383 k = s / qh;
1384 p = h * ( s - k * qh ) - k * rh;
1385
1386 while ( p < 0 )
1387 {
1388 p = p + m;
1389 }
1390 }
1391 else
1392 {
1393 p = 0;
1394 }
1395
1396 if ( a1 != 0 )
1397 {
1398 q = m / a1;
1399 k = s / q;
1400 p = p - k * ( m - a1 * q );
1401
1402 if ( 0 < p )
1403 {
1404 p = p - m;
1405 }
1406
1407 p = p + a1 * ( s - k * q );
1408
1409 while ( p < 0 )
1410 {
1411 p = p + m;
1412 }
1413 }
1414
1415 k = p / qh;
1416 p = h * ( p - k * qh ) - k * rh;
1417
1418 while ( p < 0 )
1419 {
1420 p = p + m;
1421 }
1422 }
1423
1424 if ( a0 != 0 )
1425 {
1426 q = m / a0;
1427 k = s / q;
1428 p = p - k * ( m - a0 * q );
1429
1430 if ( 0 < p )
1431 {
1432 p = p - m;
1433 }
1434
1435 p = p + a0 * ( s - k * q );
1436
1437 while ( p < 0 )
1438 {
1439 p = p + m;
1440 }
1441 }
1442 return p;
1443 }
1444 //****************************************************************************80
1445
r4_uniform_01()1446 float r4_uniform_01 ( )
1447
1448 //****************************************************************************80
1449 //
1450 // Purpose:
1451 //
1452 // R4_UNIFORM_01 returns a uniform random real number in [0,1].
1453 //
1454 // Discussion:
1455 //
1456 // This procedure returns a random floating point number from a uniform
1457 // distribution over (0,1), not including the endpoint values, using the
1458 // current random number generator.
1459 //
1460 // Licensing:
1461 //
1462 // This code is distributed under the GNU LGPL license.
1463 //
1464 // Modified:
1465 //
1466 // 10 April 2013
1467 //
1468 // Author:
1469 //
1470 // Original Pascal version by Pierre L'Ecuyer, Serge Cote.
1471 // C++ version by John Burkardt.
1472 //
1473 // Reference:
1474 //
1475 // Pierre LEcuyer, Serge Cote,
1476 // Implementing a Random Number Package with Splitting Facilities,
1477 // ACM Transactions on Mathematical Software,
1478 // Volume 17, Number 1, March 1991, pages 98-111.
1479 //
1480 // Parameters:
1481 //
1482 // Output, float R4_UNIFORM_01, a uniform random value in [0,1].
1483 //
1484 {
1485 int i;
1486 float value;
1487 //
1488 // Check whether the package must be initialized.
1489 //
1490 if ( ! initialized_get ( ) )
1491 {
1492 cout << "\n";
1493 cout << "R4_UNIFORM_01 - Note:\n";
1494 cout << " Initializing RNGLIB package.\n";
1495 initialize ( );
1496 }
1497 //
1498 // Get a random integer.
1499 //
1500 i = i4_uniform ( );
1501 //
1502 // Scale it to [0,1].
1503 //
1504 value = ( float ) ( i ) * 4.656613057E-10;
1505
1506 return value;
1507 }
1508 //****************************************************************************80
1509
r8_uniform_01()1510 double r8_uniform_01 ( )
1511
1512 //****************************************************************************80
1513 //
1514 // Purpose:
1515 //
1516 // R8_UNIFORM_01 returns a uniform random double in [0,1].
1517 //
1518 // Discussion:
1519 //
1520 // This procedure returns a random valule from a uniform
1521 // distribution over (0,1), not including the endpoint values, using the
1522 // current random number generator.
1523 //
1524 // Licensing:
1525 //
1526 // This code is distributed under the GNU LGPL license.
1527 //
1528 // Modified:
1529 //
1530 // 10 April 2013
1531 //
1532 // Author:
1533 //
1534 // Original Pascal version by Pierre L'Ecuyer, Serge Cote.
1535 // C++ version by John Burkardt.
1536 //
1537 // Reference:
1538 //
1539 // Pierre LEcuyer, Serge Cote,
1540 // Implementing a Random Number Package with Splitting Facilities,
1541 // ACM Transactions on Mathematical Software,
1542 // Volume 17, Number 1, March 1991, pages 98-111.
1543 //
1544 // Parameters:
1545 //
1546 // Output, double R8_UNIFORM_01, a uniform random value in [0,1].
1547 //
1548 {
1549 int i;
1550 double value;
1551 //
1552 // Check whether the package must be initialized.
1553 //
1554 if ( ! initialized_get ( ) )
1555 {
1556 cout << "\n";
1557 cout << "R8_UNIFORM_01 - Note:\n";
1558 cout << " Initializing RNGLIB package.\n";
1559 initialize ( );
1560 }
1561 //
1562 // Get a random integer.
1563 //
1564 i = i4_uniform ( );
1565 //
1566 // Scale it to [0,1].
1567 //
1568 value = ( double ) ( i ) * 4.656613057E-10;
1569
1570 return value;
1571 }
1572 //****************************************************************************80
1573
set_initial_seed(int ig1,int ig2)1574 void set_initial_seed ( int ig1, int ig2 )
1575
1576 //****************************************************************************80
1577 //
1578 // Purpose:
1579 //
1580 // SET_INITIAL_SEED resets the initial seed and state for all generators.
1581 //
1582 // Licensing:
1583 //
1584 // This code is distributed under the GNU LGPL license.
1585 //
1586 // Modified:
1587 //
1588 // 27 March 2013
1589 //
1590 // Author:
1591 //
1592 // Original Pascal version by Pierre L'Ecuyer, Serge Cote.
1593 // C++ version by John Burkardt.
1594 //
1595 // Reference:
1596 //
1597 // Pierre LEcuyer, Serge Cote,
1598 // Implementing a Random Number Package with Splitting Facilities,
1599 // ACM Transactions on Mathematical Software,
1600 // Volume 17, Number 1, March 1991, pages 98-111.
1601 //
1602 // Parameters:
1603 //
1604 // Input, int IG1, IG2, the initial seed values
1605 // for the first generator.
1606 // 1 <= IG1 < 2147483563
1607 // 1 <= IG2 < 2147483399
1608 //
1609 {
1610 const int a1_vw = 2082007225;
1611 const int a2_vw = 784306273;
1612 int g;
1613 const int g_max = 32;
1614 int i;
1615 const int m1 = 2147483563;
1616 const int m2 = 2147483399;
1617 int t;
1618
1619 if ( ig1 < 1 || m1 <= ig1 )
1620 {
1621 cerr << "\n";
1622 cerr << "SET_INITIAL_SEED - Fatal error!\n";
1623 cerr << " Input parameter IG1 out of bounds.\n";
1624 exit ( 1 );
1625 }
1626
1627 if ( ig2 < 1 || m2 <= ig2 )
1628 {
1629 cerr << "\n";
1630 cerr << "SET_INITIAL_SEED - Fatal error!\n";
1631 cerr << " Input parameter IG2 out of bounds.\n";
1632 exit ( 1 );
1633 }
1634 //
1635 // Because INITIALIZE calls SET_INITIAL_SEED, it's not easy to correct
1636 // the error that arises if SET_INITIAL_SEED is called before INITIALIZE.
1637 // So don't bother trying.
1638 //
1639 if ( ! initialized_get ( ) )
1640 {
1641 cout << "\n";
1642 cout << "SET_INITIAL_SEED - Fatal error!\n";
1643 cout << " The RNGLIB package has not been initialized.\n";
1644 exit ( 1 );
1645 }
1646 //
1647 // Set the initial seed, then initialize the first generator.
1648 //
1649 g = 0;
1650 cgn_set ( g );
1651
1652 ig_set ( g, ig1, ig2 );
1653
1654 t = 0;
1655 init_generator ( t );
1656 //
1657 // Now do similar operations for the other generators.
1658 //
1659 for ( g = 1; g < g_max; g++ )
1660 {
1661 cgn_set ( g );
1662 ig1 = multmod ( a1_vw, ig1, m1 );
1663 ig2 = multmod ( a2_vw, ig2, m2 );
1664 ig_set ( g, ig1, ig2 );
1665 init_generator ( t );
1666 }
1667 //
1668 // Now choose the first generator.
1669 //
1670 g = 0;
1671 cgn_set ( g );
1672
1673 return;
1674 }
1675 //****************************************************************************80
1676
set_seed(int cg1,int cg2)1677 void set_seed ( int cg1, int cg2 )
1678
1679 //****************************************************************************80
1680 //
1681 // Purpose:
1682 //
1683 // SET_SEED resets the initial seed and the state of generator G.
1684 //
1685 // Licensing:
1686 //
1687 // This code is distributed under the GNU LGPL license.
1688 //
1689 // Modified:
1690 //
1691 // 02 April 2013
1692 //
1693 // Author:
1694 //
1695 // Original Pascal version by Pierre L'Ecuyer, Serge Cote.
1696 // C++ version by John Burkardt.
1697 //
1698 // Reference:
1699 //
1700 // Pierre LEcuyer, Serge Cote,
1701 // Implementing a Random Number Package with Splitting Facilities,
1702 // ACM Transactions on Mathematical Software,
1703 // Volume 17, Number 1, March 1991, pages 98-111.
1704 //
1705 // Parameters:
1706 //
1707 // Input, int CG1, CG2, the CG values for generator G.
1708 // 1 <= CG1 < 2147483563
1709 // 1 <= CG2 < 2147483399
1710 //
1711 {
1712 int g;
1713 int i;
1714 const int m1 = 2147483563;
1715 const int m2 = 2147483399;
1716 int t;
1717
1718 if ( cg1 < 1 || m1 <= cg1 )
1719 {
1720 cerr << "\n";
1721 cerr << "SET_SEED - Fatal error!\n";
1722 cerr << " Input parameter CG1 out of bounds.\n";
1723 exit ( 1 );
1724 }
1725
1726 if ( cg2 < 1 || m2 <= cg2 )
1727 {
1728 cerr << "\n";
1729 cerr << "SET_SEED - Fatal error!\n";
1730 cerr << " Input parameter CG2 out of bounds.\n";
1731 exit ( 1 );
1732 }
1733 //
1734 // Check whether the package must be initialized.
1735 //
1736 if ( ! initialized_get ( ) )
1737 {
1738 cout << "\n";
1739 cout << "SET_SEED - Note:\n";
1740 cout << " Initializing RNGLIB package.\n";
1741 initialize ( );
1742 }
1743 //
1744 // Retrieve the current generator index.
1745 //
1746 g = cgn_get ( );
1747 //
1748 // Set the seeds.
1749 //
1750 cg_set ( g, cg1, cg2 );
1751 //
1752 // Initialize the generator.
1753 //
1754 t = 0;
1755 init_generator ( t );
1756
1757 return;
1758 }
1759 //****************************************************************************80
1760
timestamp()1761 void timestamp ( )
1762
1763 //****************************************************************************80
1764 //
1765 // Purpose:
1766 //
1767 // TIMESTAMP prints the current YMDHMS date as a time stamp.
1768 //
1769 // Example:
1770 //
1771 // 31 May 2001 09:45:54 AM
1772 //
1773 // Licensing:
1774 //
1775 // This code is distributed under the GNU LGPL license.
1776 //
1777 // Modified:
1778 //
1779 // 08 July 2009
1780 //
1781 // Author:
1782 //
1783 // John Burkardt
1784 //
1785 // Parameters:
1786 //
1787 // None
1788 //
1789 {
1790 # define TIME_SIZE 40
1791
1792 static char time_buffer[TIME_SIZE];
1793 const struct std::tm *tm_ptr;
1794 size_t len;
1795 std::time_t now;
1796
1797 now = std::time ( NULL );
1798 tm_ptr = std::localtime ( &now );
1799
1800 len = std::strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr );
1801
1802 std::cout << time_buffer << "\n";
1803
1804 return;
1805 # undef TIME_SIZE
1806 }
1807
1808 }
1809