1 // SPDX-License-Identifier: Apache-2.0
2 //
3 // Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au)
4 // Copyright 2008-2016 National ICT Australia (NICTA)
5 //
6 // Licensed under the Apache License, Version 2.0 (the "License");
7 // you may not use this file except in compliance with the License.
8 // You may obtain a copy of the License at
9 // http://www.apache.org/licenses/LICENSE-2.0
10 //
11 // Unless required by applicable law or agreed to in writing, software
12 // distributed under the License is distributed on an "AS IS" BASIS,
13 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 // See the License for the specific language governing permissions and
15 // limitations under the License.
16 // ------------------------------------------------------------------------
17
18
19 //! \addtogroup arma_rng
20 //! @{
21
22
23 #if defined(ARMA_RNG_ALT)
24 #undef ARMA_USE_EXTERN_RNG
25 #endif
26
27
28 // NOTE: mt19937_64_instance_warmup is used as a workaround
29 // NOTE: for thread_local issue on macOS 11 and/or AppleClang 12.0
30 // NOTE: see https://gitlab.com/conradsnicta/armadillo-code/-/issues/173
31 // NOTE: if this workaround causes problems, please report it and
32 // NOTE: disable the workaround by uncommenting the code block below:
33
34 // #if defined(__APPLE__) || defined(__apple_build_version__)
35 // #if !defined(ARMA_DONT_DISABLE_EXTERN_RNG)
36 // #undef ARMA_USE_EXTERN_RNG
37 // #endif
38 // #endif
39
40
41 // NOTE: workaround for another thread_local issue on macOS
42 // NOTE: where GCC (not Clang) may not have support for thread_local
43
44 #if (defined(__APPLE__) && defined(__GNUG__) && !defined(__clang__))
45 #if !defined(ARMA_DONT_DISABLE_EXTERN_RNG)
46 #undef ARMA_USE_EXTERN_RNG
47 #endif
48 #endif
49
50
51
52 #if defined(ARMA_USE_EXTERN_RNG)
53 extern thread_local std::mt19937_64 mt19937_64_instance;
54
55 #if defined(__APPLE__) || defined(__apple_build_version__)
56 namespace
57 {
58 struct mt19937_64_instance_warmup
59 {
mt19937_64_instance_warmup__anon48fe36020111::mt19937_64_instance_warmup60 inline mt19937_64_instance_warmup()
61 {
62 typename std::mt19937_64::result_type junk = mt19937_64_instance();
63 arma_ignore(junk);
64 }
65 };
66
67 static mt19937_64_instance_warmup mt19937_64_instance_warmup_run;
68 }
69 #endif
70 #endif
71
72
73
74 class arma_rng
75 {
76 public:
77
78 #if defined(ARMA_RNG_ALT)
79 typedef arma_rng_alt::seed_type seed_type;
80 #elif defined(ARMA_USE_EXTERN_RNG)
81 typedef std::mt19937_64::result_type seed_type;
82 #else
83 typedef arma_rng_cxx98::seed_type seed_type;
84 #endif
85
86 #if defined(ARMA_RNG_ALT)
87 static constexpr int rng_method = 2;
88 #elif defined(ARMA_USE_EXTERN_RNG)
89 static constexpr int rng_method = 1;
90 #else
91 static constexpr int rng_method = 0;
92 #endif
93
94 inline static void set_seed(const seed_type val);
95 inline static void set_seed_random();
96
97 template<typename eT> struct randi;
98 template<typename eT> struct randu;
99 template<typename eT> struct randn;
100 template<typename eT> struct randg;
101 };
102
103
104
105 inline
106 void
set_seed(const arma_rng::seed_type val)107 arma_rng::set_seed(const arma_rng::seed_type val)
108 {
109 #if defined(ARMA_RNG_ALT)
110 {
111 arma_rng_alt::set_seed(val);
112 }
113 #elif defined(ARMA_USE_EXTERN_RNG)
114 {
115 mt19937_64_instance.seed(val);
116 }
117 #else
118 {
119 arma_rng_cxx98::set_seed(val);
120 }
121 #endif
122 }
123
124
125
126 arma_cold
127 inline
128 void
set_seed_random()129 arma_rng::set_seed_random()
130 {
131 seed_type seed1 = seed_type(0);
132 seed_type seed2 = seed_type(0);
133 seed_type seed3 = seed_type(0);
134 seed_type seed4 = seed_type(0);
135
136 bool have_seed = false;
137
138 try
139 {
140 std::random_device rd;
141
142 if(rd.entropy() > double(0)) { seed1 = static_cast<seed_type>( rd() ); }
143
144 if(seed1 != seed_type(0)) { have_seed = true; }
145 }
146 catch(...) {}
147
148
149 if(have_seed == false)
150 {
151 try
152 {
153 union
154 {
155 seed_type a;
156 unsigned char b[sizeof(seed_type)];
157 } tmp;
158
159 tmp.a = seed_type(0);
160
161 std::ifstream f("/dev/urandom", std::ifstream::binary);
162
163 if(f.good()) { f.read((char*)(&(tmp.b[0])), sizeof(seed_type)); }
164
165 if(f.good())
166 {
167 seed2 = tmp.a;
168
169 if(seed2 != seed_type(0)) { have_seed = true; }
170 }
171 }
172 catch(...) {}
173 }
174
175
176 if(have_seed == false)
177 {
178 // get better-than-nothing seeds in case reading /dev/urandom failed
179
180 const std::chrono::system_clock::time_point tp_now = std::chrono::system_clock::now();
181
182 auto since_epoch_usec = std::chrono::duration_cast<std::chrono::microseconds>(tp_now.time_since_epoch()).count();
183
184 seed3 = static_cast<seed_type>( since_epoch_usec & 0xFFFF );
185
186 union
187 {
188 uword* a;
189 unsigned char b[sizeof(uword*)];
190 } tmp;
191
192 tmp.a = (uword*)malloc(sizeof(uword));
193
194 if(tmp.a != nullptr)
195 {
196 for(size_t i=0; i<sizeof(uword*); ++i) { seed4 += seed_type(tmp.b[i]); }
197
198 free(tmp.a);
199 }
200 }
201
202 arma_rng::set_seed( seed1 + seed2 + seed3 + seed4 );
203 }
204
205
206
207 //
208
209
210
211 template<typename eT>
212 struct arma_rng::randi
213 {
214 inline
operator eTarma_rng::randi215 operator eT ()
216 {
217 #if defined(ARMA_RNG_ALT)
218 {
219 return eT( arma_rng_alt::randi_val() );
220 }
221 #elif defined(ARMA_USE_EXTERN_RNG)
222 {
223 constexpr double scale = double(std::numeric_limits<int>::max()) / double(std::mt19937_64::max());
224
225 return eT( double(mt19937_64_instance()) * scale );
226 }
227 #else
228 {
229 return eT( arma_rng_cxx98::randi_val() );
230 }
231 #endif
232 }
233
234
235 inline
236 static
237 int
max_valarma_rng::randi238 max_val()
239 {
240 #if defined(ARMA_RNG_ALT)
241 {
242 return arma_rng_alt::randi_max_val();
243 }
244 #elif defined(ARMA_USE_EXTERN_RNG)
245 {
246 return std::numeric_limits<int>::max();
247 }
248 #else
249 {
250 return arma_rng_cxx98::randi_max_val();
251 }
252 #endif
253 }
254
255
256 inline
257 static
258 void
fillarma_rng::randi259 fill(eT* mem, const uword N, const int a, const int b)
260 {
261 #if defined(ARMA_RNG_ALT)
262 {
263 arma_rng_alt::randi_fill(mem, N, a, b);
264 }
265 #elif defined(ARMA_USE_EXTERN_RNG)
266 {
267 std::uniform_int_distribution<int> local_i_distr(a, b);
268
269 for(uword i=0; i<N; ++i) { mem[i] = eT(local_i_distr(mt19937_64_instance)); }
270 }
271 #else
272 {
273 if(N == uword(1)) { arma_rng_cxx98::randi_fill(mem, uword(1), a, b); return; }
274
275 typedef typename std::mt19937_64::result_type local_seed_type;
276
277 std::mt19937_64 local_engine;
278 std::uniform_int_distribution<int> local_i_distr(a, b);
279
280 local_engine.seed( local_seed_type(std::rand()) );
281
282 for(uword i=0; i<N; ++i) { mem[i] = eT(local_i_distr(local_engine)); }
283 }
284 #endif
285 }
286 };
287
288
289
290 //
291
292
293
294 template<typename eT>
295 struct arma_rng::randu
296 {
297 inline
operator eTarma_rng::randu298 operator eT ()
299 {
300 #if defined(ARMA_RNG_ALT)
301 {
302 return eT( arma_rng_alt::randu_val() );
303 }
304 #elif defined(ARMA_USE_EXTERN_RNG)
305 {
306 constexpr double scale = double(1.0) / double(std::mt19937_64::max());
307
308 return eT( double(mt19937_64_instance()) * scale );
309 }
310 #else
311 {
312 return eT( arma_rng_cxx98::randu_val() );
313 }
314 #endif
315 }
316
317
318 inline
319 static
320 void
fillarma_rng::randu321 fill(eT* mem, const uword N)
322 {
323 #if defined(ARMA_RNG_ALT)
324 {
325 for(uword i=0; i < N; ++i) { mem[i] = eT( arma_rng_alt::randu_val() ); }
326 }
327 #elif defined(ARMA_USE_EXTERN_RNG)
328 {
329 std::uniform_real_distribution<double> local_u_distr;
330
331 for(uword i=0; i < N; ++i) { mem[i] = eT( local_u_distr(mt19937_64_instance) ); }
332 }
333 #else
334 {
335 if(N == uword(1)) { mem[0] = eT( arma_rng_cxx98::randu_val() ); return; }
336
337 typedef typename std::mt19937_64::result_type local_seed_type;
338
339 std::mt19937_64 local_engine;
340 std::uniform_real_distribution<double> local_u_distr;
341
342 local_engine.seed( local_seed_type(std::rand()) );
343
344 for(uword i=0; i < N; ++i) { mem[i] = eT( local_u_distr(local_engine) ); }
345 }
346 #endif
347 }
348 };
349
350
351
352 template<typename T>
353 struct arma_rng::randu< std::complex<T> >
354 {
355 arma_inline
operator std::complex<T>arma_rng::randu356 operator std::complex<T> ()
357 {
358 #if defined(ARMA_RNG_ALT)
359 {
360 const T a = T( arma_rng_alt::randu_val() );
361 const T b = T( arma_rng_alt::randu_val() );
362
363 return std::complex<T>(a, b);
364 }
365 #elif defined(ARMA_USE_EXTERN_RNG)
366 {
367 std::uniform_real_distribution<double> local_u_distr;
368
369 const T a = T( local_u_distr(mt19937_64_instance) );
370 const T b = T( local_u_distr(mt19937_64_instance) );
371
372 return std::complex<T>(a, b);
373 }
374 #else
375 {
376 const T a = T( arma_rng_cxx98::randu_val() );
377 const T b = T( arma_rng_cxx98::randu_val() );
378
379 return std::complex<T>(a, b);
380 }
381 #endif
382 }
383
384
385 inline
386 static
387 void
fillarma_rng::randu388 fill(std::complex<T>* mem, const uword N)
389 {
390 #if defined(ARMA_RNG_ALT)
391 {
392 for(uword i=0; i < N; ++i)
393 {
394 const T a = T( arma_rng_alt::randu_val() );
395 const T b = T( arma_rng_alt::randu_val() );
396
397 mem[i] = std::complex<T>(a, b);
398 }
399 }
400 #elif defined(ARMA_USE_EXTERN_RNG)
401 {
402 std::uniform_real_distribution<double> local_u_distr;
403
404 for(uword i=0; i < N; ++i)
405 {
406 const T a = T( local_u_distr(mt19937_64_instance) );
407 const T b = T( local_u_distr(mt19937_64_instance) );
408
409 mem[i] = std::complex<T>(a, b);
410 }
411 }
412 #else
413 {
414 if(N == uword(1))
415 {
416 const T a = T( arma_rng_cxx98::randu_val() );
417 const T b = T( arma_rng_cxx98::randu_val() );
418
419 mem[0] = std::complex<T>(a, b);
420
421 return;
422 }
423
424 typedef typename std::mt19937_64::result_type local_seed_type;
425
426 std::mt19937_64 local_engine;
427 std::uniform_real_distribution<double> local_u_distr;
428
429 local_engine.seed( local_seed_type(std::rand()) );
430
431 for(uword i=0; i < N; ++i)
432 {
433 const T a = T( local_u_distr(local_engine) );
434 const T b = T( local_u_distr(local_engine) );
435
436 mem[i] = std::complex<T>(a, b);
437 }
438 }
439 #endif
440 }
441 };
442
443
444
445 //
446
447
448
449 template<typename eT>
450 struct arma_rng::randn
451 {
452 inline
operator eTarma_rng::randn453 operator eT () const
454 {
455 #if defined(ARMA_RNG_ALT)
456 {
457 return eT( arma_rng_alt::randn_val() );
458 }
459 #elif defined(ARMA_USE_EXTERN_RNG)
460 {
461 std::normal_distribution<double> local_n_distr;
462
463 return eT( local_n_distr(mt19937_64_instance) );
464 }
465 #else
466 {
467 return eT( arma_rng_cxx98::randn_val() );
468 }
469 #endif
470 }
471
472
473 inline
474 static
475 void
dual_valarma_rng::randn476 dual_val(eT& out1, eT& out2)
477 {
478 #if defined(ARMA_RNG_ALT)
479 {
480 arma_rng_alt::randn_dual_val(out1, out2);
481 }
482 #elif defined(ARMA_USE_EXTERN_RNG)
483 {
484 std::normal_distribution<double> local_n_distr;
485
486 out1 = eT( local_n_distr(mt19937_64_instance) );
487 out2 = eT( local_n_distr(mt19937_64_instance) );
488 }
489 #else
490 {
491 arma_rng_cxx98::randn_dual_val(out1, out2);
492 }
493 #endif
494 }
495
496
497 inline
498 static
499 void
fill_simplearma_rng::randn500 fill_simple(eT* mem, const uword N)
501 {
502 #if defined(ARMA_RNG_ALT)
503 {
504 // NOTE: old method to avoid regressions in user code that assumes specific sequence
505
506 uword i, j;
507
508 for(i=0, j=1; j < N; i+=2, j+=2) { arma_rng_alt::randn_dual_val( mem[i], mem[j] ); }
509
510 if(i < N) { mem[i] = eT( arma_rng_alt::randn_val() ); }
511 }
512 #elif defined(ARMA_USE_EXTERN_RNG)
513 {
514 std::normal_distribution<double> local_n_distr;
515
516 for(uword i=0; i < N; ++i) { mem[i] = eT( local_n_distr(mt19937_64_instance) ); }
517 }
518 #else
519 {
520 if(N == uword(1)) { mem[0] = eT( arma_rng_cxx98::randn_val() ); return; }
521
522 typedef typename std::mt19937_64::result_type local_seed_type;
523
524 std::mt19937_64 local_engine;
525 std::normal_distribution<double> local_n_distr;
526
527 local_engine.seed( local_seed_type(std::rand()) );
528
529 for(uword i=0; i < N; ++i) { mem[i] = eT( local_n_distr(local_engine) ); }
530 }
531 #endif
532 }
533
534
535 inline
536 static
537 void
fillarma_rng::randn538 fill(eT* mem, const uword N)
539 {
540 #if defined(ARMA_USE_OPENMP)
541 {
542 if((N < 1024) || omp_in_parallel()) { arma_rng::randn<eT>::fill_simple(mem, N); return; }
543
544 typedef typename std::mt19937_64::result_type local_seed_type;
545
546 const uword n_threads = uword( mp_thread_limit::get() );
547
548 std::vector< std::mt19937_64 > engine(n_threads);
549 std::vector< std::normal_distribution<double> > distr(n_threads);
550
551 for(uword t=0; t < n_threads; ++t)
552 {
553 std::mt19937_64& t_engine = engine[t];
554
555 t_engine.seed( local_seed_type(t) + local_seed_type(arma_rng::randi<local_seed_type>()) );
556 }
557
558 const uword chunk_size = N / n_threads;
559
560 #pragma omp parallel for schedule(static) num_threads(int(n_threads))
561 for(uword t=0; t < n_threads; ++t)
562 {
563 const uword start = (t+0) * chunk_size;
564 const uword endp1 = (t+1) * chunk_size;
565
566 std::mt19937_64& t_engine = engine[t];
567 std::normal_distribution<double>& t_distr = distr[t];
568
569 for(uword i=start; i < endp1; ++i) { mem[i] = eT( t_distr(t_engine)); }
570 }
571
572 std::mt19937_64& t0_engine = engine[0];
573 std::normal_distribution<double>& t0_distr = distr[0];
574
575 for(uword i=(n_threads*chunk_size); i < N; ++i) { mem[i] = eT( t0_distr(t0_engine)); }
576 }
577 #else
578 {
579 arma_rng::randn<eT>::fill_simple(mem, N);
580 }
581 #endif
582 }
583
584 };
585
586
587
588 template<typename T>
589 struct arma_rng::randn< std::complex<T> >
590 {
591 inline
operator std::complex<T>arma_rng::randn592 operator std::complex<T> () const
593 {
594 #if defined(_MSC_VER)
595 // attempt at workaround for MSVC bug
596 // does MS even test their so-called compilers before release?
597 T a;
598 T b;
599 #else
600 T a(0);
601 T b(0);
602 #endif
603
604 arma_rng::randn<T>::dual_val(a, b);
605
606 return std::complex<T>(a, b);
607 }
608
609
610 inline
611 static
612 void
dual_valarma_rng::randn613 dual_val(std::complex<T>& out1, std::complex<T>& out2)
614 {
615 #if defined(_MSC_VER)
616 T a;
617 T b;
618 #else
619 T a(0);
620 T b(0);
621 #endif
622
623 arma_rng::randn<T>::dual_val(a,b);
624 out1 = std::complex<T>(a,b);
625
626 arma_rng::randn<T>::dual_val(a,b);
627 out2 = std::complex<T>(a,b);
628 }
629
630
631 inline
632 static
633 void
fill_simplearma_rng::randn634 fill_simple(std::complex<T>* mem, const uword N)
635 {
636 #if defined(ARMA_RNG_ALT)
637 {
638 for(uword i=0; i < N; ++i) { mem[i] = std::complex<T>( arma_rng::randn< std::complex<T> >() ); }
639 }
640 #elif defined(ARMA_USE_EXTERN_RNG)
641 {
642 std::normal_distribution<double> local_n_distr;
643
644 for(uword i=0; i < N; ++i)
645 {
646 const T a = T( local_n_distr(mt19937_64_instance) );
647 const T b = T( local_n_distr(mt19937_64_instance) );
648
649 mem[i] = std::complex<T>(a,b);
650 }
651 }
652 #else
653 {
654 if(N == uword(1))
655 {
656 T a = T(0);
657 T b = T(0);
658
659 arma_rng_cxx98::randn_dual_val(a,b);
660
661 mem[0] = std::complex<T>(a,b);
662
663 return;
664 }
665
666 typedef typename std::mt19937_64::result_type local_seed_type;
667
668 std::mt19937_64 local_engine;
669 std::normal_distribution<double> local_n_distr;
670
671 local_engine.seed( local_seed_type(std::rand()) );
672
673 for(uword i=0; i < N; ++i)
674 {
675 const T a = T( local_n_distr(local_engine) );
676 const T b = T( local_n_distr(local_engine) );
677
678 mem[i] = std::complex<T>(a,b);
679 }
680 }
681 #endif
682 }
683
684
685 inline
686 static
687 void
fillarma_rng::randn688 fill(std::complex<T>* mem, const uword N)
689 {
690 #if defined(ARMA_USE_OPENMP)
691 {
692 if((N < 512) || omp_in_parallel()) { arma_rng::randn< std::complex<T> >::fill_simple(mem, N); return; }
693
694 typedef typename std::mt19937_64::result_type local_seed_type;
695
696 const uword n_threads = uword( mp_thread_limit::get() );
697
698 std::vector< std::mt19937_64 > engine(n_threads);
699 std::vector< std::normal_distribution<double> > distr(n_threads);
700
701 for(uword t=0; t < n_threads; ++t)
702 {
703 std::mt19937_64& t_engine = engine[t];
704
705 t_engine.seed( local_seed_type(t) + local_seed_type(arma_rng::randi<local_seed_type>()) );
706 }
707
708 const uword chunk_size = N / n_threads;
709
710 #pragma omp parallel for schedule(static) num_threads(int(n_threads))
711 for(uword t=0; t < n_threads; ++t)
712 {
713 const uword start = (t+0) * chunk_size;
714 const uword endp1 = (t+1) * chunk_size;
715
716 std::mt19937_64& t_engine = engine[t];
717 std::normal_distribution<double>& t_distr = distr[t];
718
719 for(uword i=start; i < endp1; ++i)
720 {
721 const T val1 = T( t_distr(t_engine) );
722 const T val2 = T( t_distr(t_engine) );
723
724 mem[i] = std::complex<T>(val1, val2);
725 }
726 }
727
728 std::mt19937_64& t0_engine = engine[0];
729 std::normal_distribution<double>& t0_distr = distr[0];
730
731 for(uword i=(n_threads*chunk_size); i < N; ++i)
732 {
733 const T val1 = T( t0_distr(t0_engine) );
734 const T val2 = T( t0_distr(t0_engine) );
735
736 mem[i] = std::complex<T>(val1, val2);
737 }
738 }
739 #else
740 {
741 arma_rng::randn< std::complex<T> >::fill_simple(mem, N);
742 }
743 #endif
744 }
745 };
746
747
748
749 //
750
751
752
753 template<typename eT>
754 struct arma_rng::randg
755 {
756 inline
757 static
758 void
fill_simplearma_rng::randg759 fill_simple(eT* mem, const uword N, const double a, const double b)
760 {
761 #if defined(ARMA_USE_EXTERN_RNG)
762 {
763 std::gamma_distribution<double> local_g_distr(a,b);
764
765 for(uword i=0; i<N; ++i) { mem[i] = eT(local_g_distr(mt19937_64_instance)); }
766 }
767 #else
768 {
769 typedef typename std::mt19937_64::result_type local_seed_type;
770
771 std::mt19937_64 local_engine;
772 std::gamma_distribution<double> local_g_distr(a,b);
773
774 local_engine.seed( local_seed_type(arma_rng::randi<local_seed_type>()) );
775
776 for(uword i=0; i<N; ++i) { mem[i] = eT(local_g_distr(local_engine)); }
777 }
778 #endif
779 }
780
781
782 inline
783 static
784 void
fillarma_rng::randg785 fill(eT* mem, const uword N, const double a, const double b)
786 {
787 #if defined(ARMA_USE_OPENMP)
788 {
789 if((N < 512) || omp_in_parallel()) { arma_rng::randg<eT>::fill_simple(mem, N, a, b); return; }
790
791 typedef std::mt19937_64 motor_type;
792 typedef std::mt19937_64::result_type ovum_type;
793 typedef std::gamma_distribution<double> distr_type;
794
795 const uword n_threads = uword( mp_thread_limit::get() );
796
797 std::vector<motor_type> g_motor(n_threads);
798 std::vector<distr_type> g_distr(n_threads);
799
800 const distr_type g_distr_base(a,b);
801
802 for(uword t=0; t < n_threads; ++t)
803 {
804 motor_type& g_motor_t = g_motor[t];
805 distr_type& g_distr_t = g_distr[t];
806
807 g_motor_t.seed( ovum_type(t) + ovum_type(arma_rng::randi<ovum_type>()) );
808
809 g_distr_t.param( g_distr_base.param() );
810 }
811
812 const uword chunk_size = N / n_threads;
813
814 #pragma omp parallel for schedule(static) num_threads(int(n_threads))
815 for(uword t=0; t < n_threads; ++t)
816 {
817 const uword start = (t+0) * chunk_size;
818 const uword endp1 = (t+1) * chunk_size;
819
820 motor_type& g_motor_t = g_motor[t];
821 distr_type& g_distr_t = g_distr[t];
822
823 for(uword i=start; i < endp1; ++i) { mem[i] = eT( g_distr_t(g_motor_t)); }
824 }
825
826 motor_type& g_motor_0 = g_motor[0];
827 distr_type& g_distr_0 = g_distr[0];
828
829 for(uword i=(n_threads*chunk_size); i < N; ++i) { mem[i] = eT( g_distr_0(g_motor_0)); }
830 }
831 #else
832 {
833 arma_rng::randg<eT>::fill_simple(mem, N, a, b);
834 }
835 #endif
836 }
837
838 };
839
840
841
842 //! @}
843