libstdc++
ext/random.tcc
Go to the documentation of this file.
1 // Random number extensions -*- C++ -*-
2 
3 // Copyright (C) 2012-2017 Free Software Foundation, Inc.
4 //
5 // This file is part of the GNU ISO C++ Library. This library is free
6 // software; you can redistribute it and/or modify it under the
7 // terms of the GNU General Public License as published by the
8 // Free Software Foundation; either version 3, or (at your option)
9 // any later version.
10 
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU General Public License for more details.
15 
16 // Under Section 7 of GPL version 3, you are granted additional
17 // permissions described in the GCC Runtime Library Exception, version
18 // 3.1, as published by the Free Software Foundation.
19 
20 // You should have received a copy of the GNU General Public License and
21 // a copy of the GCC Runtime Library Exception along with this program;
22 // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
23 // <http://www.gnu.org/licenses/>.
24 
25 /** @file ext/random.tcc
26  * This is an internal header file, included by other library headers.
27  * Do not attempt to use it directly. @headername{ext/random}
28  */
29 
30 #ifndef _EXT_RANDOM_TCC
31 #define _EXT_RANDOM_TCC 1
32 
33 #pragma GCC system_header
34 
35 namespace __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>
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(' ');
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>
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();
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>
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(' ');
476  __os.fill(__space);
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>
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();
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>
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(' ');
673  __os.fill(__space);
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>
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();
701 
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  {
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>
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(' ');
759  __os.fill(__space);
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>
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();
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(' ');
826  __os.fill(__space);
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>
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();
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(' ');
891  __os.fill(__space);
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>
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();
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  {
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>
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(' ');
987  __os.fill(__space);
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>
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();
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>
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(' ');
1057  __os.fill(__space);
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>
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();
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>
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(' ');
1161  __os.fill(__space);
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>
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();
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(' ');
1226  __os.fill(__space);
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>
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();
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(' ');
1328  __os.fill(__space);
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>
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();
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(' ');
1435  __os.fill(__space);
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>
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();
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>
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(' ');
1523  __os.fill(__space);
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>
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();
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>
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(' ');
1811  __os.fill(__space);
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>
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();
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
complex< _Tp > sqrt(const complex< _Tp > &)
Return complex square root of z.
Definition: complex:901
_OI fill_n(_OI __first, _Size __n, const _Tp &__value)
Fills the range [first,first+n) with copies of value.
Definition: stl_algobase.h:784
constexpr const _Tp * begin(initializer_list< _Tp > __ils) noexcept
Return an iterator pointing to the first element of the initializer_list.
Template class basic_ostream.
Definition: iosfwd:86
fmtflags flags() const
Access to format flags.
Definition: ios_base.h:619
complex< _Tp > cos(const complex< _Tp > &)
Return complex cosine of z.
Definition: complex:709
complex< _Tp > exp(const complex< _Tp > &)
Return complex base e exponential of z.
Definition: complex:765
streamsize precision() const
Flags access.
Definition: ios_base.h:689
ios_base & skipws(ios_base &__base)
Calls base.setf(ios_base::skipws).
Definition: ios_base.h:942
std::complex< _Tp > acos(const std::complex< _Tp > &)
acos(__z) [8.1.2].
Definition: complex:1617
constexpr const _Tp * end(initializer_list< _Tp > __ils) noexcept
Return an iterator pointing to one past the last element of the initializer_list. ...
ios_base & left(ios_base &__base)
Calls base.setf(ios_base::left, ios_base::adjustfield).
Definition: ios_base.h:999
ios_base & dec(ios_base &__base)
Calls base.setf(ios_base::dec, ios_base::basefield).
Definition: ios_base.h:1016
ISO C++ entities toplevel namespace is std.
Template class basic_istream.
Definition: iosfwd:83
char_type widen(char __c) const
Widens characters.
Definition: basic_ios.h:449
void generate(_ForwardIterator __first, _ForwardIterator __last, _Generator __gen)
Assign the result of a function object to each value in a sequence.
Definition: stl_algo.h:4414
ios_base & fixed(ios_base &__base)
Calls base.setf(ios_base::fixed, ios_base::floatfield).
Definition: ios_base.h:1041
std::basic_istream< _CharT, _Traits > & operator>>(std::basic_istream< _CharT, _Traits > &__is, bitset< _Nb > &__x)
Global I/O operators for bitsets.
Definition: bitset:1462
_OutputIterator transform(_InputIterator1 __first1, _InputIterator1 __last1, _InputIterator2 __first2, _OutputIterator __result, _BinaryOperation __binary_op)
Perform an operation on corresponding elements of two sequences.
Definition: stl_algo.h:4318
ios_base & scientific(ios_base &__base)
Calls base.setf(ios_base::scientific, ios_base::floatfield).
Definition: ios_base.h:1049
_GLIBCXX17_CONSTEXPR void advance(_InputIterator &__i, _Distance __n)
A generalization of pointer arithmetic.
complex< _Tp > sin(const complex< _Tp > &)
Return complex sine of z.
Definition: complex:827
complex< _Tp > pow(const complex< _Tp > &, int)
Return x to the y&#39;th power.
Definition: complex:987
ptrdiff_t streamsize
Integral type for I/O operation counts and buffer sizes.
Definition: postypes.h:98
complex< _Tp > log(const complex< _Tp > &)
Return complex natural logarithm of z.
Definition: complex:792
A standard container for storing a fixed size sequence of elements.
Definition: array:94
GNU extensions for public use.
char_type fill() const
Retrieves the empty character.
Definition: basic_ios.h:370
Properties of fundamental types.
Definition: limits:312