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