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