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