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