1// Random number extensions -*- C++ -*-
2
3// Copyright (C) 2012-2018 Free Software Foundation, Inc.
4//
5// This file is part of the GNU ISO C++ Library.  This library is free
6// software; you can redistribute it and/or modify it under the
7// terms of the GNU General Public License as published by the
8// Free Software Foundation; either version 3, or (at your option)
9// any later version.
10
11// This library is distributed in the hope that it will be useful,
12// but WITHOUT ANY WARRANTY; without even the implied warranty of
13// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14// GNU General Public License for more details.
15
16// Under Section 7 of GPL version 3, you are granted additional
17// permissions described in the GCC Runtime Library Exception, version
18// 3.1, as published by the Free Software Foundation.
19
20// You should have received a copy of the GNU General Public License and
21// a copy of the GCC Runtime Library Exception along with this program;
22// see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
23// <http://www.gnu.org/licenses/>.
24
25/** @file ext/random
26 *  This file is a GNU extension to the Standard C++ Library.
27 */
28
29#ifndef _EXT_RANDOM
30#define _EXT_RANDOM 1
31
32#pragma GCC system_header
33
34#if __cplusplus < 201103L
35# include <bits/c++0x_warning.h>
36#else
37
38#include <random>
39#include <algorithm>
40#include <array>
41#include <ext/cmath>
42#ifdef __SSE2__
43# include <emmintrin.h>
44#endif
45
46#if defined(_GLIBCXX_USE_C99_STDINT_TR1) && defined(UINT32_C)
47
48namespace __gnu_cxx _GLIBCXX_VISIBILITY(default)
49{
50_GLIBCXX_BEGIN_NAMESPACE_VERSION
51
52#if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
53
54  /* Mersenne twister implementation optimized for vector operations.
55   *
56   * Reference: http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/
57   */
58  template<typename _UIntType, size_t __m,
59	   size_t __pos1, size_t __sl1, size_t __sl2,
60	   size_t __sr1, size_t __sr2,
61	   uint32_t __msk1, uint32_t __msk2,
62	   uint32_t __msk3, uint32_t __msk4,
63	   uint32_t __parity1, uint32_t __parity2,
64	   uint32_t __parity3, uint32_t __parity4>
65    class simd_fast_mersenne_twister_engine
66    {
67      static_assert(std::is_unsigned<_UIntType>::value, "template argument "
68		    "substituting _UIntType not an unsigned integral type");
69      static_assert(__sr1 < 32, "first right shift too large");
70      static_assert(__sr2 < 16, "second right shift too large");
71      static_assert(__sl1 < 32, "first left shift too large");
72      static_assert(__sl2 < 16, "second left shift too large");
73
74    public:
75      typedef _UIntType result_type;
76
77    private:
78      static constexpr size_t m_w = sizeof(result_type) * 8;
79      static constexpr size_t _M_nstate = __m / 128 + 1;
80      static constexpr size_t _M_nstate32 = _M_nstate * 4;
81
82      static_assert(std::is_unsigned<_UIntType>::value, "template argument "
83		    "substituting _UIntType not an unsigned integral type");
84      static_assert(__pos1 < _M_nstate, "POS1 not smaller than state size");
85      static_assert(16 % sizeof(_UIntType) == 0,
86		    "UIntType size must divide 16");
87
88    public:
89      static constexpr size_t state_size = _M_nstate * (16
90							/ sizeof(result_type));
91      static constexpr result_type default_seed = 5489u;
92
93      // constructors and member function
94      explicit
95      simd_fast_mersenne_twister_engine(result_type __sd = default_seed)
96      { seed(__sd); }
97
98      template<typename _Sseq, typename = typename
99	std::enable_if<!std::is_same<_Sseq,
100				     simd_fast_mersenne_twister_engine>::value>
101	       ::type>
102	explicit
103	simd_fast_mersenne_twister_engine(_Sseq& __q)
104	{ seed(__q); }
105
106      void
107      seed(result_type __sd = default_seed);
108
109      template<typename _Sseq>
110	typename std::enable_if<std::is_class<_Sseq>::value>::type
111	seed(_Sseq& __q);
112
113      static constexpr result_type
114      min()
115      { return 0; }
116
117      static constexpr result_type
118      max()
119      { return std::numeric_limits<result_type>::max(); }
120
121      void
122      discard(unsigned long long __z);
123
124      result_type
125      operator()()
126      {
127	if (__builtin_expect(_M_pos >= state_size, 0))
128	  _M_gen_rand();
129
130	return _M_stateT[_M_pos++];
131      }
132
133      template<typename _UIntType_2, size_t __m_2,
134	       size_t __pos1_2, size_t __sl1_2, size_t __sl2_2,
135	       size_t __sr1_2, size_t __sr2_2,
136	       uint32_t __msk1_2, uint32_t __msk2_2,
137	       uint32_t __msk3_2, uint32_t __msk4_2,
138	       uint32_t __parity1_2, uint32_t __parity2_2,
139	       uint32_t __parity3_2, uint32_t __parity4_2>
140	friend bool
141	operator==(const simd_fast_mersenne_twister_engine<_UIntType_2,
142		   __m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2,
143		   __msk1_2, __msk2_2, __msk3_2, __msk4_2,
144		   __parity1_2, __parity2_2, __parity3_2, __parity4_2>& __lhs,
145		   const simd_fast_mersenne_twister_engine<_UIntType_2,
146		   __m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2,
147		   __msk1_2, __msk2_2, __msk3_2, __msk4_2,
148		   __parity1_2, __parity2_2, __parity3_2, __parity4_2>& __rhs);
149
150      template<typename _UIntType_2, size_t __m_2,
151	       size_t __pos1_2, size_t __sl1_2, size_t __sl2_2,
152	       size_t __sr1_2, size_t __sr2_2,
153	       uint32_t __msk1_2, uint32_t __msk2_2,
154	       uint32_t __msk3_2, uint32_t __msk4_2,
155	       uint32_t __parity1_2, uint32_t __parity2_2,
156	       uint32_t __parity3_2, uint32_t __parity4_2,
157	       typename _CharT, typename _Traits>
158	friend std::basic_ostream<_CharT, _Traits>&
159	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
160		   const __gnu_cxx::simd_fast_mersenne_twister_engine
161		   <_UIntType_2,
162		   __m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2,
163		   __msk1_2, __msk2_2, __msk3_2, __msk4_2,
164		   __parity1_2, __parity2_2, __parity3_2, __parity4_2>& __x);
165
166      template<typename _UIntType_2, size_t __m_2,
167	       size_t __pos1_2, size_t __sl1_2, size_t __sl2_2,
168	       size_t __sr1_2, size_t __sr2_2,
169	       uint32_t __msk1_2, uint32_t __msk2_2,
170	       uint32_t __msk3_2, uint32_t __msk4_2,
171	       uint32_t __parity1_2, uint32_t __parity2_2,
172	       uint32_t __parity3_2, uint32_t __parity4_2,
173	       typename _CharT, typename _Traits>
174	friend std::basic_istream<_CharT, _Traits>&
175	operator>>(std::basic_istream<_CharT, _Traits>& __is,
176		   __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType_2,
177		   __m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2,
178		   __msk1_2, __msk2_2, __msk3_2, __msk4_2,
179		   __parity1_2, __parity2_2, __parity3_2, __parity4_2>& __x);
180
181    private:
182      union
183      {
184#ifdef __SSE2__
185	__m128i _M_state[_M_nstate];
186#endif
187#ifdef __ARM_NEON
188#ifdef __aarch64__
189	__Uint32x4_t _M_state[_M_nstate];
190#endif
191#endif
192	uint32_t _M_state32[_M_nstate32];
193	result_type _M_stateT[state_size];
194      } __attribute__ ((__aligned__ (16)));
195      size_t _M_pos;
196
197      void _M_gen_rand(void);
198      void _M_period_certification();
199  };
200
201
202  template<typename _UIntType, size_t __m,
203	   size_t __pos1, size_t __sl1, size_t __sl2,
204	   size_t __sr1, size_t __sr2,
205	   uint32_t __msk1, uint32_t __msk2,
206	   uint32_t __msk3, uint32_t __msk4,
207	   uint32_t __parity1, uint32_t __parity2,
208	   uint32_t __parity3, uint32_t __parity4>
209    inline bool
210    operator!=(const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
211	       __m, __pos1, __sl1, __sl2, __sr1, __sr2, __msk1, __msk2, __msk3,
212	       __msk4, __parity1, __parity2, __parity3, __parity4>& __lhs,
213	       const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
214	       __m, __pos1, __sl1, __sl2, __sr1, __sr2, __msk1, __msk2, __msk3,
215	       __msk4, __parity1, __parity2, __parity3, __parity4>& __rhs)
216    { return !(__lhs == __rhs); }
217
218
219  /* Definitions for the SIMD-oriented Fast Mersenne Twister as defined
220   * in the C implementation by Daito and Matsumoto, as both a 32-bit
221   * and 64-bit version.
222   */
223  typedef simd_fast_mersenne_twister_engine<uint32_t, 607, 2,
224					    15, 3, 13, 3,
225					    0xfdff37ffU, 0xef7f3f7dU,
226					    0xff777b7dU, 0x7ff7fb2fU,
227					    0x00000001U, 0x00000000U,
228					    0x00000000U, 0x5986f054U>
229    sfmt607;
230
231  typedef simd_fast_mersenne_twister_engine<uint64_t, 607, 2,
232					    15, 3, 13, 3,
233					    0xfdff37ffU, 0xef7f3f7dU,
234					    0xff777b7dU, 0x7ff7fb2fU,
235					    0x00000001U, 0x00000000U,
236					    0x00000000U, 0x5986f054U>
237    sfmt607_64;
238
239
240  typedef simd_fast_mersenne_twister_engine<uint32_t, 1279, 7,
241					    14, 3, 5, 1,
242					    0xf7fefffdU, 0x7fefcfffU,
243					    0xaff3ef3fU, 0xb5ffff7fU,
244					    0x00000001U, 0x00000000U,
245					    0x00000000U, 0x20000000U>
246    sfmt1279;
247
248  typedef simd_fast_mersenne_twister_engine<uint64_t, 1279, 7,
249					    14, 3, 5, 1,
250					    0xf7fefffdU, 0x7fefcfffU,
251					    0xaff3ef3fU, 0xb5ffff7fU,
252					    0x00000001U, 0x00000000U,
253					    0x00000000U, 0x20000000U>
254    sfmt1279_64;
255
256
257  typedef simd_fast_mersenne_twister_engine<uint32_t, 2281, 12,
258					    19, 1, 5, 1,
259					    0xbff7ffbfU, 0xfdfffffeU,
260					    0xf7ffef7fU, 0xf2f7cbbfU,
261					    0x00000001U, 0x00000000U,
262					    0x00000000U, 0x41dfa600U>
263    sfmt2281;
264
265  typedef simd_fast_mersenne_twister_engine<uint64_t, 2281, 12,
266					    19, 1, 5, 1,
267					    0xbff7ffbfU, 0xfdfffffeU,
268					    0xf7ffef7fU, 0xf2f7cbbfU,
269					    0x00000001U, 0x00000000U,
270					    0x00000000U, 0x41dfa600U>
271    sfmt2281_64;
272
273
274  typedef simd_fast_mersenne_twister_engine<uint32_t, 4253, 17,
275					    20, 1, 7, 1,
276					    0x9f7bffffU, 0x9fffff5fU,
277					    0x3efffffbU, 0xfffff7bbU,
278					    0xa8000001U, 0xaf5390a3U,
279					    0xb740b3f8U, 0x6c11486dU>
280    sfmt4253;
281
282  typedef simd_fast_mersenne_twister_engine<uint64_t, 4253, 17,
283					    20, 1, 7, 1,
284					    0x9f7bffffU, 0x9fffff5fU,
285					    0x3efffffbU, 0xfffff7bbU,
286					    0xa8000001U, 0xaf5390a3U,
287					    0xb740b3f8U, 0x6c11486dU>
288    sfmt4253_64;
289
290
291  typedef simd_fast_mersenne_twister_engine<uint32_t, 11213, 68,
292					    14, 3, 7, 3,
293					    0xeffff7fbU, 0xffffffefU,
294					    0xdfdfbfffU, 0x7fffdbfdU,
295					    0x00000001U, 0x00000000U,
296					    0xe8148000U, 0xd0c7afa3U>
297    sfmt11213;
298
299  typedef simd_fast_mersenne_twister_engine<uint64_t, 11213, 68,
300					    14, 3, 7, 3,
301					    0xeffff7fbU, 0xffffffefU,
302					    0xdfdfbfffU, 0x7fffdbfdU,
303					    0x00000001U, 0x00000000U,
304					    0xe8148000U, 0xd0c7afa3U>
305    sfmt11213_64;
306
307
308  typedef simd_fast_mersenne_twister_engine<uint32_t, 19937, 122,
309					    18, 1, 11, 1,
310					    0xdfffffefU, 0xddfecb7fU,
311					    0xbffaffffU, 0xbffffff6U,
312					    0x00000001U, 0x00000000U,
313					    0x00000000U, 0x13c9e684U>
314    sfmt19937;
315
316  typedef simd_fast_mersenne_twister_engine<uint64_t, 19937, 122,
317					    18, 1, 11, 1,
318					    0xdfffffefU, 0xddfecb7fU,
319					    0xbffaffffU, 0xbffffff6U,
320					    0x00000001U, 0x00000000U,
321					    0x00000000U, 0x13c9e684U>
322    sfmt19937_64;
323
324
325  typedef simd_fast_mersenne_twister_engine<uint32_t, 44497, 330,
326					    5, 3, 9, 3,
327					    0xeffffffbU, 0xdfbebfffU,
328					    0xbfbf7befU, 0x9ffd7bffU,
329					    0x00000001U, 0x00000000U,
330					    0xa3ac4000U, 0xecc1327aU>
331    sfmt44497;
332
333  typedef simd_fast_mersenne_twister_engine<uint64_t, 44497, 330,
334					    5, 3, 9, 3,
335					    0xeffffffbU, 0xdfbebfffU,
336					    0xbfbf7befU, 0x9ffd7bffU,
337					    0x00000001U, 0x00000000U,
338					    0xa3ac4000U, 0xecc1327aU>
339    sfmt44497_64;
340
341
342  typedef simd_fast_mersenne_twister_engine<uint32_t, 86243, 366,
343					    6, 7, 19, 1,
344					    0xfdbffbffU, 0xbff7ff3fU,
345					    0xfd77efffU, 0xbf9ff3ffU,
346					    0x00000001U, 0x00000000U,
347					    0x00000000U, 0xe9528d85U>
348    sfmt86243;
349
350  typedef simd_fast_mersenne_twister_engine<uint64_t, 86243, 366,
351					    6, 7, 19, 1,
352					    0xfdbffbffU, 0xbff7ff3fU,
353					    0xfd77efffU, 0xbf9ff3ffU,
354					    0x00000001U, 0x00000000U,
355					    0x00000000U, 0xe9528d85U>
356    sfmt86243_64;
357
358
359  typedef simd_fast_mersenne_twister_engine<uint32_t, 132049, 110,
360					    19, 1, 21, 1,
361					    0xffffbb5fU, 0xfb6ebf95U,
362					    0xfffefffaU, 0xcff77fffU,
363					    0x00000001U, 0x00000000U,
364					    0xcb520000U, 0xc7e91c7dU>
365    sfmt132049;
366
367  typedef simd_fast_mersenne_twister_engine<uint64_t, 132049, 110,
368					    19, 1, 21, 1,
369					    0xffffbb5fU, 0xfb6ebf95U,
370					    0xfffefffaU, 0xcff77fffU,
371					    0x00000001U, 0x00000000U,
372					    0xcb520000U, 0xc7e91c7dU>
373    sfmt132049_64;
374
375
376  typedef simd_fast_mersenne_twister_engine<uint32_t, 216091, 627,
377					    11, 3, 10, 1,
378					    0xbff7bff7U, 0xbfffffffU,
379					    0xbffffa7fU, 0xffddfbfbU,
380					    0xf8000001U, 0x89e80709U,
381					    0x3bd2b64bU, 0x0c64b1e4U>
382    sfmt216091;
383
384  typedef simd_fast_mersenne_twister_engine<uint64_t, 216091, 627,
385					    11, 3, 10, 1,
386					    0xbff7bff7U, 0xbfffffffU,
387					    0xbffffa7fU, 0xffddfbfbU,
388					    0xf8000001U, 0x89e80709U,
389					    0x3bd2b64bU, 0x0c64b1e4U>
390    sfmt216091_64;
391
392#endif // __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
393
394  /**
395   * @brief A beta continuous distribution for random numbers.
396   *
397   * The formula for the beta probability density function is:
398   * @f[
399   *     p(x|\alpha,\beta) = \frac{1}{B(\alpha,\beta)}
400   *                         x^{\alpha - 1} (1 - x)^{\beta - 1}
401   * @f]
402   */
403  template<typename _RealType = double>
404    class beta_distribution
405    {
406      static_assert(std::is_floating_point<_RealType>::value,
407		    "template argument not a floating point type");
408
409    public:
410      /** The type of the range of the distribution. */
411      typedef _RealType result_type;
412
413      /** Parameter type. */
414      struct param_type
415      {
416	typedef beta_distribution<_RealType> distribution_type;
417	friend class beta_distribution<_RealType>;
418
419	explicit
420	param_type(_RealType __alpha_val = _RealType(1),
421		   _RealType __beta_val = _RealType(1))
422	: _M_alpha(__alpha_val), _M_beta(__beta_val)
423	{
424	  __glibcxx_assert(_M_alpha > _RealType(0));
425	  __glibcxx_assert(_M_beta > _RealType(0));
426	}
427
428	_RealType
429	alpha() const
430	{ return _M_alpha; }
431
432	_RealType
433	beta() const
434	{ return _M_beta; }
435
436	friend bool
437	operator==(const param_type& __p1, const param_type& __p2)
438	{ return (__p1._M_alpha == __p2._M_alpha
439		  && __p1._M_beta == __p2._M_beta); }
440
441	friend bool
442	operator!=(const param_type& __p1, const param_type& __p2)
443	{ return !(__p1 == __p2); }
444
445      private:
446	void
447	_M_initialize();
448
449	_RealType _M_alpha;
450	_RealType _M_beta;
451      };
452
453    public:
454      /**
455       * @brief Constructs a beta distribution with parameters
456       * @f$\alpha@f$ and @f$\beta@f$.
457       */
458      explicit
459      beta_distribution(_RealType __alpha_val = _RealType(1),
460			_RealType __beta_val = _RealType(1))
461      : _M_param(__alpha_val, __beta_val)
462      { }
463
464      explicit
465      beta_distribution(const param_type& __p)
466      : _M_param(__p)
467      { }
468
469      /**
470       * @brief Resets the distribution state.
471       */
472      void
473      reset()
474      { }
475
476      /**
477       * @brief Returns the @f$\alpha@f$ of the distribution.
478       */
479      _RealType
480      alpha() const
481      { return _M_param.alpha(); }
482
483      /**
484       * @brief Returns the @f$\beta@f$ of the distribution.
485       */
486      _RealType
487      beta() const
488      { return _M_param.beta(); }
489
490      /**
491       * @brief Returns the parameter set of the distribution.
492       */
493      param_type
494      param() const
495      { return _M_param; }
496
497      /**
498       * @brief Sets the parameter set of the distribution.
499       * @param __param The new parameter set of the distribution.
500       */
501      void
502      param(const param_type& __param)
503      { _M_param = __param; }
504
505      /**
506       * @brief Returns the greatest lower bound value of the distribution.
507       */
508      result_type
509      min() const
510      { return result_type(0); }
511
512      /**
513       * @brief Returns the least upper bound value of the distribution.
514       */
515      result_type
516      max() const
517      { return result_type(1); }
518
519      /**
520       * @brief Generating functions.
521       */
522      template<typename _UniformRandomNumberGenerator>
523	result_type
524	operator()(_UniformRandomNumberGenerator& __urng)
525	{ return this->operator()(__urng, _M_param); }
526
527      template<typename _UniformRandomNumberGenerator>
528	result_type
529	operator()(_UniformRandomNumberGenerator& __urng,
530		   const param_type& __p);
531
532      template<typename _ForwardIterator,
533	       typename _UniformRandomNumberGenerator>
534	void
535	__generate(_ForwardIterator __f, _ForwardIterator __t,
536		   _UniformRandomNumberGenerator& __urng)
537	{ this->__generate(__f, __t, __urng, _M_param); }
538
539      template<typename _ForwardIterator,
540	       typename _UniformRandomNumberGenerator>
541	void
542	__generate(_ForwardIterator __f, _ForwardIterator __t,
543		   _UniformRandomNumberGenerator& __urng,
544		   const param_type& __p)
545	{ this->__generate_impl(__f, __t, __urng, __p); }
546
547      template<typename _UniformRandomNumberGenerator>
548	void
549	__generate(result_type* __f, result_type* __t,
550		   _UniformRandomNumberGenerator& __urng,
551		   const param_type& __p)
552	{ this->__generate_impl(__f, __t, __urng, __p); }
553
554      /**
555       * @brief Return true if two beta distributions have the same
556       *        parameters and the sequences that would be generated
557       *        are equal.
558       */
559      friend bool
560      operator==(const beta_distribution& __d1,
561		 const beta_distribution& __d2)
562      { return __d1._M_param == __d2._M_param; }
563
564      /**
565       * @brief Inserts a %beta_distribution random number distribution
566       * @p __x into the output stream @p __os.
567       *
568       * @param __os An output stream.
569       * @param __x  A %beta_distribution random number distribution.
570       *
571       * @returns The output stream with the state of @p __x inserted or in
572       * an error state.
573       */
574      template<typename _RealType1, typename _CharT, typename _Traits>
575	friend std::basic_ostream<_CharT, _Traits>&
576	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
577		   const __gnu_cxx::beta_distribution<_RealType1>& __x);
578
579      /**
580       * @brief Extracts a %beta_distribution random number distribution
581       * @p __x from the input stream @p __is.
582       *
583       * @param __is An input stream.
584       * @param __x  A %beta_distribution random number generator engine.
585       *
586       * @returns The input stream with @p __x extracted or in an error state.
587       */
588      template<typename _RealType1, typename _CharT, typename _Traits>
589	friend std::basic_istream<_CharT, _Traits>&
590	operator>>(std::basic_istream<_CharT, _Traits>& __is,
591		   __gnu_cxx::beta_distribution<_RealType1>& __x);
592
593    private:
594      template<typename _ForwardIterator,
595	       typename _UniformRandomNumberGenerator>
596	void
597	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
598			_UniformRandomNumberGenerator& __urng,
599			const param_type& __p);
600
601      param_type _M_param;
602    };
603
604  /**
605   * @brief Return true if two beta distributions are different.
606   */
607  template<typename _RealType>
608    inline bool
609    operator!=(const __gnu_cxx::beta_distribution<_RealType>& __d1,
610	       const __gnu_cxx::beta_distribution<_RealType>& __d2)
611    { return !(__d1 == __d2); }
612
613
614  /**
615   * @brief A multi-variate normal continuous distribution for random numbers.
616   *
617   * The formula for the normal probability density function is
618   * @f[
619   *     p(\overrightarrow{x}|\overrightarrow{\mu },\Sigma) =
620   *       \frac{1}{\sqrt{(2\pi )^k\det(\Sigma))}}
621   *       e^{-\frac{1}{2}(\overrightarrow{x}-\overrightarrow{\mu})^\text{T}
622   *          \Sigma ^{-1}(\overrightarrow{x}-\overrightarrow{\mu})}
623   * @f]
624   *
625   * where @f$\overrightarrow{x}@f$ and @f$\overrightarrow{\mu}@f$ are
626   * vectors of dimension @f$k@f$ and @f$\Sigma@f$ is the covariance
627   * matrix (which must be positive-definite).
628   */
629  template<std::size_t _Dimen, typename _RealType = double>
630    class normal_mv_distribution
631    {
632      static_assert(std::is_floating_point<_RealType>::value,
633		    "template argument not a floating point type");
634      static_assert(_Dimen != 0, "dimension is zero");
635
636    public:
637      /** The type of the range of the distribution. */
638      typedef std::array<_RealType, _Dimen> result_type;
639      /** Parameter type. */
640      class param_type
641      {
642	static constexpr size_t _M_t_size = _Dimen * (_Dimen + 1) / 2;
643
644      public:
645	typedef normal_mv_distribution<_Dimen, _RealType> distribution_type;
646	friend class normal_mv_distribution<_Dimen, _RealType>;
647
648	param_type()
649	{
650	  std::fill(_M_mean.begin(), _M_mean.end(), _RealType(0));
651	  auto __it = _M_t.begin();
652	  for (size_t __i = 0; __i < _Dimen; ++__i)
653	    {
654	      std::fill_n(__it, __i, _RealType(0));
655	      __it += __i;
656	      *__it++ = _RealType(1);
657	    }
658	}
659
660	template<typename _ForwardIterator1, typename _ForwardIterator2>
661	  param_type(_ForwardIterator1 __meanbegin,
662		     _ForwardIterator1 __meanend,
663		     _ForwardIterator2 __varcovbegin,
664		     _ForwardIterator2 __varcovend)
665	{
666	  __glibcxx_function_requires(_ForwardIteratorConcept<
667				      _ForwardIterator1>)
668	  __glibcxx_function_requires(_ForwardIteratorConcept<
669				      _ForwardIterator2>)
670	  _GLIBCXX_DEBUG_ASSERT(std::distance(__meanbegin, __meanend)
671				<= _Dimen);
672	  const auto __dist = std::distance(__varcovbegin, __varcovend);
673	  _GLIBCXX_DEBUG_ASSERT(__dist == _Dimen * _Dimen
674				|| __dist == _Dimen * (_Dimen + 1) / 2
675				|| __dist == _Dimen);
676
677	  if (__dist == _Dimen * _Dimen)
678	    _M_init_full(__meanbegin, __meanend, __varcovbegin, __varcovend);
679	  else if (__dist == _Dimen * (_Dimen + 1) / 2)
680	    _M_init_lower(__meanbegin, __meanend, __varcovbegin, __varcovend);
681	  else
682	    {
683	      __glibcxx_assert(__dist == _Dimen);
684	      _M_init_diagonal(__meanbegin, __meanend,
685			       __varcovbegin, __varcovend);
686	    }
687	}
688
689	param_type(std::initializer_list<_RealType> __mean,
690		   std::initializer_list<_RealType> __varcov)
691	{
692	  _GLIBCXX_DEBUG_ASSERT(__mean.size() <= _Dimen);
693	  _GLIBCXX_DEBUG_ASSERT(__varcov.size() == _Dimen * _Dimen
694				|| __varcov.size() == _Dimen * (_Dimen + 1) / 2
695				|| __varcov.size() == _Dimen);
696
697	  if (__varcov.size() == _Dimen * _Dimen)
698	    _M_init_full(__mean.begin(), __mean.end(),
699			 __varcov.begin(), __varcov.end());
700	  else if (__varcov.size() == _Dimen * (_Dimen + 1) / 2)
701	    _M_init_lower(__mean.begin(), __mean.end(),
702			  __varcov.begin(), __varcov.end());
703	  else
704	    {
705	      __glibcxx_assert(__varcov.size() == _Dimen);
706	      _M_init_diagonal(__mean.begin(), __mean.end(),
707			       __varcov.begin(), __varcov.end());
708	    }
709	}
710
711	std::array<_RealType, _Dimen>
712	mean() const
713	{ return _M_mean; }
714
715	std::array<_RealType, _M_t_size>
716	varcov() const
717	{ return _M_t; }
718
719	friend bool
720	operator==(const param_type& __p1, const param_type& __p2)
721	{ return __p1._M_mean == __p2._M_mean && __p1._M_t == __p2._M_t; }
722
723	friend bool
724	operator!=(const param_type& __p1, const param_type& __p2)
725	{ return !(__p1 == __p2); }
726
727      private:
728	template <typename _InputIterator1, typename _InputIterator2>
729	  void _M_init_full(_InputIterator1 __meanbegin,
730			    _InputIterator1 __meanend,
731			    _InputIterator2 __varcovbegin,
732			    _InputIterator2 __varcovend);
733	template <typename _InputIterator1, typename _InputIterator2>
734	  void _M_init_lower(_InputIterator1 __meanbegin,
735			     _InputIterator1 __meanend,
736			     _InputIterator2 __varcovbegin,
737			     _InputIterator2 __varcovend);
738	template <typename _InputIterator1, typename _InputIterator2>
739	  void _M_init_diagonal(_InputIterator1 __meanbegin,
740				_InputIterator1 __meanend,
741				_InputIterator2 __varbegin,
742				_InputIterator2 __varend);
743
744	std::array<_RealType, _Dimen> _M_mean;
745	std::array<_RealType, _M_t_size> _M_t;
746      };
747
748    public:
749      normal_mv_distribution()
750      : _M_param(), _M_nd()
751      { }
752
753      template<typename _ForwardIterator1, typename _ForwardIterator2>
754	normal_mv_distribution(_ForwardIterator1 __meanbegin,
755			       _ForwardIterator1 __meanend,
756			       _ForwardIterator2 __varcovbegin,
757			       _ForwardIterator2 __varcovend)
758	: _M_param(__meanbegin, __meanend, __varcovbegin, __varcovend),
759	  _M_nd()
760	{ }
761
762      normal_mv_distribution(std::initializer_list<_RealType> __mean,
763			     std::initializer_list<_RealType> __varcov)
764      : _M_param(__mean, __varcov), _M_nd()
765      { }
766
767      explicit
768      normal_mv_distribution(const param_type& __p)
769      : _M_param(__p), _M_nd()
770      { }
771
772      /**
773       * @brief Resets the distribution state.
774       */
775      void
776      reset()
777      { _M_nd.reset(); }
778
779      /**
780       * @brief Returns the mean of the distribution.
781       */
782      result_type
783      mean() const
784      { return _M_param.mean(); }
785
786      /**
787       * @brief Returns the compact form of the variance/covariance
788       * matrix of the distribution.
789       */
790      std::array<_RealType, _Dimen * (_Dimen + 1) / 2>
791      varcov() const
792      { return _M_param.varcov(); }
793
794      /**
795       * @brief Returns the parameter set of the distribution.
796       */
797      param_type
798      param() const
799      { return _M_param; }
800
801      /**
802       * @brief Sets the parameter set of the distribution.
803       * @param __param The new parameter set of the distribution.
804       */
805      void
806      param(const param_type& __param)
807      { _M_param = __param; }
808
809      /**
810       * @brief Returns the greatest lower bound value of the distribution.
811       */
812      result_type
813      min() const
814      { result_type __res;
815	__res.fill(std::numeric_limits<_RealType>::lowest());
816	return __res; }
817
818      /**
819       * @brief Returns the least upper bound value of the distribution.
820       */
821      result_type
822      max() const
823      { result_type __res;
824	__res.fill(std::numeric_limits<_RealType>::max());
825	return __res; }
826
827      /**
828       * @brief Generating functions.
829       */
830      template<typename _UniformRandomNumberGenerator>
831	result_type
832	operator()(_UniformRandomNumberGenerator& __urng)
833	{ return this->operator()(__urng, _M_param); }
834
835      template<typename _UniformRandomNumberGenerator>
836	result_type
837	operator()(_UniformRandomNumberGenerator& __urng,
838		   const param_type& __p);
839
840      template<typename _ForwardIterator,
841	       typename _UniformRandomNumberGenerator>
842	void
843	__generate(_ForwardIterator __f, _ForwardIterator __t,
844		   _UniformRandomNumberGenerator& __urng)
845	{ return this->__generate_impl(__f, __t, __urng, _M_param); }
846
847      template<typename _ForwardIterator,
848	       typename _UniformRandomNumberGenerator>
849	void
850	__generate(_ForwardIterator __f, _ForwardIterator __t,
851		   _UniformRandomNumberGenerator& __urng,
852		   const param_type& __p)
853	{ return this->__generate_impl(__f, __t, __urng, __p); }
854
855      /**
856       * @brief Return true if two multi-variant normal distributions have
857       *        the same parameters and the sequences that would
858       *        be generated are equal.
859       */
860      template<size_t _Dimen1, typename _RealType1>
861	friend bool
862	operator==(const
863		   __gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>&
864		   __d1,
865		   const
866		   __gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>&
867		   __d2);
868
869      /**
870       * @brief Inserts a %normal_mv_distribution random number distribution
871       * @p __x into the output stream @p __os.
872       *
873       * @param __os An output stream.
874       * @param __x  A %normal_mv_distribution random number distribution.
875       *
876       * @returns The output stream with the state of @p __x inserted or in
877       * an error state.
878       */
879      template<size_t _Dimen1, typename _RealType1,
880	       typename _CharT, typename _Traits>
881	friend std::basic_ostream<_CharT, _Traits>&
882	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
883		   const
884		   __gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>&
885		   __x);
886
887      /**
888       * @brief Extracts a %normal_mv_distribution random number distribution
889       * @p __x from the input stream @p __is.
890       *
891       * @param __is An input stream.
892       * @param __x  A %normal_mv_distribution random number generator engine.
893       *
894       * @returns The input stream with @p __x extracted or in an error
895       *          state.
896       */
897      template<size_t _Dimen1, typename _RealType1,
898	       typename _CharT, typename _Traits>
899	friend std::basic_istream<_CharT, _Traits>&
900	operator>>(std::basic_istream<_CharT, _Traits>& __is,
901		   __gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>&
902		   __x);
903
904    private:
905      template<typename _ForwardIterator,
906	       typename _UniformRandomNumberGenerator>
907	void
908	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
909			_UniformRandomNumberGenerator& __urng,
910			const param_type& __p);
911
912      param_type _M_param;
913      std::normal_distribution<_RealType> _M_nd;
914  };
915
916  /**
917   * @brief Return true if two multi-variate normal distributions are
918   * different.
919   */
920  template<size_t _Dimen, typename _RealType>
921    inline bool
922    operator!=(const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
923	       __d1,
924	       const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
925	       __d2)
926    { return !(__d1 == __d2); }
927
928
929  /**
930   * @brief A Rice continuous distribution for random numbers.
931   *
932   * The formula for the Rice probability density function is
933   * @f[
934   *     p(x|\nu,\sigma) = \frac{x}{\sigma^2}
935   *                       \exp\left(-\frac{x^2+\nu^2}{2\sigma^2}\right)
936   *                       I_0\left(\frac{x \nu}{\sigma^2}\right)
937   * @f]
938   * where @f$I_0(z)@f$ is the modified Bessel function of the first kind
939   * of order 0 and @f$\nu >= 0@f$ and @f$\sigma > 0@f$.
940   *
941   * <table border=1 cellpadding=10 cellspacing=0>
942   * <caption align=top>Distribution Statistics</caption>
943   * <tr><td>Mean</td><td>@f$\sqrt{\pi/2}L_{1/2}(-\nu^2/2\sigma^2)@f$</td></tr>
944   * <tr><td>Variance</td><td>@f$2\sigma^2 + \nu^2
945   *                   + (\pi\sigma^2/2)L^2_{1/2}(-\nu^2/2\sigma^2)@f$</td></tr>
946   * <tr><td>Range</td><td>@f$[0, \infty)@f$</td></tr>
947   * </table>
948   * where @f$L_{1/2}(x)@f$ is the Laguerre polynomial of order 1/2.
949   */
950  template<typename _RealType = double>
951    class
952    rice_distribution
953    {
954      static_assert(std::is_floating_point<_RealType>::value,
955		    "template argument not a floating point type");
956    public:
957      /** The type of the range of the distribution. */
958      typedef _RealType result_type;
959
960      /** Parameter type. */
961      struct param_type
962      {
963	typedef rice_distribution<result_type> distribution_type;
964
965	param_type(result_type __nu_val = result_type(0),
966		   result_type __sigma_val = result_type(1))
967	: _M_nu(__nu_val), _M_sigma(__sigma_val)
968	{
969	  __glibcxx_assert(_M_nu >= result_type(0));
970	  __glibcxx_assert(_M_sigma > result_type(0));
971	}
972
973	result_type
974	nu() const
975	{ return _M_nu; }
976
977	result_type
978	sigma() const
979	{ return _M_sigma; }
980
981	friend bool
982	operator==(const param_type& __p1, const param_type& __p2)
983	{ return __p1._M_nu == __p2._M_nu && __p1._M_sigma == __p2._M_sigma; }
984
985	friend bool
986	operator!=(const param_type& __p1, const param_type& __p2)
987	{ return !(__p1 == __p2); }
988
989      private:
990	void _M_initialize();
991
992	result_type _M_nu;
993	result_type _M_sigma;
994      };
995
996      /**
997       * @brief Constructors.
998       */
999      explicit
1000      rice_distribution(result_type __nu_val = result_type(0),
1001			result_type __sigma_val = result_type(1))
1002      : _M_param(__nu_val, __sigma_val),
1003	_M_ndx(__nu_val, __sigma_val),
1004	_M_ndy(result_type(0), __sigma_val)
1005      { }
1006
1007      explicit
1008      rice_distribution(const param_type& __p)
1009      : _M_param(__p),
1010	_M_ndx(__p.nu(), __p.sigma()),
1011	_M_ndy(result_type(0), __p.sigma())
1012      { }
1013
1014      /**
1015       * @brief Resets the distribution state.
1016       */
1017      void
1018      reset()
1019      {
1020	_M_ndx.reset();
1021	_M_ndy.reset();
1022      }
1023
1024      /**
1025       * @brief Return the parameters of the distribution.
1026       */
1027      result_type
1028      nu() const
1029      { return _M_param.nu(); }
1030
1031      result_type
1032      sigma() const
1033      { return _M_param.sigma(); }
1034
1035      /**
1036       * @brief Returns the parameter set of the distribution.
1037       */
1038      param_type
1039      param() const
1040      { return _M_param; }
1041
1042      /**
1043       * @brief Sets the parameter set of the distribution.
1044       * @param __param The new parameter set of the distribution.
1045       */
1046      void
1047      param(const param_type& __param)
1048      { _M_param = __param; }
1049
1050      /**
1051       * @brief Returns the greatest lower bound value of the distribution.
1052       */
1053      result_type
1054      min() const
1055      { return result_type(0); }
1056
1057      /**
1058       * @brief Returns the least upper bound value of the distribution.
1059       */
1060      result_type
1061      max() const
1062      { return std::numeric_limits<result_type>::max(); }
1063
1064      /**
1065       * @brief Generating functions.
1066       */
1067      template<typename _UniformRandomNumberGenerator>
1068	result_type
1069	operator()(_UniformRandomNumberGenerator& __urng)
1070	{
1071	  result_type __x = this->_M_ndx(__urng);
1072	  result_type __y = this->_M_ndy(__urng);
1073#if _GLIBCXX_USE_C99_MATH_TR1
1074	  return std::hypot(__x, __y);
1075#else
1076	  return std::sqrt(__x * __x + __y * __y);
1077#endif
1078	}
1079
1080      template<typename _UniformRandomNumberGenerator>
1081	result_type
1082	operator()(_UniformRandomNumberGenerator& __urng,
1083		   const param_type& __p)
1084	{
1085	  typename std::normal_distribution<result_type>::param_type
1086	    __px(__p.nu(), __p.sigma()), __py(result_type(0), __p.sigma());
1087	  result_type __x = this->_M_ndx(__px, __urng);
1088	  result_type __y = this->_M_ndy(__py, __urng);
1089#if _GLIBCXX_USE_C99_MATH_TR1
1090	  return std::hypot(__x, __y);
1091#else
1092	  return std::sqrt(__x * __x + __y * __y);
1093#endif
1094	}
1095
1096      template<typename _ForwardIterator,
1097	       typename _UniformRandomNumberGenerator>
1098	void
1099	__generate(_ForwardIterator __f, _ForwardIterator __t,
1100		   _UniformRandomNumberGenerator& __urng)
1101	{ this->__generate(__f, __t, __urng, _M_param); }
1102
1103      template<typename _ForwardIterator,
1104	       typename _UniformRandomNumberGenerator>
1105	void
1106	__generate(_ForwardIterator __f, _ForwardIterator __t,
1107		   _UniformRandomNumberGenerator& __urng,
1108		   const param_type& __p)
1109	{ this->__generate_impl(__f, __t, __urng, __p); }
1110
1111      template<typename _UniformRandomNumberGenerator>
1112	void
1113	__generate(result_type* __f, result_type* __t,
1114		   _UniformRandomNumberGenerator& __urng,
1115		   const param_type& __p)
1116	{ this->__generate_impl(__f, __t, __urng, __p); }
1117
1118      /**
1119       * @brief Return true if two Rice distributions have
1120       *        the same parameters and the sequences that would
1121       *        be generated are equal.
1122       */
1123      friend bool
1124      operator==(const rice_distribution& __d1,
1125		 const rice_distribution& __d2)
1126      { return (__d1._M_param == __d2._M_param
1127		&& __d1._M_ndx == __d2._M_ndx
1128		&& __d1._M_ndy == __d2._M_ndy); }
1129
1130      /**
1131       * @brief Inserts a %rice_distribution random number distribution
1132       * @p __x into the output stream @p __os.
1133       *
1134       * @param __os An output stream.
1135       * @param __x  A %rice_distribution random number distribution.
1136       *
1137       * @returns The output stream with the state of @p __x inserted or in
1138       * an error state.
1139       */
1140      template<typename _RealType1, typename _CharT, typename _Traits>
1141	friend std::basic_ostream<_CharT, _Traits>&
1142	operator<<(std::basic_ostream<_CharT, _Traits>&,
1143		   const rice_distribution<_RealType1>&);
1144
1145      /**
1146       * @brief Extracts a %rice_distribution random number distribution
1147       * @p __x from the input stream @p __is.
1148       *
1149       * @param __is An input stream.
1150       * @param __x A %rice_distribution random number
1151       *            generator engine.
1152       *
1153       * @returns The input stream with @p __x extracted or in an error state.
1154       */
1155      template<typename _RealType1, typename _CharT, typename _Traits>
1156	friend std::basic_istream<_CharT, _Traits>&
1157	operator>>(std::basic_istream<_CharT, _Traits>&,
1158		   rice_distribution<_RealType1>&);
1159
1160    private:
1161      template<typename _ForwardIterator,
1162	       typename _UniformRandomNumberGenerator>
1163	void
1164	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1165			_UniformRandomNumberGenerator& __urng,
1166			const param_type& __p);
1167
1168      param_type _M_param;
1169
1170      std::normal_distribution<result_type> _M_ndx;
1171      std::normal_distribution<result_type> _M_ndy;
1172    };
1173
1174  /**
1175   * @brief Return true if two Rice distributions are not equal.
1176   */
1177  template<typename _RealType1>
1178    inline bool
1179    operator!=(const rice_distribution<_RealType1>& __d1,
1180	       const rice_distribution<_RealType1>& __d2)
1181    { return !(__d1 == __d2); }
1182
1183
1184  /**
1185   * @brief A Nakagami continuous distribution for random numbers.
1186   *
1187   * The formula for the Nakagami probability density function is
1188   * @f[
1189   *     p(x|\mu,\omega) = \frac{2\mu^\mu}{\Gamma(\mu)\omega^\mu}
1190   *                       x^{2\mu-1}e^{-\mu x / \omega}
1191   * @f]
1192   * where @f$\Gamma(z)@f$ is the gamma function and @f$\mu >= 0.5@f$
1193   * and @f$\omega > 0@f$.
1194   */
1195  template<typename _RealType = double>
1196    class
1197    nakagami_distribution
1198    {
1199      static_assert(std::is_floating_point<_RealType>::value,
1200		    "template argument not a floating point type");
1201
1202    public:
1203      /** The type of the range of the distribution. */
1204      typedef _RealType result_type;
1205
1206      /** Parameter type. */
1207      struct param_type
1208      {
1209	typedef nakagami_distribution<result_type> distribution_type;
1210
1211	param_type(result_type __mu_val = result_type(1),
1212		   result_type __omega_val = result_type(1))
1213	: _M_mu(__mu_val), _M_omega(__omega_val)
1214	{
1215	  __glibcxx_assert(_M_mu >= result_type(0.5L));
1216	  __glibcxx_assert(_M_omega > result_type(0));
1217	}
1218
1219	result_type
1220	mu() const
1221	{ return _M_mu; }
1222
1223	result_type
1224	omega() const
1225	{ return _M_omega; }
1226
1227	friend bool
1228	operator==(const param_type& __p1, const param_type& __p2)
1229	{ return __p1._M_mu == __p2._M_mu && __p1._M_omega == __p2._M_omega; }
1230
1231	friend bool
1232	operator!=(const param_type& __p1, const param_type& __p2)
1233	{ return !(__p1 == __p2); }
1234
1235      private:
1236	void _M_initialize();
1237
1238	result_type _M_mu;
1239	result_type _M_omega;
1240      };
1241
1242      /**
1243       * @brief Constructors.
1244       */
1245      explicit
1246      nakagami_distribution(result_type __mu_val = result_type(1),
1247			    result_type __omega_val = result_type(1))
1248      : _M_param(__mu_val, __omega_val),
1249	_M_gd(__mu_val, __omega_val / __mu_val)
1250      { }
1251
1252      explicit
1253      nakagami_distribution(const param_type& __p)
1254      : _M_param(__p),
1255	_M_gd(__p.mu(), __p.omega() / __p.mu())
1256      { }
1257
1258      /**
1259       * @brief Resets the distribution state.
1260       */
1261      void
1262      reset()
1263      { _M_gd.reset(); }
1264
1265      /**
1266       * @brief Return the parameters of the distribution.
1267       */
1268      result_type
1269      mu() const
1270      { return _M_param.mu(); }
1271
1272      result_type
1273      omega() const
1274      { return _M_param.omega(); }
1275
1276      /**
1277       * @brief Returns the parameter set of the distribution.
1278       */
1279      param_type
1280      param() const
1281      { return _M_param; }
1282
1283      /**
1284       * @brief Sets the parameter set of the distribution.
1285       * @param __param The new parameter set of the distribution.
1286       */
1287      void
1288      param(const param_type& __param)
1289      { _M_param = __param; }
1290
1291      /**
1292       * @brief Returns the greatest lower bound value of the distribution.
1293       */
1294      result_type
1295      min() const
1296      { return result_type(0); }
1297
1298      /**
1299       * @brief Returns the least upper bound value of the distribution.
1300       */
1301      result_type
1302      max() const
1303      { return std::numeric_limits<result_type>::max(); }
1304
1305      /**
1306       * @brief Generating functions.
1307       */
1308      template<typename _UniformRandomNumberGenerator>
1309	result_type
1310	operator()(_UniformRandomNumberGenerator& __urng)
1311	{ return std::sqrt(this->_M_gd(__urng)); }
1312
1313      template<typename _UniformRandomNumberGenerator>
1314	result_type
1315	operator()(_UniformRandomNumberGenerator& __urng,
1316		   const param_type& __p)
1317	{
1318	  typename std::gamma_distribution<result_type>::param_type
1319	    __pg(__p.mu(), __p.omega() / __p.mu());
1320	  return std::sqrt(this->_M_gd(__pg, __urng));
1321	}
1322
1323      template<typename _ForwardIterator,
1324	       typename _UniformRandomNumberGenerator>
1325	void
1326	__generate(_ForwardIterator __f, _ForwardIterator __t,
1327		   _UniformRandomNumberGenerator& __urng)
1328	{ this->__generate(__f, __t, __urng, _M_param); }
1329
1330      template<typename _ForwardIterator,
1331	       typename _UniformRandomNumberGenerator>
1332	void
1333	__generate(_ForwardIterator __f, _ForwardIterator __t,
1334		   _UniformRandomNumberGenerator& __urng,
1335		   const param_type& __p)
1336	{ this->__generate_impl(__f, __t, __urng, __p); }
1337
1338      template<typename _UniformRandomNumberGenerator>
1339	void
1340	__generate(result_type* __f, result_type* __t,
1341		   _UniformRandomNumberGenerator& __urng,
1342		   const param_type& __p)
1343	{ this->__generate_impl(__f, __t, __urng, __p); }
1344
1345      /**
1346       * @brief Return true if two Nakagami distributions have
1347       *        the same parameters and the sequences that would
1348       *        be generated are equal.
1349       */
1350      friend bool
1351      operator==(const nakagami_distribution& __d1,
1352		 const nakagami_distribution& __d2)
1353      { return (__d1._M_param == __d2._M_param
1354		&& __d1._M_gd == __d2._M_gd); }
1355
1356      /**
1357       * @brief Inserts a %nakagami_distribution random number distribution
1358       * @p __x into the output stream @p __os.
1359       *
1360       * @param __os An output stream.
1361       * @param __x  A %nakagami_distribution random number distribution.
1362       *
1363       * @returns The output stream with the state of @p __x inserted or in
1364       * an error state.
1365       */
1366      template<typename _RealType1, typename _CharT, typename _Traits>
1367	friend std::basic_ostream<_CharT, _Traits>&
1368	operator<<(std::basic_ostream<_CharT, _Traits>&,
1369		   const nakagami_distribution<_RealType1>&);
1370
1371      /**
1372       * @brief Extracts a %nakagami_distribution random number distribution
1373       * @p __x from the input stream @p __is.
1374       *
1375       * @param __is An input stream.
1376       * @param __x A %nakagami_distribution random number
1377       *            generator engine.
1378       *
1379       * @returns The input stream with @p __x extracted or in an error state.
1380       */
1381      template<typename _RealType1, typename _CharT, typename _Traits>
1382	friend std::basic_istream<_CharT, _Traits>&
1383	operator>>(std::basic_istream<_CharT, _Traits>&,
1384		   nakagami_distribution<_RealType1>&);
1385
1386    private:
1387      template<typename _ForwardIterator,
1388	       typename _UniformRandomNumberGenerator>
1389	void
1390	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1391			_UniformRandomNumberGenerator& __urng,
1392			const param_type& __p);
1393
1394      param_type _M_param;
1395
1396      std::gamma_distribution<result_type> _M_gd;
1397    };
1398
1399  /**
1400   * @brief Return true if two Nakagami distributions are not equal.
1401   */
1402  template<typename _RealType>
1403    inline bool
1404    operator!=(const nakagami_distribution<_RealType>& __d1,
1405	       const nakagami_distribution<_RealType>& __d2)
1406    { return !(__d1 == __d2); }
1407
1408
1409  /**
1410   * @brief A Pareto continuous distribution for random numbers.
1411   *
1412   * The formula for the Pareto cumulative probability function is
1413   * @f[
1414   *     P(x|\alpha,\mu) = 1 - \left(\frac{\mu}{x}\right)^\alpha
1415   * @f]
1416   * The formula for the Pareto probability density function is
1417   * @f[
1418   *     p(x|\alpha,\mu) = \frac{\alpha + 1}{\mu}
1419   *                       \left(\frac{\mu}{x}\right)^{\alpha + 1}
1420   * @f]
1421   * where @f$x >= \mu@f$ and @f$\mu > 0@f$, @f$\alpha > 0@f$.
1422   *
1423   * <table border=1 cellpadding=10 cellspacing=0>
1424   * <caption align=top>Distribution Statistics</caption>
1425   * <tr><td>Mean</td><td>@f$\alpha \mu / (\alpha - 1)@f$
1426   *              for @f$\alpha > 1@f$</td></tr>
1427   * <tr><td>Variance</td><td>@f$\alpha \mu^2 / [(\alpha - 1)^2(\alpha - 2)]@f$
1428   *              for @f$\alpha > 2@f$</td></tr>
1429   * <tr><td>Range</td><td>@f$[\mu, \infty)@f$</td></tr>
1430   * </table>
1431   */
1432  template<typename _RealType = double>
1433    class
1434    pareto_distribution
1435    {
1436      static_assert(std::is_floating_point<_RealType>::value,
1437		    "template argument not a floating point type");
1438
1439    public:
1440      /** The type of the range of the distribution. */
1441      typedef _RealType result_type;
1442
1443      /** Parameter type. */
1444      struct param_type
1445      {
1446	typedef pareto_distribution<result_type> distribution_type;
1447
1448	param_type(result_type __alpha_val = result_type(1),
1449		   result_type __mu_val = result_type(1))
1450	: _M_alpha(__alpha_val), _M_mu(__mu_val)
1451	{
1452	  __glibcxx_assert(_M_alpha > result_type(0));
1453	  __glibcxx_assert(_M_mu > result_type(0));
1454	}
1455
1456	result_type
1457	alpha() const
1458	{ return _M_alpha; }
1459
1460	result_type
1461	mu() const
1462	{ return _M_mu; }
1463
1464	friend bool
1465	operator==(const param_type& __p1, const param_type& __p2)
1466	{ return __p1._M_alpha == __p2._M_alpha && __p1._M_mu == __p2._M_mu; }
1467
1468	friend bool
1469	operator!=(const param_type& __p1, const param_type& __p2)
1470	{ return !(__p1 == __p2); }
1471
1472      private:
1473	void _M_initialize();
1474
1475	result_type _M_alpha;
1476	result_type _M_mu;
1477      };
1478
1479      /**
1480       * @brief Constructors.
1481       */
1482      explicit
1483      pareto_distribution(result_type __alpha_val = result_type(1),
1484			  result_type __mu_val = result_type(1))
1485      : _M_param(__alpha_val, __mu_val),
1486	_M_ud()
1487      { }
1488
1489      explicit
1490      pareto_distribution(const param_type& __p)
1491      : _M_param(__p),
1492	_M_ud()
1493      { }
1494
1495      /**
1496       * @brief Resets the distribution state.
1497       */
1498      void
1499      reset()
1500      {
1501	_M_ud.reset();
1502      }
1503
1504      /**
1505       * @brief Return the parameters of the distribution.
1506       */
1507      result_type
1508      alpha() const
1509      { return _M_param.alpha(); }
1510
1511      result_type
1512      mu() const
1513      { return _M_param.mu(); }
1514
1515      /**
1516       * @brief Returns the parameter set of the distribution.
1517       */
1518      param_type
1519      param() const
1520      { return _M_param; }
1521
1522      /**
1523       * @brief Sets the parameter set of the distribution.
1524       * @param __param The new parameter set of the distribution.
1525       */
1526      void
1527      param(const param_type& __param)
1528      { _M_param = __param; }
1529
1530      /**
1531       * @brief Returns the greatest lower bound value of the distribution.
1532       */
1533      result_type
1534      min() const
1535      { return this->mu(); }
1536
1537      /**
1538       * @brief Returns the least upper bound value of the distribution.
1539       */
1540      result_type
1541      max() const
1542      { return std::numeric_limits<result_type>::max(); }
1543
1544      /**
1545       * @brief Generating functions.
1546       */
1547      template<typename _UniformRandomNumberGenerator>
1548	result_type
1549	operator()(_UniformRandomNumberGenerator& __urng)
1550	{
1551	  return this->mu() * std::pow(this->_M_ud(__urng),
1552				       -result_type(1) / this->alpha());
1553	}
1554
1555      template<typename _UniformRandomNumberGenerator>
1556	result_type
1557	operator()(_UniformRandomNumberGenerator& __urng,
1558		   const param_type& __p)
1559	{
1560	  return __p.mu() * std::pow(this->_M_ud(__urng),
1561					   -result_type(1) / __p.alpha());
1562	}
1563
1564      template<typename _ForwardIterator,
1565	       typename _UniformRandomNumberGenerator>
1566	void
1567	__generate(_ForwardIterator __f, _ForwardIterator __t,
1568		   _UniformRandomNumberGenerator& __urng)
1569	{ this->__generate(__f, __t, __urng, _M_param); }
1570
1571      template<typename _ForwardIterator,
1572	       typename _UniformRandomNumberGenerator>
1573	void
1574	__generate(_ForwardIterator __f, _ForwardIterator __t,
1575		   _UniformRandomNumberGenerator& __urng,
1576		   const param_type& __p)
1577	{ this->__generate_impl(__f, __t, __urng, __p); }
1578
1579      template<typename _UniformRandomNumberGenerator>
1580	void
1581	__generate(result_type* __f, result_type* __t,
1582		   _UniformRandomNumberGenerator& __urng,
1583		   const param_type& __p)
1584	{ this->__generate_impl(__f, __t, __urng, __p); }
1585
1586      /**
1587       * @brief Return true if two Pareto distributions have
1588       *        the same parameters and the sequences that would
1589       *        be generated are equal.
1590       */
1591      friend bool
1592      operator==(const pareto_distribution& __d1,
1593		 const pareto_distribution& __d2)
1594      { return (__d1._M_param == __d2._M_param
1595		&& __d1._M_ud == __d2._M_ud); }
1596
1597      /**
1598       * @brief Inserts a %pareto_distribution random number distribution
1599       * @p __x into the output stream @p __os.
1600       *
1601       * @param __os An output stream.
1602       * @param __x  A %pareto_distribution random number distribution.
1603       *
1604       * @returns The output stream with the state of @p __x inserted or in
1605       * an error state.
1606       */
1607      template<typename _RealType1, typename _CharT, typename _Traits>
1608	friend std::basic_ostream<_CharT, _Traits>&
1609	operator<<(std::basic_ostream<_CharT, _Traits>&,
1610		   const pareto_distribution<_RealType1>&);
1611
1612      /**
1613       * @brief Extracts a %pareto_distribution random number distribution
1614       * @p __x from the input stream @p __is.
1615       *
1616       * @param __is An input stream.
1617       * @param __x A %pareto_distribution random number
1618       *            generator engine.
1619       *
1620       * @returns The input stream with @p __x extracted or in an error state.
1621       */
1622      template<typename _RealType1, typename _CharT, typename _Traits>
1623	friend std::basic_istream<_CharT, _Traits>&
1624	operator>>(std::basic_istream<_CharT, _Traits>&,
1625		   pareto_distribution<_RealType1>&);
1626
1627    private:
1628      template<typename _ForwardIterator,
1629	       typename _UniformRandomNumberGenerator>
1630	void
1631	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1632			_UniformRandomNumberGenerator& __urng,
1633			const param_type& __p);
1634
1635      param_type _M_param;
1636
1637      std::uniform_real_distribution<result_type> _M_ud;
1638    };
1639
1640  /**
1641   * @brief Return true if two Pareto distributions are not equal.
1642   */
1643  template<typename _RealType>
1644    inline bool
1645    operator!=(const pareto_distribution<_RealType>& __d1,
1646	       const pareto_distribution<_RealType>& __d2)
1647    { return !(__d1 == __d2); }
1648
1649
1650  /**
1651   * @brief A K continuous distribution for random numbers.
1652   *
1653   * The formula for the K probability density function is
1654   * @f[
1655   *     p(x|\lambda, \mu, \nu) = \frac{2}{x}
1656   *             \left(\frac{\lambda\nu x}{\mu}\right)^{\frac{\lambda + \nu}{2}}
1657   *             \frac{1}{\Gamma(\lambda)\Gamma(\nu)}
1658   *             K_{\nu - \lambda}\left(2\sqrt{\frac{\lambda\nu x}{\mu}}\right)
1659   * @f]
1660   * where @f$I_0(z)@f$ is the modified Bessel function of the second kind
1661   * of order @f$\nu - \lambda@f$ and @f$\lambda > 0@f$, @f$\mu > 0@f$
1662   * and @f$\nu > 0@f$.
1663   *
1664   * <table border=1 cellpadding=10 cellspacing=0>
1665   * <caption align=top>Distribution Statistics</caption>
1666   * <tr><td>Mean</td><td>@f$\mu@f$</td></tr>
1667   * <tr><td>Variance</td><td>@f$\mu^2\frac{\lambda + \nu + 1}{\lambda\nu}@f$</td></tr>
1668   * <tr><td>Range</td><td>@f$[0, \infty)@f$</td></tr>
1669   * </table>
1670   */
1671  template<typename _RealType = double>
1672    class
1673    k_distribution
1674    {
1675      static_assert(std::is_floating_point<_RealType>::value,
1676		    "template argument not a floating point type");
1677
1678    public:
1679      /** The type of the range of the distribution. */
1680      typedef _RealType result_type;
1681
1682      /** Parameter type. */
1683      struct param_type
1684      {
1685	typedef k_distribution<result_type> distribution_type;
1686
1687	param_type(result_type __lambda_val = result_type(1),
1688		   result_type __mu_val = result_type(1),
1689		   result_type __nu_val = result_type(1))
1690	: _M_lambda(__lambda_val), _M_mu(__mu_val), _M_nu(__nu_val)
1691	{
1692	  __glibcxx_assert(_M_lambda > result_type(0));
1693	  __glibcxx_assert(_M_mu > result_type(0));
1694	  __glibcxx_assert(_M_nu > result_type(0));
1695	}
1696
1697	result_type
1698	lambda() const
1699	{ return _M_lambda; }
1700
1701	result_type
1702	mu() const
1703	{ return _M_mu; }
1704
1705	result_type
1706	nu() const
1707	{ return _M_nu; }
1708
1709	friend bool
1710	operator==(const param_type& __p1, const param_type& __p2)
1711	{
1712	  return __p1._M_lambda == __p2._M_lambda
1713	      && __p1._M_mu == __p2._M_mu
1714	      && __p1._M_nu == __p2._M_nu;
1715	}
1716
1717	friend bool
1718	operator!=(const param_type& __p1, const param_type& __p2)
1719	{ return !(__p1 == __p2); }
1720
1721      private:
1722	void _M_initialize();
1723
1724	result_type _M_lambda;
1725	result_type _M_mu;
1726	result_type _M_nu;
1727      };
1728
1729      /**
1730       * @brief Constructors.
1731       */
1732      explicit
1733      k_distribution(result_type __lambda_val = result_type(1),
1734		     result_type __mu_val = result_type(1),
1735		     result_type __nu_val = result_type(1))
1736      : _M_param(__lambda_val, __mu_val, __nu_val),
1737	_M_gd1(__lambda_val, result_type(1) / __lambda_val),
1738	_M_gd2(__nu_val, __mu_val / __nu_val)
1739      { }
1740
1741      explicit
1742      k_distribution(const param_type& __p)
1743      : _M_param(__p),
1744	_M_gd1(__p.lambda(), result_type(1) / __p.lambda()),
1745	_M_gd2(__p.nu(), __p.mu() / __p.nu())
1746      { }
1747
1748      /**
1749       * @brief Resets the distribution state.
1750       */
1751      void
1752      reset()
1753      {
1754	_M_gd1.reset();
1755	_M_gd2.reset();
1756      }
1757
1758      /**
1759       * @brief Return the parameters of the distribution.
1760       */
1761      result_type
1762      lambda() const
1763      { return _M_param.lambda(); }
1764
1765      result_type
1766      mu() const
1767      { return _M_param.mu(); }
1768
1769      result_type
1770      nu() const
1771      { return _M_param.nu(); }
1772
1773      /**
1774       * @brief Returns the parameter set of the distribution.
1775       */
1776      param_type
1777      param() const
1778      { return _M_param; }
1779
1780      /**
1781       * @brief Sets the parameter set of the distribution.
1782       * @param __param The new parameter set of the distribution.
1783       */
1784      void
1785      param(const param_type& __param)
1786      { _M_param = __param; }
1787
1788      /**
1789       * @brief Returns the greatest lower bound value of the distribution.
1790       */
1791      result_type
1792      min() const
1793      { return result_type(0); }
1794
1795      /**
1796       * @brief Returns the least upper bound value of the distribution.
1797       */
1798      result_type
1799      max() const
1800      { return std::numeric_limits<result_type>::max(); }
1801
1802      /**
1803       * @brief Generating functions.
1804       */
1805      template<typename _UniformRandomNumberGenerator>
1806	result_type
1807	operator()(_UniformRandomNumberGenerator&);
1808
1809      template<typename _UniformRandomNumberGenerator>
1810	result_type
1811	operator()(_UniformRandomNumberGenerator&, const param_type&);
1812
1813      template<typename _ForwardIterator,
1814	       typename _UniformRandomNumberGenerator>
1815	void
1816	__generate(_ForwardIterator __f, _ForwardIterator __t,
1817		   _UniformRandomNumberGenerator& __urng)
1818	{ this->__generate(__f, __t, __urng, _M_param); }
1819
1820      template<typename _ForwardIterator,
1821	       typename _UniformRandomNumberGenerator>
1822	void
1823	__generate(_ForwardIterator __f, _ForwardIterator __t,
1824		   _UniformRandomNumberGenerator& __urng,
1825		   const param_type& __p)
1826	{ this->__generate_impl(__f, __t, __urng, __p); }
1827
1828      template<typename _UniformRandomNumberGenerator>
1829	void
1830	__generate(result_type* __f, result_type* __t,
1831		   _UniformRandomNumberGenerator& __urng,
1832		   const param_type& __p)
1833	{ this->__generate_impl(__f, __t, __urng, __p); }
1834
1835      /**
1836       * @brief Return true if two K distributions have
1837       *        the same parameters and the sequences that would
1838       *        be generated are equal.
1839       */
1840      friend bool
1841      operator==(const k_distribution& __d1,
1842		 const k_distribution& __d2)
1843      { return (__d1._M_param == __d2._M_param
1844		&& __d1._M_gd1 == __d2._M_gd1
1845		&& __d1._M_gd2 == __d2._M_gd2); }
1846
1847      /**
1848       * @brief Inserts a %k_distribution random number distribution
1849       * @p __x into the output stream @p __os.
1850       *
1851       * @param __os An output stream.
1852       * @param __x  A %k_distribution random number distribution.
1853       *
1854       * @returns The output stream with the state of @p __x inserted or in
1855       * an error state.
1856       */
1857      template<typename _RealType1, typename _CharT, typename _Traits>
1858	friend std::basic_ostream<_CharT, _Traits>&
1859	operator<<(std::basic_ostream<_CharT, _Traits>&,
1860		   const k_distribution<_RealType1>&);
1861
1862      /**
1863       * @brief Extracts a %k_distribution random number distribution
1864       * @p __x from the input stream @p __is.
1865       *
1866       * @param __is An input stream.
1867       * @param __x A %k_distribution random number
1868       *            generator engine.
1869       *
1870       * @returns The input stream with @p __x extracted or in an error state.
1871       */
1872      template<typename _RealType1, typename _CharT, typename _Traits>
1873	friend std::basic_istream<_CharT, _Traits>&
1874	operator>>(std::basic_istream<_CharT, _Traits>&,
1875		   k_distribution<_RealType1>&);
1876
1877    private:
1878      template<typename _ForwardIterator,
1879	       typename _UniformRandomNumberGenerator>
1880	void
1881	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1882			_UniformRandomNumberGenerator& __urng,
1883			const param_type& __p);
1884
1885      param_type _M_param;
1886
1887      std::gamma_distribution<result_type> _M_gd1;
1888      std::gamma_distribution<result_type> _M_gd2;
1889    };
1890
1891  /**
1892   * @brief Return true if two K distributions are not equal.
1893   */
1894  template<typename _RealType>
1895    inline bool
1896    operator!=(const k_distribution<_RealType>& __d1,
1897	       const k_distribution<_RealType>& __d2)
1898    { return !(__d1 == __d2); }
1899
1900
1901  /**
1902   * @brief An arcsine continuous distribution for random numbers.
1903   *
1904   * The formula for the arcsine probability density function is
1905   * @f[
1906   *     p(x|a,b) = \frac{1}{\pi \sqrt{(x - a)(b - x)}}
1907   * @f]
1908   * where @f$x >= a@f$ and @f$x <= b@f$.
1909   *
1910   * <table border=1 cellpadding=10 cellspacing=0>
1911   * <caption align=top>Distribution Statistics</caption>
1912   * <tr><td>Mean</td><td>@f$ (a + b) / 2 @f$</td></tr>
1913   * <tr><td>Variance</td><td>@f$ (b - a)^2 / 8 @f$</td></tr>
1914   * <tr><td>Range</td><td>@f$[a, b]@f$</td></tr>
1915   * </table>
1916   */
1917  template<typename _RealType = double>
1918    class
1919    arcsine_distribution
1920    {
1921      static_assert(std::is_floating_point<_RealType>::value,
1922		    "template argument not a floating point type");
1923
1924    public:
1925      /** The type of the range of the distribution. */
1926      typedef _RealType result_type;
1927
1928      /** Parameter type. */
1929      struct param_type
1930      {
1931	typedef arcsine_distribution<result_type> distribution_type;
1932
1933	param_type(result_type __a = result_type(0),
1934		   result_type __b = result_type(1))
1935	: _M_a(__a), _M_b(__b)
1936	{
1937	  __glibcxx_assert(_M_a <= _M_b);
1938	}
1939
1940	result_type
1941	a() const
1942	{ return _M_a; }
1943
1944	result_type
1945	b() const
1946	{ return _M_b; }
1947
1948	friend bool
1949	operator==(const param_type& __p1, const param_type& __p2)
1950	{ return __p1._M_a == __p2._M_a && __p1._M_b == __p2._M_b; }
1951
1952	friend bool
1953	operator!=(const param_type& __p1, const param_type& __p2)
1954	{ return !(__p1 == __p2); }
1955
1956      private:
1957	void _M_initialize();
1958
1959	result_type _M_a;
1960	result_type _M_b;
1961      };
1962
1963      /**
1964       * @brief Constructors.
1965       */
1966      explicit
1967      arcsine_distribution(result_type __a = result_type(0),
1968			   result_type __b = result_type(1))
1969      : _M_param(__a, __b),
1970	_M_ud(-1.5707963267948966192313216916397514L,
1971	      +1.5707963267948966192313216916397514L)
1972      { }
1973
1974      explicit
1975      arcsine_distribution(const param_type& __p)
1976      : _M_param(__p),
1977	_M_ud(-1.5707963267948966192313216916397514L,
1978	      +1.5707963267948966192313216916397514L)
1979      { }
1980
1981      /**
1982       * @brief Resets the distribution state.
1983       */
1984      void
1985      reset()
1986      { _M_ud.reset(); }
1987
1988      /**
1989       * @brief Return the parameters of the distribution.
1990       */
1991      result_type
1992      a() const
1993      { return _M_param.a(); }
1994
1995      result_type
1996      b() const
1997      { return _M_param.b(); }
1998
1999      /**
2000       * @brief Returns the parameter set of the distribution.
2001       */
2002      param_type
2003      param() const
2004      { return _M_param; }
2005
2006      /**
2007       * @brief Sets the parameter set of the distribution.
2008       * @param __param The new parameter set of the distribution.
2009       */
2010      void
2011      param(const param_type& __param)
2012      { _M_param = __param; }
2013
2014      /**
2015       * @brief Returns the greatest lower bound value of the distribution.
2016       */
2017      result_type
2018      min() const
2019      { return this->a(); }
2020
2021      /**
2022       * @brief Returns the least upper bound value of the distribution.
2023       */
2024      result_type
2025      max() const
2026      { return this->b(); }
2027
2028      /**
2029       * @brief Generating functions.
2030       */
2031      template<typename _UniformRandomNumberGenerator>
2032	result_type
2033	operator()(_UniformRandomNumberGenerator& __urng)
2034	{
2035	  result_type __x = std::sin(this->_M_ud(__urng));
2036	  return (__x * (this->b() - this->a())
2037		  + this->a() + this->b()) / result_type(2);
2038	}
2039
2040      template<typename _UniformRandomNumberGenerator>
2041	result_type
2042	operator()(_UniformRandomNumberGenerator& __urng,
2043		   const param_type& __p)
2044	{
2045	  result_type __x = std::sin(this->_M_ud(__urng));
2046	  return (__x * (__p.b() - __p.a())
2047		  + __p.a() + __p.b()) / result_type(2);
2048	}
2049
2050      template<typename _ForwardIterator,
2051	       typename _UniformRandomNumberGenerator>
2052	void
2053	__generate(_ForwardIterator __f, _ForwardIterator __t,
2054		   _UniformRandomNumberGenerator& __urng)
2055	{ this->__generate(__f, __t, __urng, _M_param); }
2056
2057      template<typename _ForwardIterator,
2058	       typename _UniformRandomNumberGenerator>
2059	void
2060	__generate(_ForwardIterator __f, _ForwardIterator __t,
2061		   _UniformRandomNumberGenerator& __urng,
2062		   const param_type& __p)
2063	{ this->__generate_impl(__f, __t, __urng, __p); }
2064
2065      template<typename _UniformRandomNumberGenerator>
2066	void
2067	__generate(result_type* __f, result_type* __t,
2068		   _UniformRandomNumberGenerator& __urng,
2069		   const param_type& __p)
2070	{ this->__generate_impl(__f, __t, __urng, __p); }
2071
2072      /**
2073       * @brief Return true if two arcsine distributions have
2074       *        the same parameters and the sequences that would
2075       *        be generated are equal.
2076       */
2077      friend bool
2078      operator==(const arcsine_distribution& __d1,
2079		 const arcsine_distribution& __d2)
2080      { return (__d1._M_param == __d2._M_param
2081		&& __d1._M_ud == __d2._M_ud); }
2082
2083      /**
2084       * @brief Inserts a %arcsine_distribution random number distribution
2085       * @p __x into the output stream @p __os.
2086       *
2087       * @param __os An output stream.
2088       * @param __x  A %arcsine_distribution random number distribution.
2089       *
2090       * @returns The output stream with the state of @p __x inserted or in
2091       * an error state.
2092       */
2093      template<typename _RealType1, typename _CharT, typename _Traits>
2094	friend std::basic_ostream<_CharT, _Traits>&
2095	operator<<(std::basic_ostream<_CharT, _Traits>&,
2096		   const arcsine_distribution<_RealType1>&);
2097
2098      /**
2099       * @brief Extracts a %arcsine_distribution random number distribution
2100       * @p __x from the input stream @p __is.
2101       *
2102       * @param __is An input stream.
2103       * @param __x A %arcsine_distribution random number
2104       *            generator engine.
2105       *
2106       * @returns The input stream with @p __x extracted or in an error state.
2107       */
2108      template<typename _RealType1, typename _CharT, typename _Traits>
2109	friend std::basic_istream<_CharT, _Traits>&
2110	operator>>(std::basic_istream<_CharT, _Traits>&,
2111		   arcsine_distribution<_RealType1>&);
2112
2113    private:
2114      template<typename _ForwardIterator,
2115	       typename _UniformRandomNumberGenerator>
2116	void
2117	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2118			_UniformRandomNumberGenerator& __urng,
2119			const param_type& __p);
2120
2121      param_type _M_param;
2122
2123      std::uniform_real_distribution<result_type> _M_ud;
2124    };
2125
2126  /**
2127   * @brief Return true if two arcsine distributions are not equal.
2128   */
2129  template<typename _RealType>
2130    inline bool
2131    operator!=(const arcsine_distribution<_RealType>& __d1,
2132	       const arcsine_distribution<_RealType>& __d2)
2133    { return !(__d1 == __d2); }
2134
2135
2136  /**
2137   * @brief A Hoyt continuous distribution for random numbers.
2138   *
2139   * The formula for the Hoyt probability density function is
2140   * @f[
2141   *     p(x|q,\omega) = \frac{(1 + q^2)x}{q\omega}
2142   *                     \exp\left(-\frac{(1 + q^2)^2 x^2}{4 q^2 \omega}\right)
2143   *                       I_0\left(\frac{(1 - q^4) x^2}{4 q^2 \omega}\right)
2144   * @f]
2145   * where @f$I_0(z)@f$ is the modified Bessel function of the first kind
2146   * of order 0 and @f$0 < q < 1@f$.
2147   *
2148   * <table border=1 cellpadding=10 cellspacing=0>
2149   * <caption align=top>Distribution Statistics</caption>
2150   * <tr><td>Mean</td><td>@f$ \sqrt{\frac{2}{\pi}} \sqrt{\frac{\omega}{1 + q^2}}
2151   *                       E(1 - q^2) @f$</td></tr>
2152   * <tr><td>Variance</td><td>@f$ \omega \left(1 - \frac{2E^2(1 - q^2)}
2153   *                                      {\pi (1 + q^2)}\right) @f$</td></tr>
2154   * <tr><td>Range</td><td>@f$[0, \infty)@f$</td></tr>
2155   * </table>
2156   * where @f$E(x)@f$ is the elliptic function of the second kind.
2157   */
2158  template<typename _RealType = double>
2159    class
2160    hoyt_distribution
2161    {
2162      static_assert(std::is_floating_point<_RealType>::value,
2163		    "template argument not a floating point type");
2164
2165    public:
2166      /** The type of the range of the distribution. */
2167      typedef _RealType result_type;
2168
2169      /** Parameter type. */
2170      struct param_type
2171      {
2172	typedef hoyt_distribution<result_type> distribution_type;
2173
2174	param_type(result_type __q = result_type(0.5L),
2175		   result_type __omega = result_type(1))
2176	: _M_q(__q), _M_omega(__omega)
2177	{
2178	  __glibcxx_assert(_M_q > result_type(0));
2179	  __glibcxx_assert(_M_q < result_type(1));
2180	}
2181
2182	result_type
2183	q() const
2184	{ return _M_q; }
2185
2186	result_type
2187	omega() const
2188	{ return _M_omega; }
2189
2190	friend bool
2191	operator==(const param_type& __p1, const param_type& __p2)
2192	{ return __p1._M_q == __p2._M_q && __p1._M_omega == __p2._M_omega; }
2193
2194	friend bool
2195	operator!=(const param_type& __p1, const param_type& __p2)
2196	{ return !(__p1 == __p2); }
2197
2198      private:
2199	void _M_initialize();
2200
2201	result_type _M_q;
2202	result_type _M_omega;
2203      };
2204
2205      /**
2206       * @brief Constructors.
2207       */
2208      explicit
2209      hoyt_distribution(result_type __q = result_type(0.5L),
2210			result_type __omega = result_type(1))
2211      : _M_param(__q, __omega),
2212	_M_ad(result_type(0.5L) * (result_type(1) + __q * __q),
2213	      result_type(0.5L) * (result_type(1) + __q * __q)
2214				/ (__q * __q)),
2215	_M_ed(result_type(1))
2216      { }
2217
2218      explicit
2219      hoyt_distribution(const param_type& __p)
2220      : _M_param(__p),
2221	_M_ad(result_type(0.5L) * (result_type(1) + __p.q() * __p.q()),
2222	      result_type(0.5L) * (result_type(1) + __p.q() * __p.q())
2223				/ (__p.q() * __p.q())),
2224	_M_ed(result_type(1))
2225      { }
2226
2227      /**
2228       * @brief Resets the distribution state.
2229       */
2230      void
2231      reset()
2232      {
2233	_M_ad.reset();
2234	_M_ed.reset();
2235      }
2236
2237      /**
2238       * @brief Return the parameters of the distribution.
2239       */
2240      result_type
2241      q() const
2242      { return _M_param.q(); }
2243
2244      result_type
2245      omega() const
2246      { return _M_param.omega(); }
2247
2248      /**
2249       * @brief Returns the parameter set of the distribution.
2250       */
2251      param_type
2252      param() const
2253      { return _M_param; }
2254
2255      /**
2256       * @brief Sets the parameter set of the distribution.
2257       * @param __param The new parameter set of the distribution.
2258       */
2259      void
2260      param(const param_type& __param)
2261      { _M_param = __param; }
2262
2263      /**
2264       * @brief Returns the greatest lower bound value of the distribution.
2265       */
2266      result_type
2267      min() const
2268      { return result_type(0); }
2269
2270      /**
2271       * @brief Returns the least upper bound value of the distribution.
2272       */
2273      result_type
2274      max() const
2275      { return std::numeric_limits<result_type>::max(); }
2276
2277      /**
2278       * @brief Generating functions.
2279       */
2280      template<typename _UniformRandomNumberGenerator>
2281	result_type
2282	operator()(_UniformRandomNumberGenerator& __urng);
2283
2284      template<typename _UniformRandomNumberGenerator>
2285	result_type
2286	operator()(_UniformRandomNumberGenerator& __urng,
2287		   const param_type& __p);
2288
2289      template<typename _ForwardIterator,
2290	       typename _UniformRandomNumberGenerator>
2291	void
2292	__generate(_ForwardIterator __f, _ForwardIterator __t,
2293		   _UniformRandomNumberGenerator& __urng)
2294	{ this->__generate(__f, __t, __urng, _M_param); }
2295
2296      template<typename _ForwardIterator,
2297	       typename _UniformRandomNumberGenerator>
2298	void
2299	__generate(_ForwardIterator __f, _ForwardIterator __t,
2300		   _UniformRandomNumberGenerator& __urng,
2301		   const param_type& __p)
2302	{ this->__generate_impl(__f, __t, __urng, __p); }
2303
2304      template<typename _UniformRandomNumberGenerator>
2305	void
2306	__generate(result_type* __f, result_type* __t,
2307		   _UniformRandomNumberGenerator& __urng,
2308		   const param_type& __p)
2309	{ this->__generate_impl(__f, __t, __urng, __p); }
2310
2311      /**
2312       * @brief Return true if two Hoyt distributions have
2313       *        the same parameters and the sequences that would
2314       *        be generated are equal.
2315       */
2316      friend bool
2317      operator==(const hoyt_distribution& __d1,
2318		 const hoyt_distribution& __d2)
2319      { return (__d1._M_param == __d2._M_param
2320		&& __d1._M_ad == __d2._M_ad
2321		&& __d1._M_ed == __d2._M_ed); }
2322
2323      /**
2324       * @brief Inserts a %hoyt_distribution random number distribution
2325       * @p __x into the output stream @p __os.
2326       *
2327       * @param __os An output stream.
2328       * @param __x  A %hoyt_distribution random number distribution.
2329       *
2330       * @returns The output stream with the state of @p __x inserted or in
2331       * an error state.
2332       */
2333      template<typename _RealType1, typename _CharT, typename _Traits>
2334	friend std::basic_ostream<_CharT, _Traits>&
2335	operator<<(std::basic_ostream<_CharT, _Traits>&,
2336		   const hoyt_distribution<_RealType1>&);
2337
2338      /**
2339       * @brief Extracts a %hoyt_distribution random number distribution
2340       * @p __x from the input stream @p __is.
2341       *
2342       * @param __is An input stream.
2343       * @param __x A %hoyt_distribution random number
2344       *            generator engine.
2345       *
2346       * @returns The input stream with @p __x extracted or in an error state.
2347       */
2348      template<typename _RealType1, typename _CharT, typename _Traits>
2349	friend std::basic_istream<_CharT, _Traits>&
2350	operator>>(std::basic_istream<_CharT, _Traits>&,
2351		   hoyt_distribution<_RealType1>&);
2352
2353    private:
2354      template<typename _ForwardIterator,
2355	       typename _UniformRandomNumberGenerator>
2356	void
2357	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2358			_UniformRandomNumberGenerator& __urng,
2359			const param_type& __p);
2360
2361      param_type _M_param;
2362
2363      __gnu_cxx::arcsine_distribution<result_type> _M_ad;
2364      std::exponential_distribution<result_type> _M_ed;
2365    };
2366
2367  /**
2368   * @brief Return true if two Hoyt distributions are not equal.
2369   */
2370  template<typename _RealType>
2371    inline bool
2372    operator!=(const hoyt_distribution<_RealType>& __d1,
2373	       const hoyt_distribution<_RealType>& __d2)
2374    { return !(__d1 == __d2); }
2375
2376
2377  /**
2378   * @brief A triangular distribution for random numbers.
2379   *
2380   * The formula for the triangular probability density function is
2381   * @f[
2382   *                  / 0                          for x < a
2383   *     p(x|a,b,c) = | \frac{2(x-a)}{(c-a)(b-a)}  for a <= x <= b
2384   *                  | \frac{2(c-x)}{(c-a)(c-b)}  for b < x <= c
2385   *                  \ 0                          for c < x
2386   * @f]
2387   *
2388   * <table border=1 cellpadding=10 cellspacing=0>
2389   * <caption align=top>Distribution Statistics</caption>
2390   * <tr><td>Mean</td><td>@f$ \frac{a+b+c}{2} @f$</td></tr>
2391   * <tr><td>Variance</td><td>@f$ \frac{a^2+b^2+c^2-ab-ac-bc}
2392   *                                   {18}@f$</td></tr>
2393   * <tr><td>Range</td><td>@f$[a, c]@f$</td></tr>
2394   * </table>
2395   */
2396  template<typename _RealType = double>
2397    class triangular_distribution
2398    {
2399      static_assert(std::is_floating_point<_RealType>::value,
2400		    "template argument not a floating point type");
2401
2402    public:
2403      /** The type of the range of the distribution. */
2404      typedef _RealType result_type;
2405
2406      /** Parameter type. */
2407      struct param_type
2408      {
2409	friend class triangular_distribution<_RealType>;
2410
2411	explicit
2412	param_type(_RealType __a = _RealType(0),
2413		   _RealType __b = _RealType(0.5),
2414		   _RealType __c = _RealType(1))
2415	: _M_a(__a), _M_b(__b), _M_c(__c)
2416	{
2417	  __glibcxx_assert(_M_a <= _M_b);
2418	  __glibcxx_assert(_M_b <= _M_c);
2419	  __glibcxx_assert(_M_a < _M_c);
2420
2421	  _M_r_ab = (_M_b - _M_a) / (_M_c - _M_a);
2422	  _M_f_ab_ac = (_M_b - _M_a) * (_M_c - _M_a);
2423	  _M_f_bc_ac = (_M_c - _M_b) * (_M_c - _M_a);
2424	}
2425
2426	_RealType
2427	a() const
2428	{ return _M_a; }
2429
2430	_RealType
2431	b() const
2432	{ return _M_b; }
2433
2434	_RealType
2435	c() const
2436	{ return _M_c; }
2437
2438	friend bool
2439	operator==(const param_type& __p1, const param_type& __p2)
2440	{
2441	  return (__p1._M_a == __p2._M_a && __p1._M_b == __p2._M_b
2442		  && __p1._M_c == __p2._M_c);
2443	}
2444
2445	friend bool
2446	operator!=(const param_type& __p1, const param_type& __p2)
2447	{ return !(__p1 == __p2); }
2448
2449      private:
2450
2451	_RealType _M_a;
2452	_RealType _M_b;
2453	_RealType _M_c;
2454	_RealType _M_r_ab;
2455	_RealType _M_f_ab_ac;
2456	_RealType _M_f_bc_ac;
2457      };
2458
2459      /**
2460       * @brief Constructs a triangle distribution with parameters
2461       * @f$ a @f$, @f$ b @f$ and @f$ c @f$.
2462       */
2463      explicit
2464      triangular_distribution(result_type __a = result_type(0),
2465			      result_type __b = result_type(0.5),
2466			      result_type __c = result_type(1))
2467      : _M_param(__a, __b, __c)
2468      { }
2469
2470      explicit
2471      triangular_distribution(const param_type& __p)
2472      : _M_param(__p)
2473      { }
2474
2475      /**
2476       * @brief Resets the distribution state.
2477       */
2478      void
2479      reset()
2480      { }
2481
2482      /**
2483       * @brief Returns the @f$ a @f$ of the distribution.
2484       */
2485      result_type
2486      a() const
2487      { return _M_param.a(); }
2488
2489      /**
2490       * @brief Returns the @f$ b @f$ of the distribution.
2491       */
2492      result_type
2493      b() const
2494      { return _M_param.b(); }
2495
2496      /**
2497       * @brief Returns the @f$ c @f$ of the distribution.
2498       */
2499      result_type
2500      c() const
2501      { return _M_param.c(); }
2502
2503      /**
2504       * @brief Returns the parameter set of the distribution.
2505       */
2506      param_type
2507      param() const
2508      { return _M_param; }
2509
2510      /**
2511       * @brief Sets the parameter set of the distribution.
2512       * @param __param The new parameter set of the distribution.
2513       */
2514      void
2515      param(const param_type& __param)
2516      { _M_param = __param; }
2517
2518      /**
2519       * @brief Returns the greatest lower bound value of the distribution.
2520       */
2521      result_type
2522      min() const
2523      { return _M_param._M_a; }
2524
2525      /**
2526       * @brief Returns the least upper bound value of the distribution.
2527       */
2528      result_type
2529      max() const
2530      { return _M_param._M_c; }
2531
2532      /**
2533       * @brief Generating functions.
2534       */
2535      template<typename _UniformRandomNumberGenerator>
2536	result_type
2537	operator()(_UniformRandomNumberGenerator& __urng)
2538	{ return this->operator()(__urng, _M_param); }
2539
2540      template<typename _UniformRandomNumberGenerator>
2541	result_type
2542	operator()(_UniformRandomNumberGenerator& __urng,
2543		   const param_type& __p)
2544	{
2545	  std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2546	    __aurng(__urng);
2547	  result_type __rnd = __aurng();
2548	  if (__rnd <= __p._M_r_ab)
2549	    return __p.a() + std::sqrt(__rnd * __p._M_f_ab_ac);
2550	  else
2551	    return __p.c() - std::sqrt((result_type(1) - __rnd)
2552				       * __p._M_f_bc_ac);
2553	}
2554
2555      template<typename _ForwardIterator,
2556	       typename _UniformRandomNumberGenerator>
2557	void
2558	__generate(_ForwardIterator __f, _ForwardIterator __t,
2559		   _UniformRandomNumberGenerator& __urng)
2560	{ this->__generate(__f, __t, __urng, _M_param); }
2561
2562      template<typename _ForwardIterator,
2563	       typename _UniformRandomNumberGenerator>
2564	void
2565	__generate(_ForwardIterator __f, _ForwardIterator __t,
2566		   _UniformRandomNumberGenerator& __urng,
2567		   const param_type& __p)
2568	{ this->__generate_impl(__f, __t, __urng, __p); }
2569
2570      template<typename _UniformRandomNumberGenerator>
2571	void
2572	__generate(result_type* __f, result_type* __t,
2573		   _UniformRandomNumberGenerator& __urng,
2574		   const param_type& __p)
2575	{ this->__generate_impl(__f, __t, __urng, __p); }
2576
2577      /**
2578       * @brief Return true if two triangle distributions have the same
2579       *        parameters and the sequences that would be generated
2580       *        are equal.
2581       */
2582      friend bool
2583      operator==(const triangular_distribution& __d1,
2584		 const triangular_distribution& __d2)
2585      { return __d1._M_param == __d2._M_param; }
2586
2587      /**
2588       * @brief Inserts a %triangular_distribution random number distribution
2589       * @p __x into the output stream @p __os.
2590       *
2591       * @param __os An output stream.
2592       * @param __x  A %triangular_distribution random number distribution.
2593       *
2594       * @returns The output stream with the state of @p __x inserted or in
2595       * an error state.
2596       */
2597      template<typename _RealType1, typename _CharT, typename _Traits>
2598	friend std::basic_ostream<_CharT, _Traits>&
2599	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2600		   const __gnu_cxx::triangular_distribution<_RealType1>& __x);
2601
2602      /**
2603       * @brief Extracts a %triangular_distribution random number distribution
2604       * @p __x from the input stream @p __is.
2605       *
2606       * @param __is An input stream.
2607       * @param __x  A %triangular_distribution random number generator engine.
2608       *
2609       * @returns The input stream with @p __x extracted or in an error state.
2610       */
2611      template<typename _RealType1, typename _CharT, typename _Traits>
2612	friend std::basic_istream<_CharT, _Traits>&
2613	operator>>(std::basic_istream<_CharT, _Traits>& __is,
2614		   __gnu_cxx::triangular_distribution<_RealType1>& __x);
2615
2616    private:
2617      template<typename _ForwardIterator,
2618	       typename _UniformRandomNumberGenerator>
2619	void
2620	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2621			_UniformRandomNumberGenerator& __urng,
2622			const param_type& __p);
2623
2624      param_type _M_param;
2625    };
2626
2627  /**
2628   * @brief Return true if two triangle distributions are different.
2629   */
2630  template<typename _RealType>
2631    inline bool
2632    operator!=(const __gnu_cxx::triangular_distribution<_RealType>& __d1,
2633	       const __gnu_cxx::triangular_distribution<_RealType>& __d2)
2634    { return !(__d1 == __d2); }
2635
2636
2637  /**
2638   * @brief A von Mises distribution for random numbers.
2639   *
2640   * The formula for the von Mises probability density function is
2641   * @f[
2642   *     p(x|\mu,\kappa) = \frac{e^{\kappa \cos(x-\mu)}}
2643   *                            {2\pi I_0(\kappa)}
2644   * @f]
2645   *
2646   * The generating functions use the method according to:
2647   *
2648   * D. J. Best and N. I. Fisher, 1979. "Efficient Simulation of the
2649   * von Mises Distribution", Journal of the Royal Statistical Society.
2650   * Series C (Applied Statistics), Vol. 28, No. 2, pp. 152-157.
2651   *
2652   * <table border=1 cellpadding=10 cellspacing=0>
2653   * <caption align=top>Distribution Statistics</caption>
2654   * <tr><td>Mean</td><td>@f$ \mu @f$</td></tr>
2655   * <tr><td>Variance</td><td>@f$ 1-I_1(\kappa)/I_0(\kappa) @f$</td></tr>
2656   * <tr><td>Range</td><td>@f$[-\pi, \pi]@f$</td></tr>
2657   * </table>
2658   */
2659  template<typename _RealType = double>
2660    class von_mises_distribution
2661    {
2662      static_assert(std::is_floating_point<_RealType>::value,
2663		    "template argument not a floating point type");
2664
2665    public:
2666      /** The type of the range of the distribution. */
2667      typedef _RealType result_type;
2668      /** Parameter type. */
2669      struct param_type
2670      {
2671	friend class von_mises_distribution<_RealType>;
2672
2673	explicit
2674	param_type(_RealType __mu = _RealType(0),
2675		   _RealType __kappa = _RealType(1))
2676	: _M_mu(__mu), _M_kappa(__kappa)
2677	{
2678	  const _RealType __pi = __gnu_cxx::__math_constants<_RealType>::__pi;
2679	  __glibcxx_assert(_M_mu >= -__pi && _M_mu <= __pi);
2680	  __glibcxx_assert(_M_kappa >= _RealType(0));
2681
2682	  auto __tau = std::sqrt(_RealType(4) * _M_kappa * _M_kappa
2683				 + _RealType(1)) + _RealType(1);
2684	  auto __rho = ((__tau - std::sqrt(_RealType(2) * __tau))
2685			/ (_RealType(2) * _M_kappa));
2686	  _M_r = (_RealType(1) + __rho * __rho) / (_RealType(2) * __rho);
2687	}
2688
2689	_RealType
2690	mu() const
2691	{ return _M_mu; }
2692
2693	_RealType
2694	kappa() const
2695	{ return _M_kappa; }
2696
2697	friend bool
2698	operator==(const param_type& __p1, const param_type& __p2)
2699	{ return __p1._M_mu == __p2._M_mu && __p1._M_kappa == __p2._M_kappa; }
2700
2701	friend bool
2702	operator!=(const param_type& __p1, const param_type& __p2)
2703	{ return !(__p1 == __p2); }
2704
2705      private:
2706	_RealType _M_mu;
2707	_RealType _M_kappa;
2708	_RealType _M_r;
2709      };
2710
2711      /**
2712       * @brief Constructs a von Mises distribution with parameters
2713       * @f$\mu@f$ and @f$\kappa@f$.
2714       */
2715      explicit
2716      von_mises_distribution(result_type __mu = result_type(0),
2717			     result_type __kappa = result_type(1))
2718	: _M_param(__mu, __kappa)
2719      { }
2720
2721      explicit
2722      von_mises_distribution(const param_type& __p)
2723      : _M_param(__p)
2724      { }
2725
2726      /**
2727       * @brief Resets the distribution state.
2728       */
2729      void
2730      reset()
2731      { }
2732
2733      /**
2734       * @brief Returns the @f$ \mu @f$ of the distribution.
2735       */
2736      result_type
2737      mu() const
2738      { return _M_param.mu(); }
2739
2740      /**
2741       * @brief Returns the @f$ \kappa @f$ of the distribution.
2742       */
2743      result_type
2744      kappa() const
2745      { return _M_param.kappa(); }
2746
2747      /**
2748       * @brief Returns the parameter set of the distribution.
2749       */
2750      param_type
2751      param() const
2752      { return _M_param; }
2753
2754      /**
2755       * @brief Sets the parameter set of the distribution.
2756       * @param __param The new parameter set of the distribution.
2757       */
2758      void
2759      param(const param_type& __param)
2760      { _M_param = __param; }
2761
2762      /**
2763       * @brief Returns the greatest lower bound value of the distribution.
2764       */
2765      result_type
2766      min() const
2767      {
2768	return -__gnu_cxx::__math_constants<result_type>::__pi;
2769      }
2770
2771      /**
2772       * @brief Returns the least upper bound value of the distribution.
2773       */
2774      result_type
2775      max() const
2776      {
2777	return __gnu_cxx::__math_constants<result_type>::__pi;
2778      }
2779
2780      /**
2781       * @brief Generating functions.
2782       */
2783      template<typename _UniformRandomNumberGenerator>
2784	result_type
2785	operator()(_UniformRandomNumberGenerator& __urng)
2786	{ return this->operator()(__urng, _M_param); }
2787
2788      template<typename _UniformRandomNumberGenerator>
2789	result_type
2790	operator()(_UniformRandomNumberGenerator& __urng,
2791		   const param_type& __p);
2792
2793      template<typename _ForwardIterator,
2794	       typename _UniformRandomNumberGenerator>
2795	void
2796	__generate(_ForwardIterator __f, _ForwardIterator __t,
2797		   _UniformRandomNumberGenerator& __urng)
2798	{ this->__generate(__f, __t, __urng, _M_param); }
2799
2800      template<typename _ForwardIterator,
2801	       typename _UniformRandomNumberGenerator>
2802	void
2803	__generate(_ForwardIterator __f, _ForwardIterator __t,
2804		   _UniformRandomNumberGenerator& __urng,
2805		   const param_type& __p)
2806	{ this->__generate_impl(__f, __t, __urng, __p); }
2807
2808      template<typename _UniformRandomNumberGenerator>
2809	void
2810	__generate(result_type* __f, result_type* __t,
2811		   _UniformRandomNumberGenerator& __urng,
2812		   const param_type& __p)
2813	{ this->__generate_impl(__f, __t, __urng, __p); }
2814
2815      /**
2816       * @brief Return true if two von Mises distributions have the same
2817       *        parameters and the sequences that would be generated
2818       *        are equal.
2819       */
2820      friend bool
2821      operator==(const von_mises_distribution& __d1,
2822		 const von_mises_distribution& __d2)
2823      { return __d1._M_param == __d2._M_param; }
2824
2825      /**
2826       * @brief Inserts a %von_mises_distribution random number distribution
2827       * @p __x into the output stream @p __os.
2828       *
2829       * @param __os An output stream.
2830       * @param __x  A %von_mises_distribution random number distribution.
2831       *
2832       * @returns The output stream with the state of @p __x inserted or in
2833       * an error state.
2834       */
2835      template<typename _RealType1, typename _CharT, typename _Traits>
2836	friend std::basic_ostream<_CharT, _Traits>&
2837	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2838		   const __gnu_cxx::von_mises_distribution<_RealType1>& __x);
2839
2840      /**
2841       * @brief Extracts a %von_mises_distribution random number distribution
2842       * @p __x from the input stream @p __is.
2843       *
2844       * @param __is An input stream.
2845       * @param __x  A %von_mises_distribution random number generator engine.
2846       *
2847       * @returns The input stream with @p __x extracted or in an error state.
2848       */
2849      template<typename _RealType1, typename _CharT, typename _Traits>
2850	friend std::basic_istream<_CharT, _Traits>&
2851	operator>>(std::basic_istream<_CharT, _Traits>& __is,
2852		   __gnu_cxx::von_mises_distribution<_RealType1>& __x);
2853
2854    private:
2855      template<typename _ForwardIterator,
2856	       typename _UniformRandomNumberGenerator>
2857	void
2858	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2859			_UniformRandomNumberGenerator& __urng,
2860			const param_type& __p);
2861
2862      param_type _M_param;
2863    };
2864
2865  /**
2866   * @brief Return true if two von Mises distributions are different.
2867   */
2868  template<typename _RealType>
2869    inline bool
2870    operator!=(const __gnu_cxx::von_mises_distribution<_RealType>& __d1,
2871	       const __gnu_cxx::von_mises_distribution<_RealType>& __d2)
2872    { return !(__d1 == __d2); }
2873
2874
2875  /**
2876   * @brief A discrete hypergeometric random number distribution.
2877   *
2878   * The hypergeometric distribution is a discrete probability distribution
2879   * that describes the probability of @p k successes in @p n draws @a without
2880   * replacement from a finite population of size @p N containing exactly @p K
2881   * successes.
2882   *
2883   * The formula for the hypergeometric probability density function is
2884   * @f[
2885   *   p(k|N,K,n) = \frac{\binom{K}{k} \binom{N-K}{n-k}}{\binom{N}{n}}
2886   * @f]
2887   * where @f$N@f$ is the total population of the distribution,
2888   * @f$K@f$ is the total population of the distribution.
2889   *
2890   * <table border=1 cellpadding=10 cellspacing=0>
2891   * <caption align=top>Distribution Statistics</caption>
2892   * <tr><td>Mean</td><td>@f$ n\frac{K}{N} @f$</td></tr>
2893   * <tr><td>Variance</td><td>@f$ n\frac{K}{N}\frac{N-K}{N}\frac{N-n}{N-1}
2894   *   @f$</td></tr>
2895   * <tr><td>Range</td><td>@f$[max(0, n+K-N), min(K, n)]@f$</td></tr>
2896   * </table>
2897   */
2898  template<typename _UIntType = unsigned int>
2899    class hypergeometric_distribution
2900    {
2901      static_assert(std::is_unsigned<_UIntType>::value, "template argument "
2902		    "substituting _UIntType not an unsigned integral type");
2903
2904    public:
2905      /** The type of the range of the distribution. */
2906      typedef _UIntType  result_type;
2907
2908      /** Parameter type. */
2909      struct param_type
2910      {
2911	typedef hypergeometric_distribution<_UIntType> distribution_type;
2912	friend class hypergeometric_distribution<_UIntType>;
2913
2914	explicit
2915	param_type(result_type __N = 10, result_type __K = 5,
2916		   result_type __n = 1)
2917	: _M_N{__N}, _M_K{__K}, _M_n{__n}
2918	{
2919	  __glibcxx_assert(_M_N >= _M_K);
2920	  __glibcxx_assert(_M_N >= _M_n);
2921	}
2922
2923	result_type
2924	total_size() const
2925	{ return _M_N; }
2926
2927	result_type
2928	successful_size() const
2929	{ return _M_K; }
2930
2931	result_type
2932	unsuccessful_size() const
2933	{ return _M_N - _M_K; }
2934
2935	result_type
2936	total_draws() const
2937	{ return _M_n; }
2938
2939	friend bool
2940	operator==(const param_type& __p1, const param_type& __p2)
2941	{ return (__p1._M_N == __p2._M_N)
2942	      && (__p1._M_K == __p2._M_K)
2943	      && (__p1._M_n == __p2._M_n); }
2944
2945	friend bool
2946	operator!=(const param_type& __p1, const param_type& __p2)
2947	{ return !(__p1 == __p2); }
2948
2949      private:
2950
2951	result_type _M_N;
2952	result_type _M_K;
2953	result_type _M_n;
2954      };
2955
2956      // constructors and member function
2957      explicit
2958      hypergeometric_distribution(result_type __N = 10, result_type __K = 5,
2959				  result_type __n = 1)
2960      : _M_param{__N, __K, __n}
2961      { }
2962
2963      explicit
2964      hypergeometric_distribution(const param_type& __p)
2965      : _M_param{__p}
2966      { }
2967
2968      /**
2969       * @brief Resets the distribution state.
2970       */
2971      void
2972      reset()
2973      { }
2974
2975      /**
2976       * @brief Returns the distribution parameter @p N,
2977       *	the total number of items.
2978       */
2979      result_type
2980      total_size() const
2981      { return this->_M_param.total_size(); }
2982
2983      /**
2984       * @brief Returns the distribution parameter @p K,
2985       *	the total number of successful items.
2986       */
2987      result_type
2988      successful_size() const
2989      { return this->_M_param.successful_size(); }
2990
2991      /**
2992       * @brief Returns the total number of unsuccessful items @f$ N - K @f$.
2993       */
2994      result_type
2995      unsuccessful_size() const
2996      { return this->_M_param.unsuccessful_size(); }
2997
2998      /**
2999       * @brief Returns the distribution parameter @p n,
3000       *	the total number of draws.
3001       */
3002      result_type
3003      total_draws() const
3004      { return this->_M_param.total_draws(); }
3005
3006      /**
3007       * @brief Returns the parameter set of the distribution.
3008       */
3009      param_type
3010      param() const
3011      { return this->_M_param; }
3012
3013      /**
3014       * @brief Sets the parameter set of the distribution.
3015       * @param __param The new parameter set of the distribution.
3016       */
3017      void
3018      param(const param_type& __param)
3019      { this->_M_param = __param; }
3020
3021      /**
3022       * @brief Returns the greatest lower bound value of the distribution.
3023       */
3024      result_type
3025      min() const
3026      {
3027	using _IntType = typename std::make_signed<result_type>::type;
3028	return static_cast<result_type>(std::max(static_cast<_IntType>(0),
3029				static_cast<_IntType>(this->total_draws()
3030						- this->unsuccessful_size())));
3031      }
3032
3033      /**
3034       * @brief Returns the least upper bound value of the distribution.
3035       */
3036      result_type
3037      max() const
3038      { return std::min(this->successful_size(), this->total_draws()); }
3039
3040      /**
3041       * @brief Generating functions.
3042       */
3043      template<typename _UniformRandomNumberGenerator>
3044	result_type
3045	operator()(_UniformRandomNumberGenerator& __urng)
3046	{ return this->operator()(__urng, this->_M_param); }
3047
3048      template<typename _UniformRandomNumberGenerator>
3049	result_type
3050	operator()(_UniformRandomNumberGenerator& __urng,
3051		   const param_type& __p);
3052
3053      template<typename _ForwardIterator,
3054	       typename _UniformRandomNumberGenerator>
3055	void
3056	__generate(_ForwardIterator __f, _ForwardIterator __t,
3057		   _UniformRandomNumberGenerator& __urng)
3058	{ this->__generate(__f, __t, __urng, this->_M_param); }
3059
3060      template<typename _ForwardIterator,
3061	       typename _UniformRandomNumberGenerator>
3062	void
3063	__generate(_ForwardIterator __f, _ForwardIterator __t,
3064		   _UniformRandomNumberGenerator& __urng,
3065		   const param_type& __p)
3066	{ this->__generate_impl(__f, __t, __urng, __p); }
3067
3068      template<typename _UniformRandomNumberGenerator>
3069	void
3070	__generate(result_type* __f, result_type* __t,
3071		   _UniformRandomNumberGenerator& __urng,
3072		   const param_type& __p)
3073	{ this->__generate_impl(__f, __t, __urng, __p); }
3074
3075       /**
3076	* @brief Return true if two hypergeometric distributions have the same
3077	*        parameters and the sequences that would be generated
3078	*        are equal.
3079	*/
3080      friend bool
3081      operator==(const hypergeometric_distribution& __d1,
3082		 const hypergeometric_distribution& __d2)
3083      { return __d1._M_param == __d2._M_param; }
3084
3085      /**
3086       * @brief Inserts a %hypergeometric_distribution random number
3087       * distribution @p __x into the output stream @p __os.
3088       *
3089       * @param __os An output stream.
3090       * @param __x  A %hypergeometric_distribution random number
3091       *             distribution.
3092       *
3093       * @returns The output stream with the state of @p __x inserted or in
3094       * an error state.
3095       */
3096      template<typename _UIntType1, typename _CharT, typename _Traits>
3097	friend std::basic_ostream<_CharT, _Traits>&
3098	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
3099		   const __gnu_cxx::hypergeometric_distribution<_UIntType1>&
3100		   __x);
3101
3102      /**
3103       * @brief Extracts a %hypergeometric_distribution random number
3104       * distribution @p __x from the input stream @p __is.
3105       *
3106       * @param __is An input stream.
3107       * @param __x  A %hypergeometric_distribution random number generator
3108       *             distribution.
3109       *
3110       * @returns The input stream with @p __x extracted or in an error
3111       *          state.
3112       */
3113      template<typename _UIntType1, typename _CharT, typename _Traits>
3114	friend std::basic_istream<_CharT, _Traits>&
3115	operator>>(std::basic_istream<_CharT, _Traits>& __is,
3116		   __gnu_cxx::hypergeometric_distribution<_UIntType1>& __x);
3117
3118    private:
3119
3120      template<typename _ForwardIterator,
3121	       typename _UniformRandomNumberGenerator>
3122	void
3123	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
3124			_UniformRandomNumberGenerator& __urng,
3125			const param_type& __p);
3126
3127      param_type _M_param;
3128    };
3129
3130  /**
3131   * @brief Return true if two hypergeometric distributions are different.
3132   */
3133  template<typename _UIntType>
3134    inline bool
3135    operator!=(const __gnu_cxx::hypergeometric_distribution<_UIntType>& __d1,
3136	       const __gnu_cxx::hypergeometric_distribution<_UIntType>& __d2)
3137    { return !(__d1 == __d2); }
3138
3139  /**
3140   * @brief A logistic continuous distribution for random numbers.
3141   *
3142   * The formula for the logistic probability density function is
3143   * @f[
3144   *     p(x|\a,\b) = \frac{e^{(x - a)/b}}{b[1 + e^{(x - a)/b}]^2}
3145   * @f]
3146   * where @f$b > 0@f$.
3147   *
3148   * The formula for the logistic probability function is
3149   * @f[
3150   *     cdf(x|\a,\b) = \frac{e^{(x - a)/b}}{1 + e^{(x - a)/b}}
3151   * @f]
3152   * where @f$b > 0@f$.
3153   *
3154   * <table border=1 cellpadding=10 cellspacing=0>
3155   * <caption align=top>Distribution Statistics</caption>
3156   * <tr><td>Mean</td><td>@f$a@f$</td></tr>
3157   * <tr><td>Variance</td><td>@f$b^2\pi^2/3@f$</td></tr>
3158   * <tr><td>Range</td><td>@f$[0, \infty)@f$</td></tr>
3159   * </table>
3160   */
3161  template<typename _RealType = double>
3162    class
3163    logistic_distribution
3164    {
3165      static_assert(std::is_floating_point<_RealType>::value,
3166		    "template argument not a floating point type");
3167
3168    public:
3169      /** The type of the range of the distribution. */
3170      typedef _RealType result_type;
3171
3172      /** Parameter type. */
3173      struct param_type
3174      {
3175	typedef logistic_distribution<result_type> distribution_type;
3176
3177	param_type(result_type __a = result_type(0),
3178		   result_type __b = result_type(1))
3179	: _M_a(__a), _M_b(__b)
3180	{
3181	  __glibcxx_assert(_M_b > result_type(0));
3182	}
3183
3184	result_type
3185	a() const
3186	{ return _M_a; }
3187
3188	result_type
3189	b() const
3190	{ return _M_b; }
3191
3192	friend bool
3193	operator==(const param_type& __p1, const param_type& __p2)
3194	{ return __p1._M_a == __p2._M_a && __p1._M_b == __p2._M_b; }
3195
3196	friend bool
3197	operator!=(const param_type& __p1, const param_type& __p2)
3198	{ return !(__p1 == __p2); }
3199
3200      private:
3201	void _M_initialize();
3202
3203	result_type _M_a;
3204	result_type _M_b;
3205      };
3206
3207      /**
3208       * @brief Constructors.
3209       */
3210      explicit
3211      logistic_distribution(result_type __a = result_type(0),
3212			    result_type __b = result_type(1))
3213      : _M_param(__a, __b)
3214      { }
3215
3216      explicit
3217      logistic_distribution(const param_type& __p)
3218      : _M_param(__p)
3219      { }
3220
3221      /**
3222       * @brief Resets the distribution state.
3223       */
3224      void
3225      reset()
3226      { }
3227
3228      /**
3229       * @brief Return the parameters of the distribution.
3230       */
3231      result_type
3232      a() const
3233      { return _M_param.a(); }
3234
3235      result_type
3236      b() const
3237      { return _M_param.b(); }
3238
3239      /**
3240       * @brief Returns the parameter set of the distribution.
3241       */
3242      param_type
3243      param() const
3244      { return _M_param; }
3245
3246      /**
3247       * @brief Sets the parameter set of the distribution.
3248       * @param __param The new parameter set of the distribution.
3249       */
3250      void
3251      param(const param_type& __param)
3252      { _M_param = __param; }
3253
3254      /**
3255       * @brief Returns the greatest lower bound value of the distribution.
3256       */
3257      result_type
3258      min() const
3259      { return -std::numeric_limits<result_type>::max(); }
3260
3261      /**
3262       * @brief Returns the least upper bound value of the distribution.
3263       */
3264      result_type
3265      max() const
3266      { return std::numeric_limits<result_type>::max(); }
3267
3268      /**
3269       * @brief Generating functions.
3270       */
3271      template<typename _UniformRandomNumberGenerator>
3272	result_type
3273	operator()(_UniformRandomNumberGenerator& __urng)
3274	{ return this->operator()(__urng, this->_M_param); }
3275
3276      template<typename _UniformRandomNumberGenerator>
3277	result_type
3278	operator()(_UniformRandomNumberGenerator&,
3279		   const param_type&);
3280
3281      template<typename _ForwardIterator,
3282	       typename _UniformRandomNumberGenerator>
3283	void
3284	__generate(_ForwardIterator __f, _ForwardIterator __t,
3285		   _UniformRandomNumberGenerator& __urng)
3286	{ this->__generate(__f, __t, __urng, this->param()); }
3287
3288      template<typename _ForwardIterator,
3289	       typename _UniformRandomNumberGenerator>
3290	void
3291	__generate(_ForwardIterator __f, _ForwardIterator __t,
3292		   _UniformRandomNumberGenerator& __urng,
3293		   const param_type& __p)
3294	{ this->__generate_impl(__f, __t, __urng, __p); }
3295
3296      template<typename _UniformRandomNumberGenerator>
3297	void
3298	__generate(result_type* __f, result_type* __t,
3299		   _UniformRandomNumberGenerator& __urng,
3300		   const param_type& __p)
3301	{ this->__generate_impl(__f, __t, __urng, __p); }
3302
3303      /**
3304       * @brief Return true if two logistic distributions have
3305       *        the same parameters and the sequences that would
3306       *        be generated are equal.
3307       */
3308      template<typename _RealType1>
3309	friend bool
3310	operator==(const logistic_distribution<_RealType1>& __d1,
3311		   const logistic_distribution<_RealType1>& __d2)
3312	{ return __d1.param() == __d2.param(); }
3313
3314      /**
3315       * @brief Inserts a %logistic_distribution random number distribution
3316       * @p __x into the output stream @p __os.
3317       *
3318       * @param __os An output stream.
3319       * @param __x  A %logistic_distribution random number distribution.
3320       *
3321       * @returns The output stream with the state of @p __x inserted or in
3322       * an error state.
3323       */
3324      template<typename _RealType1, typename _CharT, typename _Traits>
3325	friend std::basic_ostream<_CharT, _Traits>&
3326	operator<<(std::basic_ostream<_CharT, _Traits>&,
3327		   const logistic_distribution<_RealType1>&);
3328
3329      /**
3330       * @brief Extracts a %logistic_distribution random number distribution
3331       * @p __x from the input stream @p __is.
3332       *
3333       * @param __is An input stream.
3334       * @param __x A %logistic_distribution random number
3335       *            generator engine.
3336       *
3337       * @returns The input stream with @p __x extracted or in an error state.
3338       */
3339      template<typename _RealType1, typename _CharT, typename _Traits>
3340	friend std::basic_istream<_CharT, _Traits>&
3341	operator>>(std::basic_istream<_CharT, _Traits>&,
3342		   logistic_distribution<_RealType1>&);
3343
3344    private:
3345      template<typename _ForwardIterator,
3346	       typename _UniformRandomNumberGenerator>
3347	void
3348	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
3349			_UniformRandomNumberGenerator& __urng,
3350			const param_type& __p);
3351
3352      param_type _M_param;
3353    };
3354
3355  /**
3356   * @brief Return true if two logistic distributions are not equal.
3357   */
3358  template<typename _RealType1>
3359    inline bool
3360    operator!=(const logistic_distribution<_RealType1>& __d1,
3361	       const logistic_distribution<_RealType1>& __d2)
3362    { return !(__d1 == __d2); }
3363
3364
3365  /**
3366   * @brief A distribution for random coordinates on a unit sphere.
3367   *
3368   * The method used in the generation function is attributed by Donald Knuth
3369   * to G. W. Brown, Modern Mathematics for the Engineer (1956).
3370   */
3371  template<std::size_t _Dimen, typename _RealType = double>
3372    class uniform_on_sphere_distribution
3373    {
3374      static_assert(std::is_floating_point<_RealType>::value,
3375		    "template argument not a floating point type");
3376      static_assert(_Dimen != 0, "dimension is zero");
3377
3378    public:
3379      /** The type of the range of the distribution. */
3380      typedef std::array<_RealType, _Dimen> result_type;
3381
3382      /** Parameter type. */
3383      struct param_type
3384      {
3385	explicit
3386	param_type()
3387	{ }
3388
3389	friend bool
3390	operator==(const param_type&, const param_type&)
3391	{ return true; }
3392
3393	friend bool
3394	operator!=(const param_type&, const param_type&)
3395	{ return false; }
3396      };
3397
3398      /**
3399       * @brief Constructs a uniform on sphere distribution.
3400       */
3401      explicit
3402      uniform_on_sphere_distribution()
3403      : _M_param(), _M_nd()
3404      { }
3405
3406      explicit
3407      uniform_on_sphere_distribution(const param_type& __p)
3408      : _M_param(__p), _M_nd()
3409      { }
3410
3411      /**
3412       * @brief Resets the distribution state.
3413       */
3414      void
3415      reset()
3416      { _M_nd.reset(); }
3417
3418      /**
3419       * @brief Returns the parameter set of the distribution.
3420       */
3421      param_type
3422      param() const
3423      { return _M_param; }
3424
3425      /**
3426       * @brief Sets the parameter set of the distribution.
3427       * @param __param The new parameter set of the distribution.
3428       */
3429      void
3430      param(const param_type& __param)
3431      { _M_param = __param; }
3432
3433      /**
3434       * @brief Returns the greatest lower bound value of the distribution.
3435       * This function makes no sense for this distribution.
3436       */
3437      result_type
3438      min() const
3439      {
3440	result_type __res;
3441	__res.fill(0);
3442	return __res;
3443      }
3444
3445      /**
3446       * @brief Returns the least upper bound value of the distribution.
3447       * This function makes no sense for this distribution.
3448       */
3449      result_type
3450      max() const
3451      {
3452	result_type __res;
3453	__res.fill(0);
3454	return __res;
3455      }
3456
3457      /**
3458       * @brief Generating functions.
3459       */
3460      template<typename _UniformRandomNumberGenerator>
3461	result_type
3462	operator()(_UniformRandomNumberGenerator& __urng)
3463	{ return this->operator()(__urng, _M_param); }
3464
3465      template<typename _UniformRandomNumberGenerator>
3466	result_type
3467	operator()(_UniformRandomNumberGenerator& __urng,
3468		   const param_type& __p);
3469
3470      template<typename _ForwardIterator,
3471	       typename _UniformRandomNumberGenerator>
3472	void
3473	__generate(_ForwardIterator __f, _ForwardIterator __t,
3474		   _UniformRandomNumberGenerator& __urng)
3475	{ this->__generate(__f, __t, __urng, this->param()); }
3476
3477      template<typename _ForwardIterator,
3478	       typename _UniformRandomNumberGenerator>
3479	void
3480	__generate(_ForwardIterator __f, _ForwardIterator __t,
3481		   _UniformRandomNumberGenerator& __urng,
3482		   const param_type& __p)
3483	{ this->__generate_impl(__f, __t, __urng, __p); }
3484
3485      template<typename _UniformRandomNumberGenerator>
3486	void
3487	__generate(result_type* __f, result_type* __t,
3488		   _UniformRandomNumberGenerator& __urng,
3489		   const param_type& __p)
3490	{ this->__generate_impl(__f, __t, __urng, __p); }
3491
3492      /**
3493       * @brief Return true if two uniform on sphere distributions have
3494       *        the same parameters and the sequences that would be
3495       *        generated are equal.
3496       */
3497      friend bool
3498      operator==(const uniform_on_sphere_distribution& __d1,
3499		 const uniform_on_sphere_distribution& __d2)
3500      { return __d1._M_nd == __d2._M_nd; }
3501
3502      /**
3503       * @brief Inserts a %uniform_on_sphere_distribution random number
3504       *        distribution @p __x into the output stream @p __os.
3505       *
3506       * @param __os An output stream.
3507       * @param __x  A %uniform_on_sphere_distribution random number
3508       *             distribution.
3509       *
3510       * @returns The output stream with the state of @p __x inserted or in
3511       * an error state.
3512       */
3513      template<size_t _Dimen1, typename _RealType1, typename _CharT,
3514	       typename _Traits>
3515	friend std::basic_ostream<_CharT, _Traits>&
3516	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
3517		   const __gnu_cxx::uniform_on_sphere_distribution<_Dimen1,
3518								   _RealType1>&
3519		   __x);
3520
3521      /**
3522       * @brief Extracts a %uniform_on_sphere_distribution random number
3523       *        distribution
3524       * @p __x from the input stream @p __is.
3525       *
3526       * @param __is An input stream.
3527       * @param __x  A %uniform_on_sphere_distribution random number
3528       *             generator engine.
3529       *
3530       * @returns The input stream with @p __x extracted or in an error state.
3531       */
3532      template<std::size_t _Dimen1, typename _RealType1, typename _CharT,
3533	       typename _Traits>
3534	friend std::basic_istream<_CharT, _Traits>&
3535	operator>>(std::basic_istream<_CharT, _Traits>& __is,
3536		   __gnu_cxx::uniform_on_sphere_distribution<_Dimen1,
3537							     _RealType1>& __x);
3538
3539    private:
3540      template<typename _ForwardIterator,
3541	       typename _UniformRandomNumberGenerator>
3542	void
3543	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
3544			_UniformRandomNumberGenerator& __urng,
3545			const param_type& __p);
3546
3547      param_type _M_param;
3548      std::normal_distribution<_RealType> _M_nd;
3549    };
3550
3551  /**
3552   * @brief Return true if two uniform on sphere distributions are different.
3553   */
3554  template<std::size_t _Dimen, typename _RealType>
3555    inline bool
3556    operator!=(const __gnu_cxx::uniform_on_sphere_distribution<_Dimen,
3557	       _RealType>& __d1,
3558	       const __gnu_cxx::uniform_on_sphere_distribution<_Dimen,
3559	       _RealType>& __d2)
3560    { return !(__d1 == __d2); }
3561
3562
3563  /**
3564   * @brief A distribution for random coordinates inside a unit sphere.
3565   */
3566  template<std::size_t _Dimen, typename _RealType = double>
3567    class uniform_inside_sphere_distribution
3568    {
3569      static_assert(std::is_floating_point<_RealType>::value,
3570		    "template argument not a floating point type");
3571      static_assert(_Dimen != 0, "dimension is zero");
3572
3573    public:
3574      /** The type of the range of the distribution. */
3575      using result_type = std::array<_RealType, _Dimen>;
3576
3577      /** Parameter type. */
3578      struct param_type
3579      {
3580	using distribution_type
3581	  = uniform_inside_sphere_distribution<_Dimen, _RealType>;
3582	friend class uniform_inside_sphere_distribution<_Dimen, _RealType>;
3583
3584	explicit
3585	param_type(_RealType __radius = _RealType(1))
3586	: _M_radius(__radius)
3587	{
3588	  __glibcxx_assert(_M_radius > _RealType(0));
3589	}
3590
3591	_RealType
3592	radius() const
3593	{ return _M_radius; }
3594
3595	friend bool
3596	operator==(const param_type& __p1, const param_type& __p2)
3597	{ return __p1._M_radius == __p2._M_radius; }
3598
3599	friend bool
3600	operator!=(const param_type& __p1, const param_type& __p2)
3601	{ return !(__p1 == __p2); }
3602
3603      private:
3604	_RealType _M_radius;
3605      };
3606
3607      /**
3608       * @brief Constructors.
3609       */
3610      explicit
3611      uniform_inside_sphere_distribution(_RealType __radius = _RealType(1))
3612      : _M_param(__radius), _M_uosd()
3613      { }
3614
3615      explicit
3616      uniform_inside_sphere_distribution(const param_type& __p)
3617      : _M_param(__p), _M_uosd()
3618      { }
3619
3620      /**
3621       * @brief Resets the distribution state.
3622       */
3623      void
3624      reset()
3625      { _M_uosd.reset(); }
3626
3627      /**
3628       * @brief Returns the @f$radius@f$ of the distribution.
3629       */
3630      _RealType
3631      radius() const
3632      { return _M_param.radius(); }
3633
3634      /**
3635       * @brief Returns the parameter set of the distribution.
3636       */
3637      param_type
3638      param() const
3639      { return _M_param; }
3640
3641      /**
3642       * @brief Sets the parameter set of the distribution.
3643       * @param __param The new parameter set of the distribution.
3644       */
3645      void
3646      param(const param_type& __param)
3647      { _M_param = __param; }
3648
3649      /**
3650       * @brief Returns the greatest lower bound value of the distribution.
3651       * This function makes no sense for this distribution.
3652       */
3653      result_type
3654      min() const
3655      {
3656	result_type __res;
3657	__res.fill(0);
3658	return __res;
3659      }
3660
3661      /**
3662       * @brief Returns the least upper bound value of the distribution.
3663       * This function makes no sense for this distribution.
3664       */
3665      result_type
3666      max() const
3667      {
3668	result_type __res;
3669	__res.fill(0);
3670	return __res;
3671      }
3672
3673      /**
3674       * @brief Generating functions.
3675       */
3676      template<typename _UniformRandomNumberGenerator>
3677	result_type
3678	operator()(_UniformRandomNumberGenerator& __urng)
3679	{ return this->operator()(__urng, _M_param); }
3680
3681      template<typename _UniformRandomNumberGenerator>
3682	result_type
3683	operator()(_UniformRandomNumberGenerator& __urng,
3684		   const param_type& __p);
3685
3686      template<typename _ForwardIterator,
3687	       typename _UniformRandomNumberGenerator>
3688	void
3689	__generate(_ForwardIterator __f, _ForwardIterator __t,
3690		   _UniformRandomNumberGenerator& __urng)
3691	{ this->__generate(__f, __t, __urng, this->param()); }
3692
3693      template<typename _ForwardIterator,
3694	       typename _UniformRandomNumberGenerator>
3695	void
3696	__generate(_ForwardIterator __f, _ForwardIterator __t,
3697		   _UniformRandomNumberGenerator& __urng,
3698		   const param_type& __p)
3699	{ this->__generate_impl(__f, __t, __urng, __p); }
3700
3701      template<typename _UniformRandomNumberGenerator>
3702	void
3703	__generate(result_type* __f, result_type* __t,
3704		   _UniformRandomNumberGenerator& __urng,
3705		   const param_type& __p)
3706	{ this->__generate_impl(__f, __t, __urng, __p); }
3707
3708      /**
3709       * @brief Return true if two uniform on sphere distributions have
3710       *        the same parameters and the sequences that would be
3711       *        generated are equal.
3712       */
3713      friend bool
3714      operator==(const uniform_inside_sphere_distribution& __d1,
3715		 const uniform_inside_sphere_distribution& __d2)
3716      { return __d1._M_param == __d2._M_param && __d1._M_uosd == __d2._M_uosd; }
3717
3718      /**
3719       * @brief Inserts a %uniform_inside_sphere_distribution random number
3720       *        distribution @p __x into the output stream @p __os.
3721       *
3722       * @param __os An output stream.
3723       * @param __x  A %uniform_inside_sphere_distribution random number
3724       *             distribution.
3725       *
3726       * @returns The output stream with the state of @p __x inserted or in
3727       * an error state.
3728       */
3729      template<size_t _Dimen1, typename _RealType1, typename _CharT,
3730	       typename _Traits>
3731	friend std::basic_ostream<_CharT, _Traits>&
3732	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
3733		   const __gnu_cxx::uniform_inside_sphere_distribution<_Dimen1,
3734								   _RealType1>&
3735		   );
3736
3737      /**
3738       * @brief Extracts a %uniform_inside_sphere_distribution random number
3739       *        distribution
3740       * @p __x from the input stream @p __is.
3741       *
3742       * @param __is An input stream.
3743       * @param __x  A %uniform_inside_sphere_distribution random number
3744       *             generator engine.
3745       *
3746       * @returns The input stream with @p __x extracted or in an error state.
3747       */
3748      template<std::size_t _Dimen1, typename _RealType1, typename _CharT,
3749	       typename _Traits>
3750	friend std::basic_istream<_CharT, _Traits>&
3751	operator>>(std::basic_istream<_CharT, _Traits>& __is,
3752		   __gnu_cxx::uniform_inside_sphere_distribution<_Dimen1,
3753								 _RealType1>&);
3754
3755    private:
3756      template<typename _ForwardIterator,
3757	       typename _UniformRandomNumberGenerator>
3758	void
3759	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
3760			_UniformRandomNumberGenerator& __urng,
3761			const param_type& __p);
3762
3763      param_type _M_param;
3764      uniform_on_sphere_distribution<_Dimen, _RealType> _M_uosd;
3765    };
3766
3767  /**
3768   * @brief Return true if two uniform on sphere distributions are different.
3769   */
3770  template<std::size_t _Dimen, typename _RealType>
3771    inline bool
3772    operator!=(const __gnu_cxx::uniform_inside_sphere_distribution<_Dimen,
3773	       _RealType>& __d1,
3774	       const __gnu_cxx::uniform_inside_sphere_distribution<_Dimen,
3775	       _RealType>& __d2)
3776    { return !(__d1 == __d2); }
3777
3778_GLIBCXX_END_NAMESPACE_VERSION
3779} // namespace __gnu_cxx
3780
3781#include <ext/opt_random.h>
3782#include <ext/random.tcc>
3783
3784#endif // _GLIBCXX_USE_C99_STDINT_TR1 && UINT32_C
3785
3786#endif // C++11
3787
3788#endif // _EXT_RANDOM
3789