libstdc++
ext/random.tcc
Go to the documentation of this file.
1// Random number extensions -*- C++ -*-
2
3// Copyright (C) 2012-2023 Free Software Foundation, Inc.
4//
5// This file is part of the GNU ISO C++ Library. This library is free
6// software; you can redistribute it and/or modify it under the
7// terms of the GNU General Public License as published by the
8// Free Software Foundation; either version 3, or (at your option)
9// any later version.
10
11// This library is distributed in the hope that it will be useful,
12// but WITHOUT ANY WARRANTY; without even the implied warranty of
13// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14// GNU General Public License for more details.
15
16// Under Section 7 of GPL version 3, you are granted additional
17// permissions described in the GCC Runtime Library Exception, version
18// 3.1, as published by the Free Software Foundation.
19
20// You should have received a copy of the GNU General Public License and
21// a copy of the GCC Runtime Library Exception along with this program;
22// see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
23// <http://www.gnu.org/licenses/>.
24
25/** @file ext/random.tcc
26 * This is an internal header file, included by other library headers.
27 * Do not attempt to use it directly. @headername{ext/random}
28 */
29
30#ifndef _EXT_RANDOM_TCC
31#define _EXT_RANDOM_TCC 1
32
33#pragma GCC system_header
34
35#include <bits/requires_hosted.h> // GNU extensions are currently omitted
36
37namespace __gnu_cxx _GLIBCXX_VISIBILITY(default)
38{
39_GLIBCXX_BEGIN_NAMESPACE_VERSION
40
41#if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
42
43 template<typename _UIntType, size_t __m,
44 size_t __pos1, size_t __sl1, size_t __sl2,
45 size_t __sr1, size_t __sr2,
46 uint32_t __msk1, uint32_t __msk2,
47 uint32_t __msk3, uint32_t __msk4,
48 uint32_t __parity1, uint32_t __parity2,
49 uint32_t __parity3, uint32_t __parity4>
50 void simd_fast_mersenne_twister_engine<_UIntType, __m,
51 __pos1, __sl1, __sl2, __sr1, __sr2,
52 __msk1, __msk2, __msk3, __msk4,
53 __parity1, __parity2, __parity3,
54 __parity4>::
55 seed(_UIntType __seed)
56 {
57 _M_state32[0] = static_cast<uint32_t>(__seed);
58 for (size_t __i = 1; __i < _M_nstate32; ++__i)
59 _M_state32[__i] = (1812433253UL
60 * (_M_state32[__i - 1] ^ (_M_state32[__i - 1] >> 30))
61 + __i);
62 _M_pos = state_size;
63 _M_period_certification();
64 }
65
66
67 namespace {
68
69 inline uint32_t _Func1(uint32_t __x)
70 {
71 return (__x ^ (__x >> 27)) * UINT32_C(1664525);
72 }
73
74 inline uint32_t _Func2(uint32_t __x)
75 {
76 return (__x ^ (__x >> 27)) * UINT32_C(1566083941);
77 }
78
79 }
80
81
82 template<typename _UIntType, size_t __m,
83 size_t __pos1, size_t __sl1, size_t __sl2,
84 size_t __sr1, size_t __sr2,
85 uint32_t __msk1, uint32_t __msk2,
86 uint32_t __msk3, uint32_t __msk4,
87 uint32_t __parity1, uint32_t __parity2,
88 uint32_t __parity3, uint32_t __parity4>
89 template<typename _Sseq>
90 auto
91 simd_fast_mersenne_twister_engine<_UIntType, __m,
92 __pos1, __sl1, __sl2, __sr1, __sr2,
93 __msk1, __msk2, __msk3, __msk4,
94 __parity1, __parity2, __parity3,
95 __parity4>::
96 seed(_Sseq& __q)
97 -> _If_seed_seq<_Sseq>
98 {
99 size_t __lag;
100
101 if (_M_nstate32 >= 623)
102 __lag = 11;
103 else if (_M_nstate32 >= 68)
104 __lag = 7;
105 else if (_M_nstate32 >= 39)
106 __lag = 5;
107 else
108 __lag = 3;
109 const size_t __mid = (_M_nstate32 - __lag) / 2;
110
111 std::fill(_M_state32, _M_state32 + _M_nstate32, UINT32_C(0x8b8b8b8b));
112 uint32_t __arr[_M_nstate32];
113 __q.generate(__arr + 0, __arr + _M_nstate32);
114
115 uint32_t __r = _Func1(_M_state32[0] ^ _M_state32[__mid]
116 ^ _M_state32[_M_nstate32 - 1]);
117 _M_state32[__mid] += __r;
118 __r += _M_nstate32;
119 _M_state32[__mid + __lag] += __r;
120 _M_state32[0] = __r;
121
122 for (size_t __i = 1, __j = 0; __j < _M_nstate32; ++__j)
123 {
124 __r = _Func1(_M_state32[__i]
125 ^ _M_state32[(__i + __mid) % _M_nstate32]
126 ^ _M_state32[(__i + _M_nstate32 - 1) % _M_nstate32]);
127 _M_state32[(__i + __mid) % _M_nstate32] += __r;
128 __r += __arr[__j] + __i;
129 _M_state32[(__i + __mid + __lag) % _M_nstate32] += __r;
130 _M_state32[__i] = __r;
131 __i = (__i + 1) % _M_nstate32;
132 }
133 for (size_t __j = 0; __j < _M_nstate32; ++__j)
134 {
135 const size_t __i = (__j + 1) % _M_nstate32;
136 __r = _Func2(_M_state32[__i]
137 + _M_state32[(__i + __mid) % _M_nstate32]
138 + _M_state32[(__i + _M_nstate32 - 1) % _M_nstate32]);
139 _M_state32[(__i + __mid) % _M_nstate32] ^= __r;
140 __r -= __i;
141 _M_state32[(__i + __mid + __lag) % _M_nstate32] ^= __r;
142 _M_state32[__i] = __r;
143 }
144
145 _M_pos = state_size;
146 _M_period_certification();
147 }
148
149
150 template<typename _UIntType, size_t __m,
151 size_t __pos1, size_t __sl1, size_t __sl2,
152 size_t __sr1, size_t __sr2,
153 uint32_t __msk1, uint32_t __msk2,
154 uint32_t __msk3, uint32_t __msk4,
155 uint32_t __parity1, uint32_t __parity2,
156 uint32_t __parity3, uint32_t __parity4>
157 void simd_fast_mersenne_twister_engine<_UIntType, __m,
158 __pos1, __sl1, __sl2, __sr1, __sr2,
159 __msk1, __msk2, __msk3, __msk4,
160 __parity1, __parity2, __parity3,
161 __parity4>::
162 _M_period_certification(void)
163 {
164 static const uint32_t __parity[4] = { __parity1, __parity2,
165 __parity3, __parity4 };
166 uint32_t __inner = 0;
167 for (size_t __i = 0; __i < 4; ++__i)
168 if (__parity[__i] != 0)
169 __inner ^= _M_state32[__i] & __parity[__i];
170
171 if (__builtin_parity(__inner) & 1)
172 return;
173 for (size_t __i = 0; __i < 4; ++__i)
174 if (__parity[__i] != 0)
175 {
176 _M_state32[__i] ^= 1 << (__builtin_ffs(__parity[__i]) - 1);
177 return;
178 }
179 __builtin_unreachable();
180 }
181
182
183 template<typename _UIntType, size_t __m,
184 size_t __pos1, size_t __sl1, size_t __sl2,
185 size_t __sr1, size_t __sr2,
186 uint32_t __msk1, uint32_t __msk2,
187 uint32_t __msk3, uint32_t __msk4,
188 uint32_t __parity1, uint32_t __parity2,
189 uint32_t __parity3, uint32_t __parity4>
190 void simd_fast_mersenne_twister_engine<_UIntType, __m,
191 __pos1, __sl1, __sl2, __sr1, __sr2,
192 __msk1, __msk2, __msk3, __msk4,
193 __parity1, __parity2, __parity3,
194 __parity4>::
195 discard(unsigned long long __z)
196 {
197 while (__z > state_size - _M_pos)
198 {
199 __z -= state_size - _M_pos;
200
201 _M_gen_rand();
202 }
203
204 _M_pos += __z;
205 }
206
207
208#ifndef _GLIBCXX_OPT_HAVE_RANDOM_SFMT_GEN_READ
209
210 namespace {
211
212 template<size_t __shift>
213 inline void __rshift(uint32_t *__out, const uint32_t *__in)
214 {
215 uint64_t __th = ((static_cast<uint64_t>(__in[3]) << 32)
216 | static_cast<uint64_t>(__in[2]));
217 uint64_t __tl = ((static_cast<uint64_t>(__in[1]) << 32)
218 | static_cast<uint64_t>(__in[0]));
219
220 uint64_t __oh = __th >> (__shift * 8);
221 uint64_t __ol = __tl >> (__shift * 8);
222 __ol |= __th << (64 - __shift * 8);
223 __out[1] = static_cast<uint32_t>(__ol >> 32);
224 __out[0] = static_cast<uint32_t>(__ol);
225 __out[3] = static_cast<uint32_t>(__oh >> 32);
226 __out[2] = static_cast<uint32_t>(__oh);
227 }
228
229
230 template<size_t __shift>
231 inline void __lshift(uint32_t *__out, const uint32_t *__in)
232 {
233 uint64_t __th = ((static_cast<uint64_t>(__in[3]) << 32)
234 | static_cast<uint64_t>(__in[2]));
235 uint64_t __tl = ((static_cast<uint64_t>(__in[1]) << 32)
236 | static_cast<uint64_t>(__in[0]));
237
238 uint64_t __oh = __th << (__shift * 8);
239 uint64_t __ol = __tl << (__shift * 8);
240 __oh |= __tl >> (64 - __shift * 8);
241 __out[1] = static_cast<uint32_t>(__ol >> 32);
242 __out[0] = static_cast<uint32_t>(__ol);
243 __out[3] = static_cast<uint32_t>(__oh >> 32);
244 __out[2] = static_cast<uint32_t>(__oh);
245 }
246
247
248 template<size_t __sl1, size_t __sl2, size_t __sr1, size_t __sr2,
249 uint32_t __msk1, uint32_t __msk2, uint32_t __msk3, uint32_t __msk4>
250 inline void __recursion(uint32_t *__r,
251 const uint32_t *__a, const uint32_t *__b,
252 const uint32_t *__c, const uint32_t *__d)
253 {
254 uint32_t __x[4];
255 uint32_t __y[4];
256
257 __lshift<__sl2>(__x, __a);
258 __rshift<__sr2>(__y, __c);
259 __r[0] = (__a[0] ^ __x[0] ^ ((__b[0] >> __sr1) & __msk1)
260 ^ __y[0] ^ (__d[0] << __sl1));
261 __r[1] = (__a[1] ^ __x[1] ^ ((__b[1] >> __sr1) & __msk2)
262 ^ __y[1] ^ (__d[1] << __sl1));
263 __r[2] = (__a[2] ^ __x[2] ^ ((__b[2] >> __sr1) & __msk3)
264 ^ __y[2] ^ (__d[2] << __sl1));
265 __r[3] = (__a[3] ^ __x[3] ^ ((__b[3] >> __sr1) & __msk4)
266 ^ __y[3] ^ (__d[3] << __sl1));
267 }
268
269 }
270
271
272 template<typename _UIntType, size_t __m,
273 size_t __pos1, size_t __sl1, size_t __sl2,
274 size_t __sr1, size_t __sr2,
275 uint32_t __msk1, uint32_t __msk2,
276 uint32_t __msk3, uint32_t __msk4,
277 uint32_t __parity1, uint32_t __parity2,
278 uint32_t __parity3, uint32_t __parity4>
279 void simd_fast_mersenne_twister_engine<_UIntType, __m,
280 __pos1, __sl1, __sl2, __sr1, __sr2,
281 __msk1, __msk2, __msk3, __msk4,
282 __parity1, __parity2, __parity3,
283 __parity4>::
284 _M_gen_rand(void)
285 {
286 const uint32_t *__r1 = &_M_state32[_M_nstate32 - 8];
287 const uint32_t *__r2 = &_M_state32[_M_nstate32 - 4];
288 static constexpr size_t __pos1_32 = __pos1 * 4;
289
290 size_t __i;
291 for (__i = 0; __i < _M_nstate32 - __pos1_32; __i += 4)
292 {
293 __recursion<__sl1, __sl2, __sr1, __sr2,
294 __msk1, __msk2, __msk3, __msk4>
295 (&_M_state32[__i], &_M_state32[__i],
296 &_M_state32[__i + __pos1_32], __r1, __r2);
297 __r1 = __r2;
298 __r2 = &_M_state32[__i];
299 }
300
301 for (; __i < _M_nstate32; __i += 4)
302 {
303 __recursion<__sl1, __sl2, __sr1, __sr2,
304 __msk1, __msk2, __msk3, __msk4>
305 (&_M_state32[__i], &_M_state32[__i],
306 &_M_state32[__i + __pos1_32 - _M_nstate32], __r1, __r2);
307 __r1 = __r2;
308 __r2 = &_M_state32[__i];
309 }
310
311 _M_pos = 0;
312 }
313
314#endif
315
316#ifndef _GLIBCXX_OPT_HAVE_RANDOM_SFMT_OPERATOREQUAL
317 template<typename _UIntType, size_t __m,
318 size_t __pos1, size_t __sl1, size_t __sl2,
319 size_t __sr1, size_t __sr2,
320 uint32_t __msk1, uint32_t __msk2,
321 uint32_t __msk3, uint32_t __msk4,
322 uint32_t __parity1, uint32_t __parity2,
323 uint32_t __parity3, uint32_t __parity4>
324 bool
325 operator==(const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
326 __m, __pos1, __sl1, __sl2, __sr1, __sr2,
327 __msk1, __msk2, __msk3, __msk4,
328 __parity1, __parity2, __parity3, __parity4>& __lhs,
329 const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
330 __m, __pos1, __sl1, __sl2, __sr1, __sr2,
331 __msk1, __msk2, __msk3, __msk4,
332 __parity1, __parity2, __parity3, __parity4>& __rhs)
333 {
334 typedef __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
335 __m, __pos1, __sl1, __sl2, __sr1, __sr2,
336 __msk1, __msk2, __msk3, __msk4,
337 __parity1, __parity2, __parity3, __parity4> __engine;
338 return (std::equal(__lhs._M_stateT,
339 __lhs._M_stateT + __engine::state_size,
340 __rhs._M_stateT)
341 && __lhs._M_pos == __rhs._M_pos);
342 }
343#endif
344
345 template<typename _UIntType, size_t __m,
346 size_t __pos1, size_t __sl1, size_t __sl2,
347 size_t __sr1, size_t __sr2,
348 uint32_t __msk1, uint32_t __msk2,
349 uint32_t __msk3, uint32_t __msk4,
350 uint32_t __parity1, uint32_t __parity2,
351 uint32_t __parity3, uint32_t __parity4,
352 typename _CharT, typename _Traits>
355 const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
356 __m, __pos1, __sl1, __sl2, __sr1, __sr2,
357 __msk1, __msk2, __msk3, __msk4,
358 __parity1, __parity2, __parity3, __parity4>& __x)
359 {
360 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
361 typedef typename __ostream_type::ios_base __ios_base;
362
363 const typename __ios_base::fmtflags __flags = __os.flags();
364 const _CharT __fill = __os.fill();
365 const _CharT __space = __os.widen(' ');
366 __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
367 __os.fill(__space);
368
369 for (size_t __i = 0; __i < __x._M_nstate32; ++__i)
370 __os << __x._M_state32[__i] << __space;
371 __os << __x._M_pos;
372
373 __os.flags(__flags);
374 __os.fill(__fill);
375 return __os;
376 }
377
378
379 template<typename _UIntType, size_t __m,
380 size_t __pos1, size_t __sl1, size_t __sl2,
381 size_t __sr1, size_t __sr2,
382 uint32_t __msk1, uint32_t __msk2,
383 uint32_t __msk3, uint32_t __msk4,
384 uint32_t __parity1, uint32_t __parity2,
385 uint32_t __parity3, uint32_t __parity4,
386 typename _CharT, typename _Traits>
389 __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
390 __m, __pos1, __sl1, __sl2, __sr1, __sr2,
391 __msk1, __msk2, __msk3, __msk4,
392 __parity1, __parity2, __parity3, __parity4>& __x)
393 {
394 typedef std::basic_istream<_CharT, _Traits> __istream_type;
395 typedef typename __istream_type::ios_base __ios_base;
396
397 const typename __ios_base::fmtflags __flags = __is.flags();
398 __is.flags(__ios_base::dec | __ios_base::skipws);
399
400 for (size_t __i = 0; __i < __x._M_nstate32; ++__i)
401 __is >> __x._M_state32[__i];
402 __is >> __x._M_pos;
403
404 __is.flags(__flags);
405 return __is;
406 }
407
408#endif // __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
409
410 /**
411 * Iteration method due to M.D. J<o:>hnk.
412 *
413 * M.D. J<o:>hnk, Erzeugung von betaverteilten und gammaverteilten
414 * Zufallszahlen, Metrika, Volume 8, 1964
415 */
416 template<typename _RealType>
417 template<typename _UniformRandomNumberGenerator>
418 typename beta_distribution<_RealType>::result_type
419 beta_distribution<_RealType>::
420 operator()(_UniformRandomNumberGenerator& __urng,
421 const param_type& __param)
422 {
423 std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
424 __aurng(__urng);
425
426 result_type __x, __y;
427 do
428 {
429 __x = std::exp(std::log(__aurng()) / __param.alpha());
430 __y = std::exp(std::log(__aurng()) / __param.beta());
431 }
432 while (__x + __y > result_type(1));
433
434 return __x / (__x + __y);
435 }
436
437 template<typename _RealType>
438 template<typename _OutputIterator,
439 typename _UniformRandomNumberGenerator>
440 void
441 beta_distribution<_RealType>::
442 __generate_impl(_OutputIterator __f, _OutputIterator __t,
443 _UniformRandomNumberGenerator& __urng,
444 const param_type& __param)
445 {
446 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
447 result_type>)
448
449 std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
450 __aurng(__urng);
451
452 while (__f != __t)
453 {
454 result_type __x, __y;
455 do
456 {
457 __x = std::exp(std::log(__aurng()) / __param.alpha());
458 __y = std::exp(std::log(__aurng()) / __param.beta());
459 }
460 while (__x + __y > result_type(1));
461
462 *__f++ = __x / (__x + __y);
463 }
464 }
465
466 template<typename _RealType, typename _CharT, typename _Traits>
469 const __gnu_cxx::beta_distribution<_RealType>& __x)
470 {
471 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
472 typedef typename __ostream_type::ios_base __ios_base;
473
474 const typename __ios_base::fmtflags __flags = __os.flags();
475 const _CharT __fill = __os.fill();
476 const std::streamsize __precision = __os.precision();
477 const _CharT __space = __os.widen(' ');
478 __os.flags(__ios_base::scientific | __ios_base::left);
479 __os.fill(__space);
481
482 __os << __x.alpha() << __space << __x.beta();
483
484 __os.flags(__flags);
485 __os.fill(__fill);
486 __os.precision(__precision);
487 return __os;
488 }
489
490 template<typename _RealType, typename _CharT, typename _Traits>
493 __gnu_cxx::beta_distribution<_RealType>& __x)
494 {
495 typedef std::basic_istream<_CharT, _Traits> __istream_type;
496 typedef typename __istream_type::ios_base __ios_base;
497
498 const typename __ios_base::fmtflags __flags = __is.flags();
499 __is.flags(__ios_base::dec | __ios_base::skipws);
500
501 _RealType __alpha_val, __beta_val;
502 __is >> __alpha_val >> __beta_val;
503 __x.param(typename __gnu_cxx::beta_distribution<_RealType>::
504 param_type(__alpha_val, __beta_val));
505
506 __is.flags(__flags);
507 return __is;
508 }
509
510
511 template<std::size_t _Dimen, typename _RealType>
512 template<typename _InputIterator1, typename _InputIterator2>
513 void
514 normal_mv_distribution<_Dimen, _RealType>::param_type::
515 _M_init_full(_InputIterator1 __meanbegin, _InputIterator1 __meanend,
516 _InputIterator2 __varcovbegin, _InputIterator2 __varcovend)
517 {
518 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>)
519 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>)
520 std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()),
521 _M_mean.end(), _RealType(0));
522
523 // Perform the Cholesky decomposition
524 auto __w = _M_t.begin();
525 for (size_t __j = 0; __j < _Dimen; ++__j)
526 {
527 _RealType __sum = _RealType(0);
528
529 auto __slitbegin = __w;
530 auto __cit = _M_t.begin();
531 for (size_t __i = 0; __i < __j; ++__i)
532 {
533 auto __slit = __slitbegin;
534 _RealType __s = *__varcovbegin++;
535 for (size_t __k = 0; __k < __i; ++__k)
536 __s -= *__slit++ * *__cit++;
537
538 *__w++ = __s /= *__cit++;
539 __sum += __s * __s;
540 }
541
542 __sum = *__varcovbegin - __sum;
543 if (__builtin_expect(__sum <= _RealType(0), 0))
544 std::__throw_runtime_error(__N("normal_mv_distribution::"
545 "param_type::_M_init_full"));
546 *__w++ = std::sqrt(__sum);
547
548 std::advance(__varcovbegin, _Dimen - __j);
549 }
550 }
551
552 template<std::size_t _Dimen, typename _RealType>
553 template<typename _InputIterator1, typename _InputIterator2>
554 void
555 normal_mv_distribution<_Dimen, _RealType>::param_type::
556 _M_init_lower(_InputIterator1 __meanbegin, _InputIterator1 __meanend,
557 _InputIterator2 __varcovbegin, _InputIterator2 __varcovend)
558 {
559 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>)
560 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>)
561 std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()),
562 _M_mean.end(), _RealType(0));
563
564 // Perform the Cholesky decomposition
565 auto __w = _M_t.begin();
566 for (size_t __j = 0; __j < _Dimen; ++__j)
567 {
568 _RealType __sum = _RealType(0);
569
570 auto __slitbegin = __w;
571 auto __cit = _M_t.begin();
572 for (size_t __i = 0; __i < __j; ++__i)
573 {
574 auto __slit = __slitbegin;
575 _RealType __s = *__varcovbegin++;
576 for (size_t __k = 0; __k < __i; ++__k)
577 __s -= *__slit++ * *__cit++;
578
579 *__w++ = __s /= *__cit++;
580 __sum += __s * __s;
581 }
582
583 __sum = *__varcovbegin++ - __sum;
584 if (__builtin_expect(__sum <= _RealType(0), 0))
585 std::__throw_runtime_error(__N("normal_mv_distribution::"
586 "param_type::_M_init_lower"));
587 *__w++ = std::sqrt(__sum);
588 }
589 }
590
591 template<std::size_t _Dimen, typename _RealType>
592 template<typename _InputIterator1, typename _InputIterator2>
593 void
594 normal_mv_distribution<_Dimen, _RealType>::param_type::
595 _M_init_diagonal(_InputIterator1 __meanbegin, _InputIterator1 __meanend,
596 _InputIterator2 __varbegin, _InputIterator2 __varend)
597 {
598 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator1>)
599 __glibcxx_function_requires(_InputIteratorConcept<_InputIterator2>)
600 std::fill(std::copy(__meanbegin, __meanend, _M_mean.begin()),
601 _M_mean.end(), _RealType(0));
602
603 auto __w = _M_t.begin();
604 size_t __step = 0;
605 while (__varbegin != __varend)
606 {
607 std::fill_n(__w, __step, _RealType(0));
608 __w += __step++;
609 if (__builtin_expect(*__varbegin < _RealType(0), 0))
610 std::__throw_runtime_error(__N("normal_mv_distribution::"
611 "param_type::_M_init_diagonal"));
612 *__w++ = std::sqrt(*__varbegin++);
613 }
614 }
615
616 template<std::size_t _Dimen, typename _RealType>
617 template<typename _UniformRandomNumberGenerator>
618 typename normal_mv_distribution<_Dimen, _RealType>::result_type
619 normal_mv_distribution<_Dimen, _RealType>::
620 operator()(_UniformRandomNumberGenerator& __urng,
621 const param_type& __param)
622 {
623 result_type __ret;
624
625 _M_nd.__generate(__ret.begin(), __ret.end(), __urng);
626
627 auto __t_it = __param._M_t.crbegin();
628 for (size_t __i = _Dimen; __i > 0; --__i)
629 {
630 _RealType __sum = _RealType(0);
631 for (size_t __j = __i; __j > 0; --__j)
632 __sum += __ret[__j - 1] * *__t_it++;
633 __ret[__i - 1] = __sum;
634 }
635
636 return __ret;
637 }
638
639 template<std::size_t _Dimen, typename _RealType>
640 template<typename _ForwardIterator, typename _UniformRandomNumberGenerator>
641 void
642 normal_mv_distribution<_Dimen, _RealType>::
643 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
644 _UniformRandomNumberGenerator& __urng,
645 const param_type& __param)
646 {
647 __glibcxx_function_requires(_Mutable_ForwardIteratorConcept<
648 _ForwardIterator>)
649 while (__f != __t)
650 *__f++ = this->operator()(__urng, __param);
651 }
652
653 template<size_t _Dimen, typename _RealType>
654 bool
655 operator==(const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
656 __d1,
657 const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
658 __d2)
659 {
660 return __d1._M_param == __d2._M_param && __d1._M_nd == __d2._M_nd;
661 }
662
663 template<size_t _Dimen, typename _RealType, typename _CharT, typename _Traits>
666 const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& __x)
667 {
668 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
669 typedef typename __ostream_type::ios_base __ios_base;
670
671 const typename __ios_base::fmtflags __flags = __os.flags();
672 const _CharT __fill = __os.fill();
673 const std::streamsize __precision = __os.precision();
674 const _CharT __space = __os.widen(' ');
675 __os.flags(__ios_base::scientific | __ios_base::left);
676 __os.fill(__space);
678
679 auto __mean = __x._M_param.mean();
680 for (auto __it : __mean)
681 __os << __it << __space;
682 auto __t = __x._M_param.varcov();
683 for (auto __it : __t)
684 __os << __it << __space;
685
686 __os << __x._M_nd;
687
688 __os.flags(__flags);
689 __os.fill(__fill);
690 __os.precision(__precision);
691 return __os;
692 }
693
694 template<size_t _Dimen, typename _RealType, typename _CharT, typename _Traits>
697 __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>& __x)
698 {
699 typedef std::basic_istream<_CharT, _Traits> __istream_type;
700 typedef typename __istream_type::ios_base __ios_base;
701
702 const typename __ios_base::fmtflags __flags = __is.flags();
703 __is.flags(__ios_base::dec | __ios_base::skipws);
704
706 for (auto& __it : __mean)
707 __is >> __it;
708 std::array<_RealType, _Dimen * (_Dimen + 1) / 2> __varcov;
709 for (auto& __it : __varcov)
710 __is >> __it;
711
712 __is >> __x._M_nd;
713
714 // The param_type temporary is built with a private constructor,
715 // to skip the Cholesky decomposition that would be performed
716 // otherwise.
717 __x.param(typename normal_mv_distribution<_Dimen, _RealType>::
718 param_type(__mean, __varcov));
719
720 __is.flags(__flags);
721 return __is;
722 }
723
724
725 template<typename _RealType>
726 template<typename _OutputIterator,
727 typename _UniformRandomNumberGenerator>
728 void
729 rice_distribution<_RealType>::
730 __generate_impl(_OutputIterator __f, _OutputIterator __t,
731 _UniformRandomNumberGenerator& __urng,
732 const param_type& __p)
733 {
734 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
735 result_type>)
736
737 while (__f != __t)
738 {
740 __px(__p.nu(), __p.sigma()), __py(result_type(0), __p.sigma());
741 result_type __x = this->_M_ndx(__px, __urng);
742 result_type __y = this->_M_ndy(__py, __urng);
743#if _GLIBCXX_USE_C99_MATH_TR1
744 *__f++ = std::hypot(__x, __y);
745#else
746 *__f++ = std::sqrt(__x * __x + __y * __y);
747#endif
748 }
749 }
750
751 template<typename _RealType, typename _CharT, typename _Traits>
754 const rice_distribution<_RealType>& __x)
755 {
756 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
757 typedef typename __ostream_type::ios_base __ios_base;
758
759 const typename __ios_base::fmtflags __flags = __os.flags();
760 const _CharT __fill = __os.fill();
761 const std::streamsize __precision = __os.precision();
762 const _CharT __space = __os.widen(' ');
763 __os.flags(__ios_base::scientific | __ios_base::left);
764 __os.fill(__space);
766
767 __os << __x.nu() << __space << __x.sigma();
768 __os << __space << __x._M_ndx;
769 __os << __space << __x._M_ndy;
770
771 __os.flags(__flags);
772 __os.fill(__fill);
773 __os.precision(__precision);
774 return __os;
775 }
776
777 template<typename _RealType, typename _CharT, typename _Traits>
780 rice_distribution<_RealType>& __x)
781 {
782 typedef std::basic_istream<_CharT, _Traits> __istream_type;
783 typedef typename __istream_type::ios_base __ios_base;
784
785 const typename __ios_base::fmtflags __flags = __is.flags();
786 __is.flags(__ios_base::dec | __ios_base::skipws);
787
788 _RealType __nu_val, __sigma_val;
789 __is >> __nu_val >> __sigma_val;
790 __is >> __x._M_ndx;
791 __is >> __x._M_ndy;
792 __x.param(typename rice_distribution<_RealType>::
793 param_type(__nu_val, __sigma_val));
794
795 __is.flags(__flags);
796 return __is;
797 }
798
799
800 template<typename _RealType>
801 template<typename _OutputIterator,
802 typename _UniformRandomNumberGenerator>
803 void
804 nakagami_distribution<_RealType>::
805 __generate_impl(_OutputIterator __f, _OutputIterator __t,
806 _UniformRandomNumberGenerator& __urng,
807 const param_type& __p)
808 {
809 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
810 result_type>)
811
812 typename std::gamma_distribution<result_type>::param_type
813 __pg(__p.mu(), __p.omega() / __p.mu());
814 while (__f != __t)
815 *__f++ = std::sqrt(this->_M_gd(__pg, __urng));
816 }
817
818 template<typename _RealType, typename _CharT, typename _Traits>
819 std::basic_ostream<_CharT, _Traits>&
820 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
821 const nakagami_distribution<_RealType>& __x)
822 {
823 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
824 typedef typename __ostream_type::ios_base __ios_base;
825
826 const typename __ios_base::fmtflags __flags = __os.flags();
827 const _CharT __fill = __os.fill();
828 const std::streamsize __precision = __os.precision();
829 const _CharT __space = __os.widen(' ');
830 __os.flags(__ios_base::scientific | __ios_base::left);
831 __os.fill(__space);
833
834 __os << __x.mu() << __space << __x.omega();
835 __os << __space << __x._M_gd;
836
837 __os.flags(__flags);
838 __os.fill(__fill);
839 __os.precision(__precision);
840 return __os;
841 }
842
843 template<typename _RealType, typename _CharT, typename _Traits>
846 nakagami_distribution<_RealType>& __x)
847 {
848 typedef std::basic_istream<_CharT, _Traits> __istream_type;
849 typedef typename __istream_type::ios_base __ios_base;
850
851 const typename __ios_base::fmtflags __flags = __is.flags();
852 __is.flags(__ios_base::dec | __ios_base::skipws);
853
854 _RealType __mu_val, __omega_val;
855 __is >> __mu_val >> __omega_val;
856 __is >> __x._M_gd;
857 __x.param(typename nakagami_distribution<_RealType>::
858 param_type(__mu_val, __omega_val));
859
860 __is.flags(__flags);
861 return __is;
862 }
863
864
865 template<typename _RealType>
866 template<typename _OutputIterator,
867 typename _UniformRandomNumberGenerator>
868 void
869 pareto_distribution<_RealType>::
870 __generate_impl(_OutputIterator __f, _OutputIterator __t,
871 _UniformRandomNumberGenerator& __urng,
872 const param_type& __p)
873 {
874 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
875 result_type>)
876
877 result_type __mu_val = __p.mu();
878 result_type __malphinv = -result_type(1) / __p.alpha();
879 while (__f != __t)
880 *__f++ = __mu_val * std::pow(this->_M_ud(__urng), __malphinv);
881 }
882
883 template<typename _RealType, typename _CharT, typename _Traits>
884 std::basic_ostream<_CharT, _Traits>&
885 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
886 const pareto_distribution<_RealType>& __x)
887 {
888 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
889 typedef typename __ostream_type::ios_base __ios_base;
890
891 const typename __ios_base::fmtflags __flags = __os.flags();
892 const _CharT __fill = __os.fill();
893 const std::streamsize __precision = __os.precision();
894 const _CharT __space = __os.widen(' ');
895 __os.flags(__ios_base::scientific | __ios_base::left);
896 __os.fill(__space);
898
899 __os << __x.alpha() << __space << __x.mu();
900 __os << __space << __x._M_ud;
901
902 __os.flags(__flags);
903 __os.fill(__fill);
904 __os.precision(__precision);
905 return __os;
906 }
907
908 template<typename _RealType, typename _CharT, typename _Traits>
911 pareto_distribution<_RealType>& __x)
912 {
913 typedef std::basic_istream<_CharT, _Traits> __istream_type;
914 typedef typename __istream_type::ios_base __ios_base;
915
916 const typename __ios_base::fmtflags __flags = __is.flags();
917 __is.flags(__ios_base::dec | __ios_base::skipws);
918
919 _RealType __alpha_val, __mu_val;
920 __is >> __alpha_val >> __mu_val;
921 __is >> __x._M_ud;
922 __x.param(typename pareto_distribution<_RealType>::
923 param_type(__alpha_val, __mu_val));
924
925 __is.flags(__flags);
926 return __is;
927 }
928
929
930 template<typename _RealType>
931 template<typename _UniformRandomNumberGenerator>
932 typename k_distribution<_RealType>::result_type
933 k_distribution<_RealType>::
934 operator()(_UniformRandomNumberGenerator& __urng)
935 {
936 result_type __x = this->_M_gd1(__urng);
937 result_type __y = this->_M_gd2(__urng);
938 return std::sqrt(__x * __y);
939 }
940
941 template<typename _RealType>
942 template<typename _UniformRandomNumberGenerator>
943 typename k_distribution<_RealType>::result_type
944 k_distribution<_RealType>::
945 operator()(_UniformRandomNumberGenerator& __urng,
946 const param_type& __p)
947 {
949 __p1(__p.lambda(), result_type(1) / __p.lambda()),
950 __p2(__p.nu(), __p.mu() / __p.nu());
951 result_type __x = this->_M_gd1(__p1, __urng);
952 result_type __y = this->_M_gd2(__p2, __urng);
953 return std::sqrt(__x * __y);
954 }
955
956 template<typename _RealType>
957 template<typename _OutputIterator,
958 typename _UniformRandomNumberGenerator>
959 void
960 k_distribution<_RealType>::
961 __generate_impl(_OutputIterator __f, _OutputIterator __t,
962 _UniformRandomNumberGenerator& __urng,
963 const param_type& __p)
964 {
965 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
966 result_type>)
967
968 typename std::gamma_distribution<result_type>::param_type
969 __p1(__p.lambda(), result_type(1) / __p.lambda()),
970 __p2(__p.nu(), __p.mu() / __p.nu());
971 while (__f != __t)
972 {
973 result_type __x = this->_M_gd1(__p1, __urng);
974 result_type __y = this->_M_gd2(__p2, __urng);
975 *__f++ = std::sqrt(__x * __y);
976 }
977 }
978
979 template<typename _RealType, typename _CharT, typename _Traits>
982 const k_distribution<_RealType>& __x)
983 {
984 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
985 typedef typename __ostream_type::ios_base __ios_base;
986
987 const typename __ios_base::fmtflags __flags = __os.flags();
988 const _CharT __fill = __os.fill();
989 const std::streamsize __precision = __os.precision();
990 const _CharT __space = __os.widen(' ');
991 __os.flags(__ios_base::scientific | __ios_base::left);
992 __os.fill(__space);
994
995 __os << __x.lambda() << __space << __x.mu() << __space << __x.nu();
996 __os << __space << __x._M_gd1;
997 __os << __space << __x._M_gd2;
998
999 __os.flags(__flags);
1000 __os.fill(__fill);
1001 __os.precision(__precision);
1002 return __os;
1003 }
1004
1005 template<typename _RealType, typename _CharT, typename _Traits>
1008 k_distribution<_RealType>& __x)
1009 {
1010 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1011 typedef typename __istream_type::ios_base __ios_base;
1012
1013 const typename __ios_base::fmtflags __flags = __is.flags();
1014 __is.flags(__ios_base::dec | __ios_base::skipws);
1015
1016 _RealType __lambda_val, __mu_val, __nu_val;
1017 __is >> __lambda_val >> __mu_val >> __nu_val;
1018 __is >> __x._M_gd1;
1019 __is >> __x._M_gd2;
1020 __x.param(typename k_distribution<_RealType>::
1021 param_type(__lambda_val, __mu_val, __nu_val));
1022
1023 __is.flags(__flags);
1024 return __is;
1025 }
1026
1027
1028 template<typename _RealType>
1029 template<typename _OutputIterator,
1030 typename _UniformRandomNumberGenerator>
1031 void
1032 arcsine_distribution<_RealType>::
1033 __generate_impl(_OutputIterator __f, _OutputIterator __t,
1034 _UniformRandomNumberGenerator& __urng,
1035 const param_type& __p)
1036 {
1037 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
1038 result_type>)
1039
1040 result_type __dif = __p.b() - __p.a();
1041 result_type __sum = __p.a() + __p.b();
1042 while (__f != __t)
1043 {
1044 result_type __x = std::sin(this->_M_ud(__urng));
1045 *__f++ = (__x * __dif + __sum) / result_type(2);
1046 }
1047 }
1048
1049 template<typename _RealType, typename _CharT, typename _Traits>
1052 const arcsine_distribution<_RealType>& __x)
1053 {
1054 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1055 typedef typename __ostream_type::ios_base __ios_base;
1056
1057 const typename __ios_base::fmtflags __flags = __os.flags();
1058 const _CharT __fill = __os.fill();
1059 const std::streamsize __precision = __os.precision();
1060 const _CharT __space = __os.widen(' ');
1061 __os.flags(__ios_base::scientific | __ios_base::left);
1062 __os.fill(__space);
1064
1065 __os << __x.a() << __space << __x.b();
1066 __os << __space << __x._M_ud;
1067
1068 __os.flags(__flags);
1069 __os.fill(__fill);
1070 __os.precision(__precision);
1071 return __os;
1072 }
1073
1074 template<typename _RealType, typename _CharT, typename _Traits>
1077 arcsine_distribution<_RealType>& __x)
1078 {
1079 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1080 typedef typename __istream_type::ios_base __ios_base;
1081
1082 const typename __ios_base::fmtflags __flags = __is.flags();
1083 __is.flags(__ios_base::dec | __ios_base::skipws);
1084
1085 _RealType __a, __b;
1086 __is >> __a >> __b;
1087 __is >> __x._M_ud;
1088 __x.param(typename arcsine_distribution<_RealType>::
1089 param_type(__a, __b));
1090
1091 __is.flags(__flags);
1092 return __is;
1093 }
1094
1095
1096 template<typename _RealType>
1097 template<typename _UniformRandomNumberGenerator>
1098 typename hoyt_distribution<_RealType>::result_type
1099 hoyt_distribution<_RealType>::
1100 operator()(_UniformRandomNumberGenerator& __urng)
1101 {
1102 result_type __x = this->_M_ad(__urng);
1103 result_type __y = this->_M_ed(__urng);
1104 return (result_type(2) * this->q()
1105 / (result_type(1) + this->q() * this->q()))
1106 * std::sqrt(this->omega() * __x * __y);
1107 }
1108
1109 template<typename _RealType>
1110 template<typename _UniformRandomNumberGenerator>
1111 typename hoyt_distribution<_RealType>::result_type
1112 hoyt_distribution<_RealType>::
1113 operator()(_UniformRandomNumberGenerator& __urng,
1114 const param_type& __p)
1115 {
1116 result_type __q2 = __p.q() * __p.q();
1117 result_type __num = result_type(0.5L) * (result_type(1) + __q2);
1118 typename __gnu_cxx::arcsine_distribution<result_type>::param_type
1119 __pa(__num, __num / __q2);
1120 result_type __x = this->_M_ad(__pa, __urng);
1121 result_type __y = this->_M_ed(__urng);
1122 return (result_type(2) * __p.q() / (result_type(1) + __q2))
1123 * std::sqrt(__p.omega() * __x * __y);
1124 }
1125
1126 template<typename _RealType>
1127 template<typename _OutputIterator,
1128 typename _UniformRandomNumberGenerator>
1129 void
1130 hoyt_distribution<_RealType>::
1131 __generate_impl(_OutputIterator __f, _OutputIterator __t,
1132 _UniformRandomNumberGenerator& __urng,
1133 const param_type& __p)
1134 {
1135 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
1136 result_type>)
1137
1138 result_type __2q = result_type(2) * __p.q();
1139 result_type __q2 = __p.q() * __p.q();
1140 result_type __q2p1 = result_type(1) + __q2;
1141 result_type __num = result_type(0.5L) * __q2p1;
1142 result_type __omega = __p.omega();
1143 typename __gnu_cxx::arcsine_distribution<result_type>::param_type
1144 __pa(__num, __num / __q2);
1145 while (__f != __t)
1146 {
1147 result_type __x = this->_M_ad(__pa, __urng);
1148 result_type __y = this->_M_ed(__urng);
1149 *__f++ = (__2q / __q2p1) * std::sqrt(__omega * __x * __y);
1150 }
1151 }
1152
1153 template<typename _RealType, typename _CharT, typename _Traits>
1156 const hoyt_distribution<_RealType>& __x)
1157 {
1158 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1159 typedef typename __ostream_type::ios_base __ios_base;
1160
1161 const typename __ios_base::fmtflags __flags = __os.flags();
1162 const _CharT __fill = __os.fill();
1163 const std::streamsize __precision = __os.precision();
1164 const _CharT __space = __os.widen(' ');
1165 __os.flags(__ios_base::scientific | __ios_base::left);
1166 __os.fill(__space);
1168
1169 __os << __x.q() << __space << __x.omega();
1170 __os << __space << __x._M_ad;
1171 __os << __space << __x._M_ed;
1172
1173 __os.flags(__flags);
1174 __os.fill(__fill);
1175 __os.precision(__precision);
1176 return __os;
1177 }
1178
1179 template<typename _RealType, typename _CharT, typename _Traits>
1182 hoyt_distribution<_RealType>& __x)
1183 {
1184 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1185 typedef typename __istream_type::ios_base __ios_base;
1186
1187 const typename __ios_base::fmtflags __flags = __is.flags();
1188 __is.flags(__ios_base::dec | __ios_base::skipws);
1189
1190 _RealType __q, __omega;
1191 __is >> __q >> __omega;
1192 __is >> __x._M_ad;
1193 __is >> __x._M_ed;
1194 __x.param(typename hoyt_distribution<_RealType>::
1195 param_type(__q, __omega));
1196
1197 __is.flags(__flags);
1198 return __is;
1199 }
1200
1201
1202 template<typename _RealType>
1203 template<typename _OutputIterator,
1204 typename _UniformRandomNumberGenerator>
1205 void
1206 triangular_distribution<_RealType>::
1207 __generate_impl(_OutputIterator __f, _OutputIterator __t,
1208 _UniformRandomNumberGenerator& __urng,
1209 const param_type& __param)
1210 {
1211 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
1212 result_type>)
1213
1214 while (__f != __t)
1215 *__f++ = this->operator()(__urng, __param);
1216 }
1217
1218 template<typename _RealType, typename _CharT, typename _Traits>
1219 std::basic_ostream<_CharT, _Traits>&
1220 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1221 const __gnu_cxx::triangular_distribution<_RealType>& __x)
1222 {
1223 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1224 typedef typename __ostream_type::ios_base __ios_base;
1225
1226 const typename __ios_base::fmtflags __flags = __os.flags();
1227 const _CharT __fill = __os.fill();
1228 const std::streamsize __precision = __os.precision();
1229 const _CharT __space = __os.widen(' ');
1230 __os.flags(__ios_base::scientific | __ios_base::left);
1231 __os.fill(__space);
1233
1234 __os << __x.a() << __space << __x.b() << __space << __x.c();
1235
1236 __os.flags(__flags);
1237 __os.fill(__fill);
1238 __os.precision(__precision);
1239 return __os;
1240 }
1241
1242 template<typename _RealType, typename _CharT, typename _Traits>
1245 __gnu_cxx::triangular_distribution<_RealType>& __x)
1246 {
1247 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1248 typedef typename __istream_type::ios_base __ios_base;
1249
1250 const typename __ios_base::fmtflags __flags = __is.flags();
1251 __is.flags(__ios_base::dec | __ios_base::skipws);
1252
1253 _RealType __a, __b, __c;
1254 __is >> __a >> __b >> __c;
1255 __x.param(typename __gnu_cxx::triangular_distribution<_RealType>::
1256 param_type(__a, __b, __c));
1257
1258 __is.flags(__flags);
1259 return __is;
1260 }
1261
1262
1263 template<typename _RealType>
1264 template<typename _UniformRandomNumberGenerator>
1265 typename von_mises_distribution<_RealType>::result_type
1266 von_mises_distribution<_RealType>::
1267 operator()(_UniformRandomNumberGenerator& __urng,
1268 const param_type& __p)
1269 {
1270 const result_type __pi
1271 = __gnu_cxx::__math_constants<result_type>::__pi;
1272 std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1273 __aurng(__urng);
1274
1275 result_type __f;
1276 while (1)
1277 {
1278 result_type __rnd = std::cos(__pi * __aurng());
1279 __f = (result_type(1) + __p._M_r * __rnd) / (__p._M_r + __rnd);
1280 result_type __c = __p._M_kappa * (__p._M_r - __f);
1281
1282 result_type __rnd2 = __aurng();
1283 if (__c * (result_type(2) - __c) > __rnd2)
1284 break;
1285 if (std::log(__c / __rnd2) >= __c - result_type(1))
1286 break;
1287 }
1288
1289 result_type __res = std::acos(__f);
1290#if _GLIBCXX_USE_C99_MATH_TR1
1291 __res = std::copysign(__res, __aurng() - result_type(0.5));
1292#else
1293 if (__aurng() < result_type(0.5))
1294 __res = -__res;
1295#endif
1296 __res += __p._M_mu;
1297 if (__res > __pi)
1298 __res -= result_type(2) * __pi;
1299 else if (__res < -__pi)
1300 __res += result_type(2) * __pi;
1301 return __res;
1302 }
1303
1304 template<typename _RealType>
1305 template<typename _OutputIterator,
1306 typename _UniformRandomNumberGenerator>
1307 void
1308 von_mises_distribution<_RealType>::
1309 __generate_impl(_OutputIterator __f, _OutputIterator __t,
1310 _UniformRandomNumberGenerator& __urng,
1311 const param_type& __param)
1312 {
1313 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
1314 result_type>)
1315
1316 while (__f != __t)
1317 *__f++ = this->operator()(__urng, __param);
1318 }
1319
1320 template<typename _RealType, typename _CharT, typename _Traits>
1321 std::basic_ostream<_CharT, _Traits>&
1322 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1323 const __gnu_cxx::von_mises_distribution<_RealType>& __x)
1324 {
1325 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1326 typedef typename __ostream_type::ios_base __ios_base;
1327
1328 const typename __ios_base::fmtflags __flags = __os.flags();
1329 const _CharT __fill = __os.fill();
1330 const std::streamsize __precision = __os.precision();
1331 const _CharT __space = __os.widen(' ');
1332 __os.flags(__ios_base::scientific | __ios_base::left);
1333 __os.fill(__space);
1335
1336 __os << __x.mu() << __space << __x.kappa();
1337
1338 __os.flags(__flags);
1339 __os.fill(__fill);
1340 __os.precision(__precision);
1341 return __os;
1342 }
1343
1344 template<typename _RealType, typename _CharT, typename _Traits>
1347 __gnu_cxx::von_mises_distribution<_RealType>& __x)
1348 {
1349 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1350 typedef typename __istream_type::ios_base __ios_base;
1351
1352 const typename __ios_base::fmtflags __flags = __is.flags();
1353 __is.flags(__ios_base::dec | __ios_base::skipws);
1354
1355 _RealType __mu, __kappa;
1356 __is >> __mu >> __kappa;
1357 __x.param(typename __gnu_cxx::von_mises_distribution<_RealType>::
1358 param_type(__mu, __kappa));
1359
1360 __is.flags(__flags);
1361 return __is;
1362 }
1363
1364
1365 template<typename _UIntType>
1366 template<typename _UniformRandomNumberGenerator>
1367 typename hypergeometric_distribution<_UIntType>::result_type
1368 hypergeometric_distribution<_UIntType>::
1369 operator()(_UniformRandomNumberGenerator& __urng,
1370 const param_type& __param)
1371 {
1372 std::__detail::_Adaptor<_UniformRandomNumberGenerator, double>
1373 __aurng(__urng);
1374
1375 result_type __a = __param.successful_size();
1376 result_type __b = __param.total_size();
1377 result_type __k = 0;
1378
1379 if (__param.total_draws() < __param.total_size() / 2)
1380 {
1381 for (result_type __i = 0; __i < __param.total_draws(); ++__i)
1382 {
1383 if (__b * __aurng() < __a)
1384 {
1385 ++__k;
1386 if (__k == __param.successful_size())
1387 return __k;
1388 --__a;
1389 }
1390 --__b;
1391 }
1392 return __k;
1393 }
1394 else
1395 {
1396 for (result_type __i = 0; __i < __param.unsuccessful_size(); ++__i)
1397 {
1398 if (__b * __aurng() < __a)
1399 {
1400 ++__k;
1401 if (__k == __param.successful_size())
1402 return __param.successful_size() - __k;
1403 --__a;
1404 }
1405 --__b;
1406 }
1407 return __param.successful_size() - __k;
1408 }
1409 }
1410
1411 template<typename _UIntType>
1412 template<typename _OutputIterator,
1413 typename _UniformRandomNumberGenerator>
1414 void
1415 hypergeometric_distribution<_UIntType>::
1416 __generate_impl(_OutputIterator __f, _OutputIterator __t,
1417 _UniformRandomNumberGenerator& __urng,
1418 const param_type& __param)
1419 {
1420 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
1421 result_type>)
1422
1423 while (__f != __t)
1424 *__f++ = this->operator()(__urng);
1425 }
1426
1427 template<typename _UIntType, typename _CharT, typename _Traits>
1428 std::basic_ostream<_CharT, _Traits>&
1429 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1430 const __gnu_cxx::hypergeometric_distribution<_UIntType>& __x)
1431 {
1432 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1433 typedef typename __ostream_type::ios_base __ios_base;
1434
1435 const typename __ios_base::fmtflags __flags = __os.flags();
1436 const _CharT __fill = __os.fill();
1437 const std::streamsize __precision = __os.precision();
1438 const _CharT __space = __os.widen(' ');
1439 __os.flags(__ios_base::scientific | __ios_base::left);
1440 __os.fill(__space);
1442
1443 __os << __x.total_size() << __space << __x.successful_size() << __space
1444 << __x.total_draws();
1445
1446 __os.flags(__flags);
1447 __os.fill(__fill);
1448 __os.precision(__precision);
1449 return __os;
1450 }
1451
1452 template<typename _UIntType, typename _CharT, typename _Traits>
1455 __gnu_cxx::hypergeometric_distribution<_UIntType>& __x)
1456 {
1457 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1458 typedef typename __istream_type::ios_base __ios_base;
1459
1460 const typename __ios_base::fmtflags __flags = __is.flags();
1461 __is.flags(__ios_base::dec | __ios_base::skipws);
1462
1463 _UIntType __total_size, __successful_size, __total_draws;
1464 __is >> __total_size >> __successful_size >> __total_draws;
1465 __x.param(typename __gnu_cxx::hypergeometric_distribution<_UIntType>::
1466 param_type(__total_size, __successful_size, __total_draws));
1467
1468 __is.flags(__flags);
1469 return __is;
1470 }
1471
1472
1473 template<typename _RealType>
1474 template<typename _UniformRandomNumberGenerator>
1475 typename logistic_distribution<_RealType>::result_type
1476 logistic_distribution<_RealType>::
1477 operator()(_UniformRandomNumberGenerator& __urng,
1478 const param_type& __p)
1479 {
1480 std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1481 __aurng(__urng);
1482
1483 result_type __arg = result_type(1);
1484 while (__arg == result_type(1) || __arg == result_type(0))
1485 __arg = __aurng();
1486 return __p.a()
1487 + __p.b() * std::log(__arg / (result_type(1) - __arg));
1488 }
1489
1490 template<typename _RealType>
1491 template<typename _OutputIterator,
1492 typename _UniformRandomNumberGenerator>
1493 void
1494 logistic_distribution<_RealType>::
1495 __generate_impl(_OutputIterator __f, _OutputIterator __t,
1496 _UniformRandomNumberGenerator& __urng,
1497 const param_type& __p)
1498 {
1499 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
1500 result_type>)
1501
1502 std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
1503 __aurng(__urng);
1504
1505 while (__f != __t)
1506 {
1507 result_type __arg = result_type(1);
1508 while (__arg == result_type(1) || __arg == result_type(0))
1509 __arg = __aurng();
1510 *__f++ = __p.a()
1511 + __p.b() * std::log(__arg / (result_type(1) - __arg));
1512 }
1513 }
1514
1515 template<typename _RealType, typename _CharT, typename _Traits>
1518 const logistic_distribution<_RealType>& __x)
1519 {
1520 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1521 typedef typename __ostream_type::ios_base __ios_base;
1522
1523 const typename __ios_base::fmtflags __flags = __os.flags();
1524 const _CharT __fill = __os.fill();
1525 const std::streamsize __precision = __os.precision();
1526 const _CharT __space = __os.widen(' ');
1527 __os.flags(__ios_base::scientific | __ios_base::left);
1528 __os.fill(__space);
1530
1531 __os << __x.a() << __space << __x.b();
1532
1533 __os.flags(__flags);
1534 __os.fill(__fill);
1535 __os.precision(__precision);
1536 return __os;
1537 }
1538
1539 template<typename _RealType, typename _CharT, typename _Traits>
1542 logistic_distribution<_RealType>& __x)
1543 {
1544 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1545 typedef typename __istream_type::ios_base __ios_base;
1546
1547 const typename __ios_base::fmtflags __flags = __is.flags();
1548 __is.flags(__ios_base::dec | __ios_base::skipws);
1549
1550 _RealType __a, __b;
1551 __is >> __a >> __b;
1552 __x.param(typename logistic_distribution<_RealType>::
1553 param_type(__a, __b));
1554
1555 __is.flags(__flags);
1556 return __is;
1557 }
1558
1559
1560 namespace {
1561
1562 // Helper class for the uniform_on_sphere_distribution generation
1563 // function.
1564 template<std::size_t _Dimen, typename _RealType>
1565 class uniform_on_sphere_helper
1566 {
1567 typedef typename uniform_on_sphere_distribution<_Dimen, _RealType>::
1568 result_type result_type;
1569
1570 public:
1571 template<typename _NormalDistribution,
1572 typename _UniformRandomNumberGenerator>
1573 result_type operator()(_NormalDistribution& __nd,
1574 _UniformRandomNumberGenerator& __urng)
1575 {
1576 result_type __ret;
1577 typename result_type::value_type __norm;
1578
1579 do
1580 {
1581 auto __sum = _RealType(0);
1582
1583 std::generate(__ret.begin(), __ret.end(),
1584 [&__nd, &__urng, &__sum](){
1585 _RealType __t = __nd(__urng);
1586 __sum += __t * __t;
1587 return __t; });
1588 __norm = std::sqrt(__sum);
1589 }
1590 while (__norm == _RealType(0) || ! __builtin_isfinite(__norm));
1591
1592 std::transform(__ret.begin(), __ret.end(), __ret.begin(),
1593 [__norm](_RealType __val){ return __val / __norm; });
1594
1595 return __ret;
1596 }
1597 };
1598
1599
1600 template<typename _RealType>
1601 class uniform_on_sphere_helper<2, _RealType>
1602 {
1603 typedef typename uniform_on_sphere_distribution<2, _RealType>::
1604 result_type result_type;
1605
1606 public:
1607 template<typename _NormalDistribution,
1608 typename _UniformRandomNumberGenerator>
1609 result_type operator()(_NormalDistribution&,
1610 _UniformRandomNumberGenerator& __urng)
1611 {
1612 result_type __ret;
1613 _RealType __sq;
1614 std::__detail::_Adaptor<_UniformRandomNumberGenerator,
1615 _RealType> __aurng(__urng);
1616
1617 do
1618 {
1619 __ret[0] = _RealType(2) * __aurng() - _RealType(1);
1620 __ret[1] = _RealType(2) * __aurng() - _RealType(1);
1621
1622 __sq = __ret[0] * __ret[0] + __ret[1] * __ret[1];
1623 }
1624 while (__sq == _RealType(0) || __sq > _RealType(1));
1625
1626#if _GLIBCXX_USE_C99_MATH_TR1
1627 // Yes, we do not just use sqrt(__sq) because hypot() is more
1628 // accurate.
1629 auto __norm = std::hypot(__ret[0], __ret[1]);
1630#else
1631 auto __norm = std::sqrt(__sq);
1632#endif
1633 __ret[0] /= __norm;
1634 __ret[1] /= __norm;
1635
1636 return __ret;
1637 }
1638 };
1639
1640 }
1641
1642
1643 template<std::size_t _Dimen, typename _RealType>
1644 template<typename _UniformRandomNumberGenerator>
1645 typename uniform_on_sphere_distribution<_Dimen, _RealType>::result_type
1646 uniform_on_sphere_distribution<_Dimen, _RealType>::
1647 operator()(_UniformRandomNumberGenerator& __urng,
1648 const param_type& __p)
1649 {
1650 uniform_on_sphere_helper<_Dimen, _RealType> __helper;
1651 return __helper(_M_nd, __urng);
1652 }
1653
1654 template<std::size_t _Dimen, typename _RealType>
1655 template<typename _OutputIterator,
1656 typename _UniformRandomNumberGenerator>
1657 void
1658 uniform_on_sphere_distribution<_Dimen, _RealType>::
1659 __generate_impl(_OutputIterator __f, _OutputIterator __t,
1660 _UniformRandomNumberGenerator& __urng,
1661 const param_type& __param)
1662 {
1663 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
1664 result_type>)
1665
1666 while (__f != __t)
1667 *__f++ = this->operator()(__urng, __param);
1668 }
1669
1670 template<std::size_t _Dimen, typename _RealType, typename _CharT,
1671 typename _Traits>
1672 std::basic_ostream<_CharT, _Traits>&
1673 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1674 const __gnu_cxx::uniform_on_sphere_distribution<_Dimen,
1675 _RealType>& __x)
1676 {
1677 return __os << __x._M_nd;
1678 }
1679
1680 template<std::size_t _Dimen, typename _RealType, typename _CharT,
1681 typename _Traits>
1684 __gnu_cxx::uniform_on_sphere_distribution<_Dimen,
1685 _RealType>& __x)
1686 {
1687 return __is >> __x._M_nd;
1688 }
1689
1690
1691 namespace {
1692
1693 // Helper class for the uniform_inside_sphere_distribution generation
1694 // function.
1695 template<std::size_t _Dimen, bool _SmallDimen, typename _RealType>
1696 class uniform_inside_sphere_helper;
1697
1698 template<std::size_t _Dimen, typename _RealType>
1699 class uniform_inside_sphere_helper<_Dimen, false, _RealType>
1700 {
1701 using result_type
1702 = typename uniform_inside_sphere_distribution<_Dimen, _RealType>::
1703 result_type;
1704
1705 public:
1706 template<typename _UniformOnSphereDistribution,
1707 typename _UniformRandomNumberGenerator>
1708 result_type
1709 operator()(_UniformOnSphereDistribution& __uosd,
1710 _UniformRandomNumberGenerator& __urng,
1711 _RealType __radius)
1712 {
1713 std::__detail::_Adaptor<_UniformRandomNumberGenerator,
1714 _RealType> __aurng(__urng);
1715
1716 _RealType __pow = 1 / _RealType(_Dimen);
1717 _RealType __urt = __radius * std::pow(__aurng(), __pow);
1718 result_type __ret = __uosd(__aurng);
1719
1720 std::transform(__ret.begin(), __ret.end(), __ret.begin(),
1721 [__urt](_RealType __val)
1722 { return __val * __urt; });
1723
1724 return __ret;
1725 }
1726 };
1727
1728 // Helper class for the uniform_inside_sphere_distribution generation
1729 // function specialized for small dimensions.
1730 template<std::size_t _Dimen, typename _RealType>
1731 class uniform_inside_sphere_helper<_Dimen, true, _RealType>
1732 {
1733 using result_type
1734 = typename uniform_inside_sphere_distribution<_Dimen, _RealType>::
1735 result_type;
1736
1737 public:
1738 template<typename _UniformOnSphereDistribution,
1739 typename _UniformRandomNumberGenerator>
1740 result_type
1741 operator()(_UniformOnSphereDistribution&,
1742 _UniformRandomNumberGenerator& __urng,
1743 _RealType __radius)
1744 {
1745 result_type __ret;
1746 _RealType __sq;
1747 _RealType __radsq = __radius * __radius;
1748 std::__detail::_Adaptor<_UniformRandomNumberGenerator,
1749 _RealType> __aurng(__urng);
1750
1751 do
1752 {
1753 __sq = _RealType(0);
1754 for (int i = 0; i < _Dimen; ++i)
1755 {
1756 __ret[i] = _RealType(2) * __aurng() - _RealType(1);
1757 __sq += __ret[i] * __ret[i];
1758 }
1759 }
1760 while (__sq > _RealType(1));
1761
1762 for (int i = 0; i < _Dimen; ++i)
1763 __ret[i] *= __radius;
1764
1765 return __ret;
1766 }
1767 };
1768 } // namespace
1769
1770 //
1771 // Experiments have shown that rejection is more efficient than transform
1772 // for dimensions less than 8.
1773 //
1774 template<std::size_t _Dimen, typename _RealType>
1775 template<typename _UniformRandomNumberGenerator>
1776 typename uniform_inside_sphere_distribution<_Dimen, _RealType>::result_type
1777 uniform_inside_sphere_distribution<_Dimen, _RealType>::
1778 operator()(_UniformRandomNumberGenerator& __urng,
1779 const param_type& __p)
1780 {
1781 uniform_inside_sphere_helper<_Dimen, _Dimen < 8, _RealType> __helper;
1782 return __helper(_M_uosd, __urng, __p.radius());
1783 }
1784
1785 template<std::size_t _Dimen, typename _RealType>
1786 template<typename _OutputIterator,
1787 typename _UniformRandomNumberGenerator>
1788 void
1789 uniform_inside_sphere_distribution<_Dimen, _RealType>::
1790 __generate_impl(_OutputIterator __f, _OutputIterator __t,
1791 _UniformRandomNumberGenerator& __urng,
1792 const param_type& __param)
1793 {
1794 __glibcxx_function_requires(_OutputIteratorConcept<_OutputIterator,
1795 result_type>)
1796
1797 while (__f != __t)
1798 *__f++ = this->operator()(__urng, __param);
1799 }
1800
1801 template<std::size_t _Dimen, typename _RealType, typename _CharT,
1802 typename _Traits>
1803 std::basic_ostream<_CharT, _Traits>&
1804 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1805 const __gnu_cxx::uniform_inside_sphere_distribution<_Dimen,
1806 _RealType>& __x)
1807 {
1808 typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1809 typedef typename __ostream_type::ios_base __ios_base;
1810
1811 const typename __ios_base::fmtflags __flags = __os.flags();
1812 const _CharT __fill = __os.fill();
1813 const std::streamsize __precision = __os.precision();
1814 const _CharT __space = __os.widen(' ');
1815 __os.flags(__ios_base::scientific | __ios_base::left);
1816 __os.fill(__space);
1818
1819 __os << __x.radius() << __space << __x._M_uosd;
1820
1821 __os.flags(__flags);
1822 __os.fill(__fill);
1823 __os.precision(__precision);
1824
1825 return __os;
1826 }
1827
1828 template<std::size_t _Dimen, typename _RealType, typename _CharT,
1829 typename _Traits>
1832 __gnu_cxx::uniform_inside_sphere_distribution<_Dimen,
1833 _RealType>& __x)
1834 {
1835 typedef std::basic_istream<_CharT, _Traits> __istream_type;
1836 typedef typename __istream_type::ios_base __ios_base;
1837
1838 const typename __ios_base::fmtflags __flags = __is.flags();
1839 __is.flags(__ios_base::dec | __ios_base::skipws);
1840
1841 _RealType __radius_val;
1842 __is >> __radius_val >> __x._M_uosd;
1843 __x.param(typename uniform_inside_sphere_distribution<_Dimen, _RealType>::
1844 param_type(__radius_val));
1845
1846 __is.flags(__flags);
1847
1848 return __is;
1849 }
1850
1851_GLIBCXX_END_NAMESPACE_VERSION
1852} // namespace __gnu_cxx
1853
1854
1855#endif // _EXT_RANDOM_TCC
complex< _Tp > sin(const complex< _Tp > &)
Return complex sine of z.
Definition: complex:1120
complex< _Tp > log(const complex< _Tp > &)
Return complex natural logarithm of z.
Definition: complex:1085
complex< _Tp > exp(const complex< _Tp > &)
Return complex base e exponential of z.
Definition: complex:1058
complex< _Tp > pow(const complex< _Tp > &, int)
Return x to the y'th power.
Definition: complex:1280
complex< _Tp > cos(const complex< _Tp > &)
Return complex cosine of z.
Definition: complex:1002
complex< _Tp > sqrt(const complex< _Tp > &)
Return complex square root of z.
Definition: complex:1194
ISO C++ entities toplevel namespace is std.
ptrdiff_t streamsize
Integral type for I/O operation counts and buffer sizes.
Definition: postypes.h:68
constexpr void advance(_InputIterator &__i, _Distance __n)
A generalization of pointer arithmetic.
std::basic_istream< _CharT, _Traits > & operator>>(std::basic_istream< _CharT, _Traits > &__is, bitset< _Nb > &__x)
Global I/O operators for bitsets.
Definition: bitset:1593
std::basic_ostream< _CharT, _Traits > & operator<<(std::basic_ostream< _CharT, _Traits > &__os, const bitset< _Nb > &__x)
Global I/O operators for bitsets.
Definition: bitset:1683
std::complex< _Tp > acos(const std::complex< _Tp > &)
acos(__z) [8.1.2].
Definition: complex:2203
GNU extensions for public use.
A standard container for storing a fixed size sequence of elements.
Definition: array:95
char_type widen(char __c) const
Widens characters.
Definition: basic_ios.h:449
char_type fill() const
Retrieves the empty character.
Definition: basic_ios.h:370
Template class basic_istream.
Definition: istream:61
Template class basic_ostream.
Definition: ostream:61
Properties of fundamental types.
Definition: limits:313
streamsize precision() const
Flags access.
Definition: ios_base.h:732
fmtflags flags() const
Access to format flags.
Definition: ios_base.h:662