random revision 1.1.1.1
1// Random number extensions -*- C++ -*-
2
3// Copyright (C) 2012-2013 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 <array>
40#include <ext/cmath>
41#ifdef __SSE2__
42# include <x86intrin.h>
43#endif
44
45#ifdef _GLIBCXX_USE_C99_STDINT_TR1
46
47namespace __gnu_cxx _GLIBCXX_VISIBILITY(default)
48{
49_GLIBCXX_BEGIN_NAMESPACE_VERSION
50
51#if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
52
53  /* Mersenne twister implementation optimized for vector operations.
54   *
55   * Reference: http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/
56   */
57  template<typename _UIntType, size_t __m,
58	   size_t __pos1, size_t __sl1, size_t __sl2,
59	   size_t __sr1, size_t __sr2,
60	   uint32_t __msk1, uint32_t __msk2,
61	   uint32_t __msk3, uint32_t __msk4,
62	   uint32_t __parity1, uint32_t __parity2,
63	   uint32_t __parity3, uint32_t __parity4>
64    class simd_fast_mersenne_twister_engine
65    {
66      static_assert(std::is_unsigned<_UIntType>::value, "template argument "
67		    "substituting _UIntType not an unsigned integral type");
68      static_assert(__sr1 < 32, "first right shift too large");
69      static_assert(__sr2 < 16, "second right shift too large");
70      static_assert(__sl1 < 32, "first left shift too large");
71      static_assert(__sl2 < 16, "second left shift too large");
72
73    public:
74      typedef _UIntType result_type;
75
76    private:
77      static constexpr size_t m_w = sizeof(result_type) * 8;
78      static constexpr size_t _M_nstate = __m / 128 + 1;
79      static constexpr size_t _M_nstate32 = _M_nstate * 4;
80
81      static_assert(std::is_unsigned<_UIntType>::value, "template argument "
82		    "substituting _UIntType not an unsigned integral type");
83      static_assert(__pos1 < _M_nstate, "POS1 not smaller than state size");
84      static_assert(16 % sizeof(_UIntType) == 0,
85		    "UIntType size must divide 16");
86
87    public:
88      static constexpr size_t state_size = _M_nstate * (16
89							/ sizeof(result_type));
90      static constexpr result_type default_seed = 5489u;
91
92      // constructors and member function
93      explicit
94      simd_fast_mersenne_twister_engine(result_type __sd = default_seed)
95      { seed(__sd); }
96
97      template<typename _Sseq, typename = typename
98	std::enable_if<!std::is_same<_Sseq,
99				     simd_fast_mersenne_twister_engine>::value>
100	       ::type>
101	explicit
102	simd_fast_mersenne_twister_engine(_Sseq& __q)
103	{ seed(__q); }
104
105      void
106      seed(result_type __sd = default_seed);
107
108      template<typename _Sseq>
109	typename std::enable_if<std::is_class<_Sseq>::value>::type
110	seed(_Sseq& __q);
111
112      static constexpr result_type
113      min()
114      { return 0; };
115
116      static constexpr result_type
117      max()
118      { return std::numeric_limits<result_type>::max(); }
119
120      void
121      discard(unsigned long long __z);
122
123      result_type
124      operator()()
125      {
126	if (__builtin_expect(_M_pos >= state_size, 0))
127	  _M_gen_rand();
128
129	return _M_stateT[_M_pos++];
130      }
131
132      template<typename _UIntType_2, size_t __m_2,
133	       size_t __pos1_2, size_t __sl1_2, size_t __sl2_2,
134	       size_t __sr1_2, size_t __sr2_2,
135	       uint32_t __msk1_2, uint32_t __msk2_2,
136	       uint32_t __msk3_2, uint32_t __msk4_2,
137	       uint32_t __parity1_2, uint32_t __parity2_2,
138	       uint32_t __parity3_2, uint32_t __parity4_2>
139	friend bool
140	operator==(const simd_fast_mersenne_twister_engine<_UIntType_2,
141		   __m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2,
142		   __msk1_2, __msk2_2, __msk3_2, __msk4_2,
143		   __parity1_2, __parity2_2, __parity3_2, __parity4_2>& __lhs,
144		   const simd_fast_mersenne_twister_engine<_UIntType_2,
145		   __m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2,
146		   __msk1_2, __msk2_2, __msk3_2, __msk4_2,
147		   __parity1_2, __parity2_2, __parity3_2, __parity4_2>& __rhs);
148
149      template<typename _UIntType_2, size_t __m_2,
150	       size_t __pos1_2, size_t __sl1_2, size_t __sl2_2,
151	       size_t __sr1_2, size_t __sr2_2,
152	       uint32_t __msk1_2, uint32_t __msk2_2,
153	       uint32_t __msk3_2, uint32_t __msk4_2,
154	       uint32_t __parity1_2, uint32_t __parity2_2,
155	       uint32_t __parity3_2, uint32_t __parity4_2,
156	       typename _CharT, typename _Traits>
157	friend std::basic_ostream<_CharT, _Traits>&
158	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
159		   const __gnu_cxx::simd_fast_mersenne_twister_engine
160		   <_UIntType_2,
161		   __m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2,
162		   __msk1_2, __msk2_2, __msk3_2, __msk4_2,
163		   __parity1_2, __parity2_2, __parity3_2, __parity4_2>& __x);
164
165      template<typename _UIntType_2, size_t __m_2,
166	       size_t __pos1_2, size_t __sl1_2, size_t __sl2_2,
167	       size_t __sr1_2, size_t __sr2_2,
168	       uint32_t __msk1_2, uint32_t __msk2_2,
169	       uint32_t __msk3_2, uint32_t __msk4_2,
170	       uint32_t __parity1_2, uint32_t __parity2_2,
171	       uint32_t __parity3_2, uint32_t __parity4_2,
172	       typename _CharT, typename _Traits>
173	friend std::basic_istream<_CharT, _Traits>&
174	operator>>(std::basic_istream<_CharT, _Traits>& __is,
175		   __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType_2,
176		   __m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2,
177		   __msk1_2, __msk2_2, __msk3_2, __msk4_2,
178		   __parity1_2, __parity2_2, __parity3_2, __parity4_2>& __x);
179
180    private:
181      union
182      {
183#ifdef __SSE2__
184	__m128i _M_state[_M_nstate];
185#endif
186	uint32_t _M_state32[_M_nstate32];
187	result_type _M_stateT[state_size];
188      } __attribute__ ((__aligned__ (16)));
189      size_t _M_pos;
190
191      void _M_gen_rand(void);
192      void _M_period_certification();
193  };
194
195
196  template<typename _UIntType, size_t __m,
197	   size_t __pos1, size_t __sl1, size_t __sl2,
198	   size_t __sr1, size_t __sr2,
199	   uint32_t __msk1, uint32_t __msk2,
200	   uint32_t __msk3, uint32_t __msk4,
201	   uint32_t __parity1, uint32_t __parity2,
202	   uint32_t __parity3, uint32_t __parity4>
203    inline bool
204    operator!=(const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
205	       __m, __pos1, __sl1, __sl2, __sr1, __sr2, __msk1, __msk2, __msk3,
206	       __msk4, __parity1, __parity2, __parity3, __parity4>& __lhs,
207	       const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
208	       __m, __pos1, __sl1, __sl2, __sr1, __sr2, __msk1, __msk2, __msk3,
209	       __msk4, __parity1, __parity2, __parity3, __parity4>& __rhs)
210    { return !(__lhs == __rhs); }
211
212
213  /* Definitions for the SIMD-oriented Fast Mersenne Twister as defined
214   * in the C implementation by Daito and Matsumoto, as both a 32-bit
215   * and 64-bit version.
216   */
217  typedef simd_fast_mersenne_twister_engine<uint32_t, 607, 2,
218					    15, 3, 13, 3,
219					    0xfdff37ffU, 0xef7f3f7dU,
220					    0xff777b7dU, 0x7ff7fb2fU,
221					    0x00000001U, 0x00000000U,
222					    0x00000000U, 0x5986f054U>
223    sfmt607;
224
225  typedef simd_fast_mersenne_twister_engine<uint64_t, 607, 2,
226					    15, 3, 13, 3,
227					    0xfdff37ffU, 0xef7f3f7dU,
228					    0xff777b7dU, 0x7ff7fb2fU,
229					    0x00000001U, 0x00000000U,
230					    0x00000000U, 0x5986f054U>
231    sfmt607_64;
232
233
234  typedef simd_fast_mersenne_twister_engine<uint32_t, 1279, 7,
235					    14, 3, 5, 1,
236					    0xf7fefffdU, 0x7fefcfffU,
237					    0xaff3ef3fU, 0xb5ffff7fU,
238					    0x00000001U, 0x00000000U,
239					    0x00000000U, 0x20000000U>
240    sfmt1279;
241
242  typedef simd_fast_mersenne_twister_engine<uint64_t, 1279, 7,
243					    14, 3, 5, 1,
244					    0xf7fefffdU, 0x7fefcfffU,
245					    0xaff3ef3fU, 0xb5ffff7fU,
246					    0x00000001U, 0x00000000U,
247					    0x00000000U, 0x20000000U>
248    sfmt1279_64;
249
250
251  typedef simd_fast_mersenne_twister_engine<uint32_t, 2281, 12,
252					    19, 1, 5, 1,
253					    0xbff7ffbfU, 0xfdfffffeU,
254					    0xf7ffef7fU, 0xf2f7cbbfU,
255					    0x00000001U, 0x00000000U,
256					    0x00000000U, 0x41dfa600U>
257    sfmt2281;
258
259  typedef simd_fast_mersenne_twister_engine<uint64_t, 2281, 12,
260					    19, 1, 5, 1,
261					    0xbff7ffbfU, 0xfdfffffeU,
262					    0xf7ffef7fU, 0xf2f7cbbfU,
263					    0x00000001U, 0x00000000U,
264					    0x00000000U, 0x41dfa600U>
265    sfmt2281_64;
266
267
268  typedef simd_fast_mersenne_twister_engine<uint32_t, 4253, 17,
269					    20, 1, 7, 1,
270					    0x9f7bffffU, 0x9fffff5fU,
271					    0x3efffffbU, 0xfffff7bbU,
272					    0xa8000001U, 0xaf5390a3U,
273					    0xb740b3f8U, 0x6c11486dU>
274    sfmt4253;
275
276  typedef simd_fast_mersenne_twister_engine<uint64_t, 4253, 17,
277					    20, 1, 7, 1,
278					    0x9f7bffffU, 0x9fffff5fU,
279					    0x3efffffbU, 0xfffff7bbU,
280					    0xa8000001U, 0xaf5390a3U,
281					    0xb740b3f8U, 0x6c11486dU>
282    sfmt4253_64;
283
284
285  typedef simd_fast_mersenne_twister_engine<uint32_t, 11213, 68,
286					    14, 3, 7, 3,
287					    0xeffff7fbU, 0xffffffefU,
288					    0xdfdfbfffU, 0x7fffdbfdU,
289					    0x00000001U, 0x00000000U,
290					    0xe8148000U, 0xd0c7afa3U>
291    sfmt11213;
292
293  typedef simd_fast_mersenne_twister_engine<uint64_t, 11213, 68,
294					    14, 3, 7, 3,
295					    0xeffff7fbU, 0xffffffefU,
296					    0xdfdfbfffU, 0x7fffdbfdU,
297					    0x00000001U, 0x00000000U,
298					    0xe8148000U, 0xd0c7afa3U>
299    sfmt11213_64;
300
301
302  typedef simd_fast_mersenne_twister_engine<uint32_t, 19937, 122,
303					    18, 1, 11, 1,
304					    0xdfffffefU, 0xddfecb7fU,
305					    0xbffaffffU, 0xbffffff6U,
306					    0x00000001U, 0x00000000U,
307					    0x00000000U, 0x13c9e684U>
308    sfmt19937;
309
310  typedef simd_fast_mersenne_twister_engine<uint64_t, 19937, 122,
311					    18, 1, 11, 1,
312					    0xdfffffefU, 0xddfecb7fU,
313					    0xbffaffffU, 0xbffffff6U,
314					    0x00000001U, 0x00000000U,
315					    0x00000000U, 0x13c9e684U>
316    sfmt19937_64;
317
318
319  typedef simd_fast_mersenne_twister_engine<uint32_t, 44497, 330,
320					    5, 3, 9, 3,
321					    0xeffffffbU, 0xdfbebfffU,
322					    0xbfbf7befU, 0x9ffd7bffU,
323					    0x00000001U, 0x00000000U,
324					    0xa3ac4000U, 0xecc1327aU>
325    sfmt44497;
326
327  typedef simd_fast_mersenne_twister_engine<uint64_t, 44497, 330,
328					    5, 3, 9, 3,
329					    0xeffffffbU, 0xdfbebfffU,
330					    0xbfbf7befU, 0x9ffd7bffU,
331					    0x00000001U, 0x00000000U,
332					    0xa3ac4000U, 0xecc1327aU>
333    sfmt44497_64;
334
335
336  typedef simd_fast_mersenne_twister_engine<uint32_t, 86243, 366,
337					    6, 7, 19, 1,
338					    0xfdbffbffU, 0xbff7ff3fU,
339					    0xfd77efffU, 0xbf9ff3ffU,
340					    0x00000001U, 0x00000000U,
341					    0x00000000U, 0xe9528d85U>
342    sfmt86243;
343
344  typedef simd_fast_mersenne_twister_engine<uint64_t, 86243, 366,
345					    6, 7, 19, 1,
346					    0xfdbffbffU, 0xbff7ff3fU,
347					    0xfd77efffU, 0xbf9ff3ffU,
348					    0x00000001U, 0x00000000U,
349					    0x00000000U, 0xe9528d85U>
350    sfmt86243_64;
351
352
353  typedef simd_fast_mersenne_twister_engine<uint32_t, 132049, 110,
354					    19, 1, 21, 1,
355					    0xffffbb5fU, 0xfb6ebf95U,
356					    0xfffefffaU, 0xcff77fffU,
357					    0x00000001U, 0x00000000U,
358					    0xcb520000U, 0xc7e91c7dU>
359    sfmt132049;
360
361  typedef simd_fast_mersenne_twister_engine<uint64_t, 132049, 110,
362					    19, 1, 21, 1,
363					    0xffffbb5fU, 0xfb6ebf95U,
364					    0xfffefffaU, 0xcff77fffU,
365					    0x00000001U, 0x00000000U,
366					    0xcb520000U, 0xc7e91c7dU>
367    sfmt132049_64;
368
369
370  typedef simd_fast_mersenne_twister_engine<uint32_t, 216091, 627,
371					    11, 3, 10, 1,
372					    0xbff7bff7U, 0xbfffffffU,
373					    0xbffffa7fU, 0xffddfbfbU,
374					    0xf8000001U, 0x89e80709U,
375					    0x3bd2b64bU, 0x0c64b1e4U>
376    sfmt216091;
377
378  typedef simd_fast_mersenne_twister_engine<uint64_t, 216091, 627,
379					    11, 3, 10, 1,
380					    0xbff7bff7U, 0xbfffffffU,
381					    0xbffffa7fU, 0xffddfbfbU,
382					    0xf8000001U, 0x89e80709U,
383					    0x3bd2b64bU, 0x0c64b1e4U>
384    sfmt216091_64;
385
386#endif // __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
387
388  /**
389   * @brief A beta continuous distribution for random numbers.
390   *
391   * The formula for the beta probability density function is:
392   * @f[
393   *     p(x|\alpha,\beta) = \frac{1}{B(\alpha,\beta)}
394   *                         x^{\alpha - 1} (1 - x)^{\beta - 1}
395   * @f]
396   */
397  template<typename _RealType = double>
398    class beta_distribution
399    {
400      static_assert(std::is_floating_point<_RealType>::value,
401		    "template argument not a floating point type");
402
403    public:
404      /** The type of the range of the distribution. */
405      typedef _RealType result_type;
406      /** Parameter type. */
407      struct param_type
408      {
409	typedef beta_distribution<_RealType> distribution_type;
410	friend class beta_distribution<_RealType>;
411
412	explicit
413	param_type(_RealType __alpha_val = _RealType(1),
414		   _RealType __beta_val = _RealType(1))
415	: _M_alpha(__alpha_val), _M_beta(__beta_val)
416	{
417	  _GLIBCXX_DEBUG_ASSERT(_M_alpha > _RealType(0));
418	  _GLIBCXX_DEBUG_ASSERT(_M_beta > _RealType(0));
419	}
420
421	_RealType
422	alpha() const
423	{ return _M_alpha; }
424
425	_RealType
426	beta() const
427	{ return _M_beta; }
428
429	friend bool
430	operator==(const param_type& __p1, const param_type& __p2)
431	{ return (__p1._M_alpha == __p2._M_alpha
432		  && __p1._M_beta == __p2._M_beta); }
433
434      private:
435	void
436	_M_initialize();
437
438	_RealType _M_alpha;
439	_RealType _M_beta;
440      };
441
442    public:
443      /**
444       * @brief Constructs a beta distribution with parameters
445       * @f$\alpha@f$ and @f$\beta@f$.
446       */
447      explicit
448      beta_distribution(_RealType __alpha_val = _RealType(1),
449			_RealType __beta_val = _RealType(1))
450      : _M_param(__alpha_val, __beta_val)
451      { }
452
453      explicit
454      beta_distribution(const param_type& __p)
455      : _M_param(__p)
456      { }
457
458      /**
459       * @brief Resets the distribution state.
460       */
461      void
462      reset()
463      { }
464
465      /**
466       * @brief Returns the @f$\alpha@f$ of the distribution.
467       */
468      _RealType
469      alpha() const
470      { return _M_param.alpha(); }
471
472      /**
473       * @brief Returns the @f$\beta@f$ of the distribution.
474       */
475      _RealType
476      beta() const
477      { return _M_param.beta(); }
478
479      /**
480       * @brief Returns the parameter set of the distribution.
481       */
482      param_type
483      param() const
484      { return _M_param; }
485
486      /**
487       * @brief Sets the parameter set of the distribution.
488       * @param __param The new parameter set of the distribution.
489       */
490      void
491      param(const param_type& __param)
492      { _M_param = __param; }
493
494      /**
495       * @brief Returns the greatest lower bound value of the distribution.
496       */
497      result_type
498      min() const
499      { return result_type(0); }
500
501      /**
502       * @brief Returns the least upper bound value of the distribution.
503       */
504      result_type
505      max() const
506      { return result_type(1); }
507
508      /**
509       * @brief Generating functions.
510       */
511      template<typename _UniformRandomNumberGenerator>
512	result_type
513	operator()(_UniformRandomNumberGenerator& __urng)
514	{ return this->operator()(__urng, _M_param); }
515
516      template<typename _UniformRandomNumberGenerator>
517	result_type
518	operator()(_UniformRandomNumberGenerator& __urng,
519		   const param_type& __p);
520
521      template<typename _ForwardIterator,
522	       typename _UniformRandomNumberGenerator>
523	void
524	__generate(_ForwardIterator __f, _ForwardIterator __t,
525		   _UniformRandomNumberGenerator& __urng)
526	{ this->__generate(__f, __t, __urng, _M_param); }
527
528      template<typename _ForwardIterator,
529	       typename _UniformRandomNumberGenerator>
530	void
531	__generate(_ForwardIterator __f, _ForwardIterator __t,
532		   _UniformRandomNumberGenerator& __urng,
533		   const param_type& __p)
534	{ this->__generate_impl(__f, __t, __urng, __p); }
535
536      template<typename _UniformRandomNumberGenerator>
537	void
538	__generate(result_type* __f, result_type* __t,
539		   _UniformRandomNumberGenerator& __urng,
540		   const param_type& __p)
541	{ this->__generate_impl(__f, __t, __urng, __p); }
542
543      /**
544       * @brief Return true if two beta distributions have the same
545       *        parameters and the sequences that would be generated
546       *        are equal.
547       */
548      friend bool
549      operator==(const beta_distribution& __d1,
550		 const beta_distribution& __d2)
551      { return __d1._M_param == __d2._M_param; }
552
553      /**
554       * @brief Inserts a %beta_distribution random number distribution
555       * @p __x into the output stream @p __os.
556       *
557       * @param __os An output stream.
558       * @param __x  A %beta_distribution random number distribution.
559       *
560       * @returns The output stream with the state of @p __x inserted or in
561       * an error state.
562       */
563      template<typename _RealType1, typename _CharT, typename _Traits>
564	friend std::basic_ostream<_CharT, _Traits>&
565	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
566		   const __gnu_cxx::beta_distribution<_RealType1>& __x);
567
568      /**
569       * @brief Extracts a %beta_distribution random number distribution
570       * @p __x from the input stream @p __is.
571       *
572       * @param __is An input stream.
573       * @param __x  A %beta_distribution random number generator engine.
574       *
575       * @returns The input stream with @p __x extracted or in an error state.
576       */
577      template<typename _RealType1, typename _CharT, typename _Traits>
578	friend std::basic_istream<_CharT, _Traits>&
579	operator>>(std::basic_istream<_CharT, _Traits>& __is,
580		   __gnu_cxx::beta_distribution<_RealType1>& __x);
581
582    private:
583      template<typename _ForwardIterator,
584	       typename _UniformRandomNumberGenerator>
585	void
586	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
587			_UniformRandomNumberGenerator& __urng,
588			const param_type& __p);
589
590      param_type _M_param;
591    };
592
593  /**
594   * @brief Return true if two beta distributions are different.
595   */
596  template<typename _RealType>
597    inline bool
598    operator!=(const __gnu_cxx::beta_distribution<_RealType>& __d1,
599	       const __gnu_cxx::beta_distribution<_RealType>& __d2)
600   { return !(__d1 == __d2); }
601
602
603  /**
604   * @brief A multi-variate normal continuous distribution for random numbers.
605   *
606   * The formula for the normal probability density function is
607   * @f[
608   *     p(\overrightarrow{x}|\overrightarrow{\mu },\Sigma) =
609   *       \frac{1}{\sqrt{(2\pi )^k\det(\Sigma))}}
610   *       e^{-\frac{1}{2}(\overrightarrow{x}-\overrightarrow{\mu})^\text{T}
611   *          \Sigma ^{-1}(\overrightarrow{x}-\overrightarrow{\mu})}
612   * @f]
613   *
614   * where @f$\overrightarrow{x}@f$ and @f$\overrightarrow{\mu}@f$ are
615   * vectors of dimension @f$k@f$ and @f$\Sigma@f$ is the covariance
616   * matrix (which must be positive-definite).
617   */
618  template<std::size_t _Dimen, typename _RealType = double>
619    class normal_mv_distribution
620    {
621      static_assert(std::is_floating_point<_RealType>::value,
622		    "template argument not a floating point type");
623      static_assert(_Dimen != 0, "dimension is zero");
624
625    public:
626      /** The type of the range of the distribution. */
627      typedef std::array<_RealType, _Dimen> result_type;
628      /** Parameter type. */
629      class param_type
630      {
631	static constexpr size_t _M_t_size = _Dimen * (_Dimen + 1) / 2;
632
633      public:
634	typedef normal_mv_distribution<_Dimen, _RealType> distribution_type;
635	friend class normal_mv_distribution<_Dimen, _RealType>;
636
637	param_type()
638	{
639	  std::fill(_M_mean.begin(), _M_mean.end(), _RealType(0));
640	  auto __it = _M_t.begin();
641	  for (size_t __i = 0; __i < _Dimen; ++__i)
642	    {
643	      std::fill_n(__it, __i, _RealType(0));
644	      __it += __i;
645	      *__it++ = _RealType(1);
646	    }
647	}
648
649	template<typename _ForwardIterator1, typename _ForwardIterator2>
650	  param_type(_ForwardIterator1 __meanbegin,
651		     _ForwardIterator1 __meanend,
652		     _ForwardIterator2 __varcovbegin,
653		     _ForwardIterator2 __varcovend)
654	{
655	  __glibcxx_function_requires(_ForwardIteratorConcept<
656				      _ForwardIterator1>)
657	  __glibcxx_function_requires(_ForwardIteratorConcept<
658				      _ForwardIterator2>)
659	  _GLIBCXX_DEBUG_ASSERT(std::distance(__meanbegin, __meanend)
660				<= _Dimen);
661	  const auto __dist = std::distance(__varcovbegin, __varcovend);
662	  _GLIBCXX_DEBUG_ASSERT(__dist == _Dimen * _Dimen
663				|| __dist == _Dimen * (_Dimen + 1) / 2
664				|| __dist == _Dimen);
665
666	  if (__dist == _Dimen * _Dimen)
667	    _M_init_full(__meanbegin, __meanend, __varcovbegin, __varcovend);
668	  else if (__dist == _Dimen * (_Dimen + 1) / 2)
669	    _M_init_lower(__meanbegin, __meanend, __varcovbegin, __varcovend);
670	  else
671	    _M_init_diagonal(__meanbegin, __meanend,
672			     __varcovbegin, __varcovend);
673	}
674
675	param_type(std::initializer_list<_RealType> __mean,
676		   std::initializer_list<_RealType> __varcov)
677	{
678	  _GLIBCXX_DEBUG_ASSERT(__mean.size() <= _Dimen);
679	  _GLIBCXX_DEBUG_ASSERT(__varcov.size() == _Dimen * _Dimen
680				|| __varcov.size() == _Dimen * (_Dimen + 1) / 2
681				|| __varcov.size() == _Dimen);
682
683	  if (__varcov.size() == _Dimen * _Dimen)
684	    _M_init_full(__mean.begin(), __mean.end(),
685			 __varcov.begin(), __varcov.end());
686	  else if (__varcov.size() == _Dimen * (_Dimen + 1) / 2)
687	    _M_init_lower(__mean.begin(), __mean.end(),
688			  __varcov.begin(), __varcov.end());
689	  else
690	    _M_init_diagonal(__mean.begin(), __mean.end(),
691			     __varcov.begin(), __varcov.end());
692	}
693
694	std::array<_RealType, _Dimen>
695	mean() const
696	{ return _M_mean; }
697
698	std::array<_RealType, _M_t_size>
699	varcov() const
700	{ return _M_t; }
701
702	friend bool
703	operator==(const param_type& __p1, const param_type& __p2)
704	{ return __p1._M_mean == __p2._M_mean && __p1._M_t == __p2._M_t; }
705
706      private:
707	template <typename _InputIterator1, typename _InputIterator2>
708	  void _M_init_full(_InputIterator1 __meanbegin,
709			    _InputIterator1 __meanend,
710			    _InputIterator2 __varcovbegin,
711			    _InputIterator2 __varcovend);
712	template <typename _InputIterator1, typename _InputIterator2>
713	  void _M_init_lower(_InputIterator1 __meanbegin,
714			     _InputIterator1 __meanend,
715			     _InputIterator2 __varcovbegin,
716			     _InputIterator2 __varcovend);
717	template <typename _InputIterator1, typename _InputIterator2>
718	  void _M_init_diagonal(_InputIterator1 __meanbegin,
719				_InputIterator1 __meanend,
720				_InputIterator2 __varbegin,
721				_InputIterator2 __varend);
722
723	std::array<_RealType, _Dimen> _M_mean;
724	std::array<_RealType, _M_t_size> _M_t;
725      };
726
727    public:
728      normal_mv_distribution()
729      : _M_param(), _M_nd()
730      { }
731
732      template<typename _ForwardIterator1, typename _ForwardIterator2>
733	normal_mv_distribution(_ForwardIterator1 __meanbegin,
734			       _ForwardIterator1 __meanend,
735			       _ForwardIterator2 __varcovbegin,
736			       _ForwardIterator2 __varcovend)
737	: _M_param(__meanbegin, __meanend, __varcovbegin, __varcovend),
738	  _M_nd()
739	{ }
740
741      normal_mv_distribution(std::initializer_list<_RealType> __mean,
742			     std::initializer_list<_RealType> __varcov)
743      : _M_param(__mean, __varcov), _M_nd()
744      { }
745
746      explicit
747      normal_mv_distribution(const param_type& __p)
748      : _M_param(__p), _M_nd()
749      { }
750
751      /**
752       * @brief Resets the distribution state.
753       */
754      void
755      reset()
756      { _M_nd.reset(); }
757
758      /**
759       * @brief Returns the mean of the distribution.
760       */
761      result_type
762      mean() const
763      { return _M_param.mean(); }
764
765      /**
766       * @brief Returns the compact form of the variance/covariance
767       * matrix of the distribution.
768       */
769      std::array<_RealType, _Dimen * (_Dimen + 1) / 2>
770      varcov() const
771      { return _M_param.varcov(); }
772
773      /**
774       * @brief Returns the parameter set of the distribution.
775       */
776      param_type
777      param() const
778      { return _M_param; }
779
780      /**
781       * @brief Sets the parameter set of the distribution.
782       * @param __param The new parameter set of the distribution.
783       */
784      void
785      param(const param_type& __param)
786      { _M_param = __param; }
787
788      /**
789       * @brief Returns the greatest lower bound value of the distribution.
790       */
791      result_type
792      min() const
793      { result_type __res;
794	__res.fill(std::numeric_limits<_RealType>::lowest());
795	return __res; }
796
797      /**
798       * @brief Returns the least upper bound value of the distribution.
799       */
800      result_type
801      max() const
802      { result_type __res;
803	__res.fill(std::numeric_limits<_RealType>::max());
804	return __res; }
805
806      /**
807       * @brief Generating functions.
808       */
809      template<typename _UniformRandomNumberGenerator>
810	result_type
811	operator()(_UniformRandomNumberGenerator& __urng)
812	{ return this->operator()(__urng, _M_param); }
813
814      template<typename _UniformRandomNumberGenerator>
815	result_type
816	operator()(_UniformRandomNumberGenerator& __urng,
817		   const param_type& __p);
818
819      template<typename _ForwardIterator,
820	       typename _UniformRandomNumberGenerator>
821	void
822	__generate(_ForwardIterator __f, _ForwardIterator __t,
823		   _UniformRandomNumberGenerator& __urng)
824	{ return this->__generate_impl(__f, __t, __urng, _M_param); }
825
826      template<typename _ForwardIterator,
827	       typename _UniformRandomNumberGenerator>
828	void
829	__generate(_ForwardIterator __f, _ForwardIterator __t,
830		   _UniformRandomNumberGenerator& __urng,
831		   const param_type& __p)
832	{ return this->__generate_impl(__f, __t, __urng, __p); }
833
834      /**
835       * @brief Return true if two multi-variant normal distributions have
836       *        the same parameters and the sequences that would
837       *        be generated are equal.
838       */
839      template<size_t _Dimen1, typename _RealType1>
840	friend bool
841	operator==(const
842		   __gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>&
843		   __d1,
844		   const
845		   __gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>&
846		   __d2);
847
848      /**
849       * @brief Inserts a %normal_mv_distribution random number distribution
850       * @p __x into the output stream @p __os.
851       *
852       * @param __os An output stream.
853       * @param __x  A %normal_mv_distribution random number distribution.
854       *
855       * @returns The output stream with the state of @p __x inserted or in
856       * an error state.
857       */
858      template<size_t _Dimen1, typename _RealType1,
859	       typename _CharT, typename _Traits>
860	friend std::basic_ostream<_CharT, _Traits>&
861	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
862		   const
863		   __gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>&
864		   __x);
865
866      /**
867       * @brief Extracts a %normal_mv_distribution random number distribution
868       * @p __x from the input stream @p __is.
869       *
870       * @param __is An input stream.
871       * @param __x  A %normal_mv_distribution random number generator engine.
872       *
873       * @returns The input stream with @p __x extracted or in an error
874       *          state.
875       */
876      template<size_t _Dimen1, typename _RealType1,
877	       typename _CharT, typename _Traits>
878	friend std::basic_istream<_CharT, _Traits>&
879	operator>>(std::basic_istream<_CharT, _Traits>& __is,
880		   __gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>&
881		   __x);
882
883    private:
884      template<typename _ForwardIterator,
885	       typename _UniformRandomNumberGenerator>
886	void
887	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
888			_UniformRandomNumberGenerator& __urng,
889			const param_type& __p);
890
891      param_type _M_param;
892      std::normal_distribution<_RealType> _M_nd;
893  };
894
895  /**
896   * @brief Return true if two multi-variate normal distributions are
897   * different.
898   */
899  template<size_t _Dimen, typename _RealType>
900    inline bool
901    operator!=(const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
902	       __d1,
903	       const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
904	       __d2)
905    { return !(__d1 == __d2); }
906
907
908  /**
909   * @brief A Rice continuous distribution for random numbers.
910   *
911   * The formula for the Rice probability density function is
912   * @f[
913   *     p(x|\nu,\sigma) = \frac{x}{\sigma^2}
914   *                       \exp\left(-\frac{x^2+\nu^2}{2\sigma^2}\right)
915   *                       I_0\left(\frac{x \nu}{\sigma^2}\right)
916   * @f]
917   * where @f$I_0(z)@f$ is the modified Bessel function of the first kind
918   * of order 0 and @f$\nu >= 0@f$ and @f$\sigma > 0@f$.
919   *
920   * <table border=1 cellpadding=10 cellspacing=0>
921   * <caption align=top>Distribution Statistics</caption>
922   * <tr><td>Mean</td><td>@f$\sqrt{\pi/2}L_{1/2}(-\nu^2/2\sigma^2)@f$</td></tr>
923   * <tr><td>Variance</td><td>@f$2\sigma^2 + \nu^2
924   *                   + (\pi\sigma^2/2)L^2_{1/2}(-\nu^2/2\sigma^2)@f$</td></tr>
925   * <tr><td>Range</td><td>@f$[0, \infty)@f$</td></tr>
926   * </table>
927   * where @f$L_{1/2}(x)@f$ is the Laguerre polynomial of order 1/2.
928   */
929  template<typename _RealType = double>
930    class
931    rice_distribution
932    {
933      static_assert(std::is_floating_point<_RealType>::value,
934		    "template argument not a floating point type");
935    public:
936      /** The type of the range of the distribution. */
937      typedef _RealType result_type;
938      /** Parameter type. */
939      struct param_type
940      {
941	typedef rice_distribution<result_type> distribution_type;
942
943	param_type(result_type __nu_val = result_type(0),
944		   result_type __sigma_val = result_type(1))
945	: _M_nu(__nu_val), _M_sigma(__sigma_val)
946	{
947	  _GLIBCXX_DEBUG_ASSERT(_M_nu >= result_type(0));
948	  _GLIBCXX_DEBUG_ASSERT(_M_sigma > result_type(0));
949	}
950
951	result_type
952	nu() const
953	{ return _M_nu; }
954
955	result_type
956	sigma() const
957	{ return _M_sigma; }
958
959	friend bool
960	operator==(const param_type& __p1, const param_type& __p2)
961	{ return __p1._M_nu == __p2._M_nu
962	      && __p1._M_sigma == __p2._M_sigma; }
963
964      private:
965	void _M_initialize();
966
967	result_type _M_nu;
968	result_type _M_sigma;
969      };
970
971      /**
972       * @brief Constructors.
973       */
974      explicit
975      rice_distribution(result_type __nu_val = result_type(0),
976			result_type __sigma_val = result_type(1))
977      : _M_param(__nu_val, __sigma_val),
978	_M_ndx(__nu_val, __sigma_val),
979	_M_ndy(result_type(0), __sigma_val)
980      { }
981
982      explicit
983      rice_distribution(const param_type& __p)
984      : _M_param(__p),
985	_M_ndx(__p.nu(), __p.sigma()),
986	_M_ndy(result_type(0), __p.sigma())
987      { }
988
989      /**
990       * @brief Resets the distribution state.
991       */
992      void
993      reset()
994      {
995	_M_ndx.reset();
996	_M_ndy.reset();
997      }
998
999      /**
1000       * @brief Return the parameters of the distribution.
1001       */
1002      result_type
1003      nu() const
1004      { return _M_param.nu(); }
1005
1006      result_type
1007      sigma() const
1008      { return _M_param.sigma(); }
1009
1010      /**
1011       * @brief Returns the parameter set of the distribution.
1012       */
1013      param_type
1014      param() const
1015      { return _M_param; }
1016
1017      /**
1018       * @brief Sets the parameter set of the distribution.
1019       * @param __param The new parameter set of the distribution.
1020       */
1021      void
1022      param(const param_type& __param)
1023      { _M_param = __param; }
1024
1025      /**
1026       * @brief Returns the greatest lower bound value of the distribution.
1027       */
1028      result_type
1029      min() const
1030      { return result_type(0); }
1031
1032      /**
1033       * @brief Returns the least upper bound value of the distribution.
1034       */
1035      result_type
1036      max() const
1037      { return std::numeric_limits<result_type>::max(); }
1038
1039      /**
1040       * @brief Generating functions.
1041       */
1042      template<typename _UniformRandomNumberGenerator>
1043	result_type
1044	operator()(_UniformRandomNumberGenerator& __urng)
1045	{
1046	  result_type __x = this->_M_ndx(__urng);
1047	  result_type __y = this->_M_ndy(__urng);
1048#if _GLIBCXX_USE_C99_MATH_TR1
1049	  return std::hypot(__x, __y);
1050#else
1051	  return std::sqrt(__x * __x + __y * __y);
1052#endif
1053	}
1054
1055      template<typename _UniformRandomNumberGenerator>
1056	result_type
1057	operator()(_UniformRandomNumberGenerator& __urng,
1058		   const param_type& __p)
1059	{
1060	  typename std::normal_distribution<result_type>::param_type
1061	    __px(__p.nu(), __p.sigma()), __py(result_type(0), __p.sigma());
1062	  result_type __x = this->_M_ndx(__px, __urng);
1063	  result_type __y = this->_M_ndy(__py, __urng);
1064#if _GLIBCXX_USE_C99_MATH_TR1
1065	  return std::hypot(__x, __y);
1066#else
1067	  return std::sqrt(__x * __x + __y * __y);
1068#endif
1069	}
1070
1071      template<typename _ForwardIterator,
1072	       typename _UniformRandomNumberGenerator>
1073	void
1074	__generate(_ForwardIterator __f, _ForwardIterator __t,
1075		   _UniformRandomNumberGenerator& __urng)
1076	{ this->__generate(__f, __t, __urng, _M_param); }
1077
1078      template<typename _ForwardIterator,
1079	       typename _UniformRandomNumberGenerator>
1080	void
1081	__generate(_ForwardIterator __f, _ForwardIterator __t,
1082		   _UniformRandomNumberGenerator& __urng,
1083		   const param_type& __p)
1084	{ this->__generate_impl(__f, __t, __urng, __p); }
1085
1086      template<typename _UniformRandomNumberGenerator>
1087	void
1088	__generate(result_type* __f, result_type* __t,
1089		   _UniformRandomNumberGenerator& __urng,
1090		   const param_type& __p)
1091	{ this->__generate_impl(__f, __t, __urng, __p); }
1092
1093      /**
1094       * @brief Return true if two Rice distributions have
1095       *        the same parameters and the sequences that would
1096       *        be generated are equal.
1097       */
1098      friend bool
1099      operator==(const rice_distribution& __d1,
1100		 const rice_distribution& __d2)
1101      { return (__d1._M_param == __d2._M_param
1102		&& __d1._M_ndx == __d2._M_ndx
1103		&& __d1._M_ndy == __d2._M_ndy); }
1104
1105      /**
1106       * @brief Inserts a %rice_distribution random number distribution
1107       * @p __x into the output stream @p __os.
1108       *
1109       * @param __os An output stream.
1110       * @param __x  A %rice_distribution random number distribution.
1111       *
1112       * @returns The output stream with the state of @p __x inserted or in
1113       * an error state.
1114       */
1115      template<typename _RealType1, typename _CharT, typename _Traits>
1116	friend std::basic_ostream<_CharT, _Traits>&
1117	operator<<(std::basic_ostream<_CharT, _Traits>&,
1118		   const rice_distribution<_RealType1>&);
1119
1120      /**
1121       * @brief Extracts a %rice_distribution random number distribution
1122       * @p __x from the input stream @p __is.
1123       *
1124       * @param __is An input stream.
1125       * @param __x A %rice_distribution random number
1126       *            generator engine.
1127       *
1128       * @returns The input stream with @p __x extracted or in an error state.
1129       */
1130      template<typename _RealType1, typename _CharT, typename _Traits>
1131	friend std::basic_istream<_CharT, _Traits>&
1132	operator>>(std::basic_istream<_CharT, _Traits>&,
1133		   rice_distribution<_RealType1>&);
1134
1135    private:
1136      template<typename _ForwardIterator,
1137	       typename _UniformRandomNumberGenerator>
1138	void
1139	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1140			_UniformRandomNumberGenerator& __urng,
1141			const param_type& __p);
1142
1143      param_type _M_param;
1144
1145      std::normal_distribution<result_type> _M_ndx;
1146      std::normal_distribution<result_type> _M_ndy;
1147    };
1148
1149  /**
1150   * @brief Return true if two Rice distributions are not equal.
1151   */
1152  template<typename _RealType1>
1153    inline bool
1154    operator!=(const rice_distribution<_RealType1>& __d1,
1155	       const rice_distribution<_RealType1>& __d2)
1156    { return !(__d1 == __d2); }
1157
1158
1159  /**
1160   * @brief A Nakagami continuous distribution for random numbers.
1161   *
1162   * The formula for the Nakagami probability density function is
1163   * @f[
1164   *     p(x|\mu,\omega) = \frac{2\mu^\mu}{\Gamma(\mu)\omega^\mu}
1165   *                       x^{2\mu-1}e^{-\mu x / \omega}
1166   * @f]
1167   * where @f$\Gamma(z)@f$ is the gamma function and @f$\mu >= 0.5@f$
1168   * and @f$\omega > 0@f$.
1169   */
1170  template<typename _RealType = double>
1171    class
1172    nakagami_distribution
1173    {
1174      static_assert(std::is_floating_point<_RealType>::value,
1175		    "template argument not a floating point type");
1176
1177    public:
1178      /** The type of the range of the distribution. */
1179      typedef _RealType result_type;
1180      /** Parameter type. */
1181      struct param_type
1182      {
1183	typedef nakagami_distribution<result_type> distribution_type;
1184
1185	param_type(result_type __mu_val = result_type(1),
1186		   result_type __omega_val = result_type(1))
1187	: _M_mu(__mu_val), _M_omega(__omega_val)
1188	{
1189	  _GLIBCXX_DEBUG_ASSERT(_M_mu >= result_type(0.5L));
1190	  _GLIBCXX_DEBUG_ASSERT(_M_omega > result_type(0));
1191	}
1192
1193	result_type
1194	mu() const
1195	{ return _M_mu; }
1196
1197	result_type
1198	omega() const
1199	{ return _M_omega; }
1200
1201	friend bool
1202	operator==(const param_type& __p1, const param_type& __p2)
1203	{ return __p1._M_mu == __p2._M_mu
1204	      && __p1._M_omega == __p2._M_omega; }
1205
1206      private:
1207	void _M_initialize();
1208
1209	result_type _M_mu;
1210	result_type _M_omega;
1211      };
1212
1213      /**
1214       * @brief Constructors.
1215       */
1216      explicit
1217      nakagami_distribution(result_type __mu_val = result_type(1),
1218			    result_type __omega_val = result_type(1))
1219      : _M_param(__mu_val, __omega_val),
1220	_M_gd(__mu_val, __omega_val / __mu_val)
1221      { }
1222
1223      explicit
1224      nakagami_distribution(const param_type& __p)
1225      : _M_param(__p),
1226	_M_gd(__p.mu(), __p.omega() / __p.mu())
1227      { }
1228
1229      /**
1230       * @brief Resets the distribution state.
1231       */
1232      void
1233      reset()
1234      { _M_gd.reset(); }
1235
1236      /**
1237       * @brief Return the parameters of the distribution.
1238       */
1239      result_type
1240      mu() const
1241      { return _M_param.mu(); }
1242
1243      result_type
1244      omega() const
1245      { return _M_param.omega(); }
1246
1247      /**
1248       * @brief Returns the parameter set of the distribution.
1249       */
1250      param_type
1251      param() const
1252      { return _M_param; }
1253
1254      /**
1255       * @brief Sets the parameter set of the distribution.
1256       * @param __param The new parameter set of the distribution.
1257       */
1258      void
1259      param(const param_type& __param)
1260      { _M_param = __param; }
1261
1262      /**
1263       * @brief Returns the greatest lower bound value of the distribution.
1264       */
1265      result_type
1266      min() const
1267      { return result_type(0); }
1268
1269      /**
1270       * @brief Returns the least upper bound value of the distribution.
1271       */
1272      result_type
1273      max() const
1274      { return std::numeric_limits<result_type>::max(); }
1275
1276      /**
1277       * @brief Generating functions.
1278       */
1279      template<typename _UniformRandomNumberGenerator>
1280	result_type
1281	operator()(_UniformRandomNumberGenerator& __urng)
1282	{ return std::sqrt(this->_M_gd(__urng)); }
1283
1284      template<typename _UniformRandomNumberGenerator>
1285	result_type
1286	operator()(_UniformRandomNumberGenerator& __urng,
1287		   const param_type& __p)
1288	{
1289	  typename std::gamma_distribution<result_type>::param_type
1290	    __pg(__p.mu(), __p.omega() / __p.mu());
1291	  return std::sqrt(this->_M_gd(__pg, __urng));
1292	}
1293
1294      template<typename _ForwardIterator,
1295	       typename _UniformRandomNumberGenerator>
1296	void
1297	__generate(_ForwardIterator __f, _ForwardIterator __t,
1298		   _UniformRandomNumberGenerator& __urng)
1299	{ this->__generate(__f, __t, __urng, _M_param); }
1300
1301      template<typename _ForwardIterator,
1302	       typename _UniformRandomNumberGenerator>
1303	void
1304	__generate(_ForwardIterator __f, _ForwardIterator __t,
1305		   _UniformRandomNumberGenerator& __urng,
1306		   const param_type& __p)
1307	{ this->__generate_impl(__f, __t, __urng, __p); }
1308
1309      template<typename _UniformRandomNumberGenerator>
1310	void
1311	__generate(result_type* __f, result_type* __t,
1312		   _UniformRandomNumberGenerator& __urng,
1313		   const param_type& __p)
1314	{ this->__generate_impl(__f, __t, __urng, __p); }
1315
1316      /**
1317       * @brief Return true if two Nakagami distributions have
1318       *        the same parameters and the sequences that would
1319       *        be generated are equal.
1320       */
1321      friend bool
1322      operator==(const nakagami_distribution& __d1,
1323		 const nakagami_distribution& __d2)
1324      { return (__d1._M_param == __d2._M_param
1325		&& __d1._M_gd == __d2._M_gd); }
1326
1327      /**
1328       * @brief Inserts a %nakagami_distribution random number distribution
1329       * @p __x into the output stream @p __os.
1330       *
1331       * @param __os An output stream.
1332       * @param __x  A %nakagami_distribution random number distribution.
1333       *
1334       * @returns The output stream with the state of @p __x inserted or in
1335       * an error state.
1336       */
1337      template<typename _RealType1, typename _CharT, typename _Traits>
1338	friend std::basic_ostream<_CharT, _Traits>&
1339	operator<<(std::basic_ostream<_CharT, _Traits>&,
1340		   const nakagami_distribution<_RealType1>&);
1341
1342      /**
1343       * @brief Extracts a %nakagami_distribution random number distribution
1344       * @p __x from the input stream @p __is.
1345       *
1346       * @param __is An input stream.
1347       * @param __x A %nakagami_distribution random number
1348       *            generator engine.
1349       *
1350       * @returns The input stream with @p __x extracted or in an error state.
1351       */
1352      template<typename _RealType1, typename _CharT, typename _Traits>
1353	friend std::basic_istream<_CharT, _Traits>&
1354	operator>>(std::basic_istream<_CharT, _Traits>&,
1355		   nakagami_distribution<_RealType1>&);
1356
1357    private:
1358      template<typename _ForwardIterator,
1359	       typename _UniformRandomNumberGenerator>
1360	void
1361	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1362			_UniformRandomNumberGenerator& __urng,
1363			const param_type& __p);
1364
1365      param_type _M_param;
1366
1367      std::gamma_distribution<result_type> _M_gd;
1368    };
1369
1370  /**
1371   * @brief Return true if two Nakagami distributions are not equal.
1372   */
1373  template<typename _RealType>
1374    inline bool
1375    operator!=(const nakagami_distribution<_RealType>& __d1,
1376	       const nakagami_distribution<_RealType>& __d2)
1377    { return !(__d1 == __d2); }
1378
1379
1380  /**
1381   * @brief A Pareto continuous distribution for random numbers.
1382   *
1383   * The formula for the Pareto cumulative probability function is
1384   * @f[
1385   *     P(x|\alpha,\mu) = 1 - \left(\frac{\mu}{x}\right)^\alpha
1386   * @f]
1387   * The formula for the Pareto probability density function is
1388   * @f[
1389   *     p(x|\alpha,\mu) = \frac{\alpha + 1}{\mu}
1390   *                       \left(\frac{\mu}{x}\right)^{\alpha + 1}
1391   * @f]
1392   * where @f$x >= \mu@f$ and @f$\mu > 0@f$, @f$\alpha > 0@f$.
1393   *
1394   * <table border=1 cellpadding=10 cellspacing=0>
1395   * <caption align=top>Distribution Statistics</caption>
1396   * <tr><td>Mean</td><td>@f$\alpha \mu / (\alpha - 1)@f$
1397   *              for @f$\alpha > 1@f$</td></tr>
1398   * <tr><td>Variance</td><td>@f$\alpha \mu^2 / [(\alpha - 1)^2(\alpha - 2)]@f$
1399   *              for @f$\alpha > 2@f$</td></tr>
1400   * <tr><td>Range</td><td>@f$[\mu, \infty)@f$</td></tr>
1401   * </table>
1402   */
1403  template<typename _RealType = double>
1404    class
1405    pareto_distribution
1406    {
1407      static_assert(std::is_floating_point<_RealType>::value,
1408		    "template argument not a floating point type");
1409
1410    public:
1411      /** The type of the range of the distribution. */
1412      typedef _RealType result_type;
1413      /** Parameter type. */
1414      struct param_type
1415      {
1416	typedef pareto_distribution<result_type> distribution_type;
1417
1418	param_type(result_type __alpha_val = result_type(1),
1419		   result_type __mu_val = result_type(1))
1420	: _M_alpha(__alpha_val), _M_mu(__mu_val)
1421	{
1422	  _GLIBCXX_DEBUG_ASSERT(_M_alpha > result_type(0));
1423	  _GLIBCXX_DEBUG_ASSERT(_M_mu > result_type(0));
1424	}
1425
1426	result_type
1427	alpha() const
1428	{ return _M_alpha; }
1429
1430	result_type
1431	mu() const
1432	{ return _M_mu; }
1433
1434	friend bool
1435	operator==(const param_type& __p1, const param_type& __p2)
1436	{ return __p1._M_alpha == __p2._M_alpha && __p1._M_mu == __p2._M_mu; }
1437
1438      private:
1439	void _M_initialize();
1440
1441	result_type _M_alpha;
1442	result_type _M_mu;
1443      };
1444
1445      /**
1446       * @brief Constructors.
1447       */
1448      explicit
1449      pareto_distribution(result_type __alpha_val = result_type(1),
1450			  result_type __mu_val = result_type(1))
1451      : _M_param(__alpha_val, __mu_val),
1452	_M_ud()
1453      { }
1454
1455      explicit
1456      pareto_distribution(const param_type& __p)
1457      : _M_param(__p),
1458	_M_ud()
1459      { }
1460
1461      /**
1462       * @brief Resets the distribution state.
1463       */
1464      void
1465      reset()
1466      {
1467	_M_ud.reset();
1468      }
1469
1470      /**
1471       * @brief Return the parameters of the distribution.
1472       */
1473      result_type
1474      alpha() const
1475      { return _M_param.alpha(); }
1476
1477      result_type
1478      mu() const
1479      { return _M_param.mu(); }
1480
1481      /**
1482       * @brief Returns the parameter set of the distribution.
1483       */
1484      param_type
1485      param() const
1486      { return _M_param; }
1487
1488      /**
1489       * @brief Sets the parameter set of the distribution.
1490       * @param __param The new parameter set of the distribution.
1491       */
1492      void
1493      param(const param_type& __param)
1494      { _M_param = __param; }
1495
1496      /**
1497       * @brief Returns the greatest lower bound value of the distribution.
1498       */
1499      result_type
1500      min() const
1501      { return this->mu(); }
1502
1503      /**
1504       * @brief Returns the least upper bound value of the distribution.
1505       */
1506      result_type
1507      max() const
1508      { return std::numeric_limits<result_type>::max(); }
1509
1510      /**
1511       * @brief Generating functions.
1512       */
1513      template<typename _UniformRandomNumberGenerator>
1514	result_type
1515	operator()(_UniformRandomNumberGenerator& __urng)
1516	{
1517	  return this->mu() * std::pow(this->_M_ud(__urng),
1518				       -result_type(1) / this->alpha());
1519	}
1520
1521      template<typename _UniformRandomNumberGenerator>
1522	result_type
1523	operator()(_UniformRandomNumberGenerator& __urng,
1524		   const param_type& __p)
1525	{
1526	  return __p.mu() * std::pow(this->_M_ud(__urng),
1527					   -result_type(1) / __p.alpha());
1528	}
1529
1530      template<typename _ForwardIterator,
1531	       typename _UniformRandomNumberGenerator>
1532	void
1533	__generate(_ForwardIterator __f, _ForwardIterator __t,
1534		   _UniformRandomNumberGenerator& __urng)
1535	{ this->__generate(__f, __t, __urng, _M_param); }
1536
1537      template<typename _ForwardIterator,
1538	       typename _UniformRandomNumberGenerator>
1539	void
1540	__generate(_ForwardIterator __f, _ForwardIterator __t,
1541		   _UniformRandomNumberGenerator& __urng,
1542		   const param_type& __p)
1543	{ this->__generate_impl(__f, __t, __urng, __p); }
1544
1545      template<typename _UniformRandomNumberGenerator>
1546	void
1547	__generate(result_type* __f, result_type* __t,
1548		   _UniformRandomNumberGenerator& __urng,
1549		   const param_type& __p)
1550	{ this->__generate_impl(__f, __t, __urng, __p); }
1551
1552      /**
1553       * @brief Return true if two Pareto distributions have
1554       *        the same parameters and the sequences that would
1555       *        be generated are equal.
1556       */
1557      friend bool
1558      operator==(const pareto_distribution& __d1,
1559		 const pareto_distribution& __d2)
1560      { return (__d1._M_param == __d2._M_param
1561		&& __d1._M_ud == __d2._M_ud); }
1562
1563      /**
1564       * @brief Inserts a %pareto_distribution random number distribution
1565       * @p __x into the output stream @p __os.
1566       *
1567       * @param __os An output stream.
1568       * @param __x  A %pareto_distribution random number distribution.
1569       *
1570       * @returns The output stream with the state of @p __x inserted or in
1571       * an error state.
1572       */
1573      template<typename _RealType1, typename _CharT, typename _Traits>
1574	friend std::basic_ostream<_CharT, _Traits>&
1575	operator<<(std::basic_ostream<_CharT, _Traits>&,
1576		   const pareto_distribution<_RealType1>&);
1577
1578      /**
1579       * @brief Extracts a %pareto_distribution random number distribution
1580       * @p __x from the input stream @p __is.
1581       *
1582       * @param __is An input stream.
1583       * @param __x A %pareto_distribution random number
1584       *            generator engine.
1585       *
1586       * @returns The input stream with @p __x extracted or in an error state.
1587       */
1588      template<typename _RealType1, typename _CharT, typename _Traits>
1589	friend std::basic_istream<_CharT, _Traits>&
1590	operator>>(std::basic_istream<_CharT, _Traits>&,
1591		   pareto_distribution<_RealType1>&);
1592
1593    private:
1594      template<typename _ForwardIterator,
1595	       typename _UniformRandomNumberGenerator>
1596	void
1597	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1598			_UniformRandomNumberGenerator& __urng,
1599			const param_type& __p);
1600
1601      param_type _M_param;
1602
1603      std::uniform_real_distribution<result_type> _M_ud;
1604    };
1605
1606  /**
1607   * @brief Return true if two Pareto distributions are not equal.
1608   */
1609  template<typename _RealType>
1610    inline bool
1611    operator!=(const pareto_distribution<_RealType>& __d1,
1612	       const pareto_distribution<_RealType>& __d2)
1613    { return !(__d1 == __d2); }
1614
1615
1616  /**
1617   * @brief A K continuous distribution for random numbers.
1618   *
1619   * The formula for the K probability density function is
1620   * @f[
1621   *     p(x|\lambda, \mu, \nu) = \frac{2}{x}
1622   *             \left(\frac{\lambda\nu x}{\mu}\right)^{\frac{\lambda + \nu}{2}}
1623   *             \frac{1}{\Gamma(\lambda)\Gamma(\nu)}
1624   *             K_{\nu - \lambda}\left(2\sqrt{\frac{\lambda\nu x}{\mu}}\right)
1625   * @f]
1626   * where @f$I_0(z)@f$ is the modified Bessel function of the second kind
1627   * of order @f$\nu - \lambda@f$ and @f$\lambda > 0@f$, @f$\mu > 0@f$
1628   * and @f$\nu > 0@f$.
1629   *
1630   * <table border=1 cellpadding=10 cellspacing=0>
1631   * <caption align=top>Distribution Statistics</caption>
1632   * <tr><td>Mean</td><td>@f$\mu@f$</td></tr>
1633   * <tr><td>Variance</td><td>@f$\mu^2\frac{\lambda + \nu + 1}{\lambda\nu}@f$</td></tr>
1634   * <tr><td>Range</td><td>@f$[0, \infty)@f$</td></tr>
1635   * </table>
1636   */
1637  template<typename _RealType = double>
1638    class
1639    k_distribution
1640    {
1641      static_assert(std::is_floating_point<_RealType>::value,
1642		    "template argument not a floating point type");
1643
1644    public:
1645      /** The type of the range of the distribution. */
1646      typedef _RealType result_type;
1647      /** Parameter type. */
1648      struct param_type
1649      {
1650	typedef k_distribution<result_type> distribution_type;
1651
1652	param_type(result_type __lambda_val = result_type(1),
1653		   result_type __mu_val = result_type(1),
1654		   result_type __nu_val = result_type(1))
1655	: _M_lambda(__lambda_val), _M_mu(__mu_val), _M_nu(__nu_val)
1656	{
1657	  _GLIBCXX_DEBUG_ASSERT(_M_lambda > result_type(0));
1658	  _GLIBCXX_DEBUG_ASSERT(_M_mu > result_type(0));
1659	  _GLIBCXX_DEBUG_ASSERT(_M_nu > result_type(0));
1660	}
1661
1662	result_type
1663	lambda() const
1664	{ return _M_lambda; }
1665
1666	result_type
1667	mu() const
1668	{ return _M_mu; }
1669
1670	result_type
1671	nu() const
1672	{ return _M_nu; }
1673
1674	friend bool
1675	operator==(const param_type& __p1, const param_type& __p2)
1676	{ return __p1._M_lambda == __p2._M_lambda
1677	      && __p1._M_mu == __p2._M_mu
1678	      && __p1._M_nu == __p2._M_nu; }
1679
1680      private:
1681	void _M_initialize();
1682
1683	result_type _M_lambda;
1684	result_type _M_mu;
1685	result_type _M_nu;
1686      };
1687
1688      /**
1689       * @brief Constructors.
1690       */
1691      explicit
1692      k_distribution(result_type __lambda_val = result_type(1),
1693		     result_type __mu_val = result_type(1),
1694		     result_type __nu_val = result_type(1))
1695      : _M_param(__lambda_val, __mu_val, __nu_val),
1696	_M_gd1(__lambda_val, result_type(1) / __lambda_val),
1697	_M_gd2(__nu_val, __mu_val / __nu_val)
1698      { }
1699
1700      explicit
1701      k_distribution(const param_type& __p)
1702      : _M_param(__p),
1703	_M_gd1(__p.lambda(), result_type(1) / __p.lambda()),
1704	_M_gd2(__p.nu(), __p.mu() / __p.nu())
1705      { }
1706
1707      /**
1708       * @brief Resets the distribution state.
1709       */
1710      void
1711      reset()
1712      {
1713	_M_gd1.reset();
1714	_M_gd2.reset();
1715      }
1716
1717      /**
1718       * @brief Return the parameters of the distribution.
1719       */
1720      result_type
1721      lambda() const
1722      { return _M_param.lambda(); }
1723
1724      result_type
1725      mu() const
1726      { return _M_param.mu(); }
1727
1728      result_type
1729      nu() const
1730      { return _M_param.nu(); }
1731
1732      /**
1733       * @brief Returns the parameter set of the distribution.
1734       */
1735      param_type
1736      param() const
1737      { return _M_param; }
1738
1739      /**
1740       * @brief Sets the parameter set of the distribution.
1741       * @param __param The new parameter set of the distribution.
1742       */
1743      void
1744      param(const param_type& __param)
1745      { _M_param = __param; }
1746
1747      /**
1748       * @brief Returns the greatest lower bound value of the distribution.
1749       */
1750      result_type
1751      min() const
1752      { return result_type(0); }
1753
1754      /**
1755       * @brief Returns the least upper bound value of the distribution.
1756       */
1757      result_type
1758      max() const
1759      { return std::numeric_limits<result_type>::max(); }
1760
1761      /**
1762       * @brief Generating functions.
1763       */
1764      template<typename _UniformRandomNumberGenerator>
1765	result_type
1766	operator()(_UniformRandomNumberGenerator&);
1767
1768      template<typename _UniformRandomNumberGenerator>
1769	result_type
1770	operator()(_UniformRandomNumberGenerator&, const param_type&);
1771
1772      template<typename _ForwardIterator,
1773	       typename _UniformRandomNumberGenerator>
1774	void
1775	__generate(_ForwardIterator __f, _ForwardIterator __t,
1776		   _UniformRandomNumberGenerator& __urng)
1777	{ this->__generate(__f, __t, __urng, _M_param); }
1778
1779      template<typename _ForwardIterator,
1780	       typename _UniformRandomNumberGenerator>
1781	void
1782	__generate(_ForwardIterator __f, _ForwardIterator __t,
1783		   _UniformRandomNumberGenerator& __urng,
1784		   const param_type& __p)
1785	{ this->__generate_impl(__f, __t, __urng, __p); }
1786
1787      template<typename _UniformRandomNumberGenerator>
1788	void
1789	__generate(result_type* __f, result_type* __t,
1790		   _UniformRandomNumberGenerator& __urng,
1791		   const param_type& __p)
1792	{ this->__generate_impl(__f, __t, __urng, __p); }
1793
1794      /**
1795       * @brief Return true if two K distributions have
1796       *        the same parameters and the sequences that would
1797       *        be generated are equal.
1798       */
1799      friend bool
1800      operator==(const k_distribution& __d1,
1801		 const k_distribution& __d2)
1802      { return (__d1._M_param == __d2._M_param
1803		&& __d1._M_gd1 == __d2._M_gd1
1804		&& __d1._M_gd2 == __d2._M_gd2); }
1805
1806      /**
1807       * @brief Inserts a %k_distribution random number distribution
1808       * @p __x into the output stream @p __os.
1809       *
1810       * @param __os An output stream.
1811       * @param __x  A %k_distribution random number distribution.
1812       *
1813       * @returns The output stream with the state of @p __x inserted or in
1814       * an error state.
1815       */
1816      template<typename _RealType1, typename _CharT, typename _Traits>
1817	friend std::basic_ostream<_CharT, _Traits>&
1818	operator<<(std::basic_ostream<_CharT, _Traits>&,
1819		   const k_distribution<_RealType1>&);
1820
1821      /**
1822       * @brief Extracts a %k_distribution random number distribution
1823       * @p __x from the input stream @p __is.
1824       *
1825       * @param __is An input stream.
1826       * @param __x A %k_distribution random number
1827       *            generator engine.
1828       *
1829       * @returns The input stream with @p __x extracted or in an error state.
1830       */
1831      template<typename _RealType1, typename _CharT, typename _Traits>
1832	friend std::basic_istream<_CharT, _Traits>&
1833	operator>>(std::basic_istream<_CharT, _Traits>&,
1834		   k_distribution<_RealType1>&);
1835
1836    private:
1837      template<typename _ForwardIterator,
1838	       typename _UniformRandomNumberGenerator>
1839	void
1840	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1841			_UniformRandomNumberGenerator& __urng,
1842			const param_type& __p);
1843
1844      param_type _M_param;
1845
1846      std::gamma_distribution<result_type> _M_gd1;
1847      std::gamma_distribution<result_type> _M_gd2;
1848    };
1849
1850  /**
1851   * @brief Return true if two K distributions are not equal.
1852   */
1853  template<typename _RealType>
1854    inline bool
1855    operator!=(const k_distribution<_RealType>& __d1,
1856	       const k_distribution<_RealType>& __d2)
1857    { return !(__d1 == __d2); }
1858
1859
1860  /**
1861   * @brief An arcsine continuous distribution for random numbers.
1862   *
1863   * The formula for the arcsine probability density function is
1864   * @f[
1865   *     p(x|a,b) = \frac{1}{\pi \sqrt{(x - a)(b - x)}}
1866   * @f]
1867   * where @f$x >= a@f$ and @f$x <= b@f$.
1868   *
1869   * <table border=1 cellpadding=10 cellspacing=0>
1870   * <caption align=top>Distribution Statistics</caption>
1871   * <tr><td>Mean</td><td>@f$ (a + b) / 2 @f$</td></tr>
1872   * <tr><td>Variance</td><td>@f$ (b - a)^2 / 8 @f$</td></tr>
1873   * <tr><td>Range</td><td>@f$[a, b]@f$</td></tr>
1874   * </table>
1875   */
1876  template<typename _RealType = double>
1877    class
1878    arcsine_distribution
1879    {
1880      static_assert(std::is_floating_point<_RealType>::value,
1881		    "template argument not a floating point type");
1882
1883    public:
1884      /** The type of the range of the distribution. */
1885      typedef _RealType result_type;
1886      /** Parameter type. */
1887      struct param_type
1888      {
1889	typedef arcsine_distribution<result_type> distribution_type;
1890
1891	param_type(result_type __a = result_type(0),
1892		   result_type __b = result_type(1))
1893	: _M_a(__a), _M_b(__b)
1894	{
1895	  _GLIBCXX_DEBUG_ASSERT(_M_a <= _M_b);
1896	}
1897
1898	result_type
1899	a() const
1900	{ return _M_a; }
1901
1902	result_type
1903	b() const
1904	{ return _M_b; }
1905
1906	friend bool
1907	operator==(const param_type& __p1, const param_type& __p2)
1908	{ return __p1._M_a == __p2._M_a && __p1._M_b == __p2._M_b; }
1909
1910      private:
1911	void _M_initialize();
1912
1913	result_type _M_a;
1914	result_type _M_b;
1915      };
1916
1917      /**
1918       * @brief Constructors.
1919       */
1920      explicit
1921      arcsine_distribution(result_type __a = result_type(0),
1922			   result_type __b = result_type(1))
1923      : _M_param(__a, __b),
1924	_M_ud(-1.5707963267948966192313216916397514L,
1925	      +1.5707963267948966192313216916397514L)
1926      { }
1927
1928      explicit
1929      arcsine_distribution(const param_type& __p)
1930      : _M_param(__p),
1931	_M_ud(-1.5707963267948966192313216916397514L,
1932	      +1.5707963267948966192313216916397514L)
1933      { }
1934
1935      /**
1936       * @brief Resets the distribution state.
1937       */
1938      void
1939      reset()
1940      { _M_ud.reset(); }
1941
1942      /**
1943       * @brief Return the parameters of the distribution.
1944       */
1945      result_type
1946      a() const
1947      { return _M_param.a(); }
1948
1949      result_type
1950      b() const
1951      { return _M_param.b(); }
1952
1953      /**
1954       * @brief Returns the parameter set of the distribution.
1955       */
1956      param_type
1957      param() const
1958      { return _M_param; }
1959
1960      /**
1961       * @brief Sets the parameter set of the distribution.
1962       * @param __param The new parameter set of the distribution.
1963       */
1964      void
1965      param(const param_type& __param)
1966      { _M_param = __param; }
1967
1968      /**
1969       * @brief Returns the greatest lower bound value of the distribution.
1970       */
1971      result_type
1972      min() const
1973      { return this->a(); }
1974
1975      /**
1976       * @brief Returns the least upper bound value of the distribution.
1977       */
1978      result_type
1979      max() const
1980      { return this->b(); }
1981
1982      /**
1983       * @brief Generating functions.
1984       */
1985      template<typename _UniformRandomNumberGenerator>
1986	result_type
1987	operator()(_UniformRandomNumberGenerator& __urng)
1988	{
1989	  result_type __x = std::sin(this->_M_ud(__urng));
1990	  return (__x * (this->b() - this->a())
1991		  + this->a() + this->b()) / result_type(2);
1992	}
1993
1994      template<typename _UniformRandomNumberGenerator>
1995	result_type
1996	operator()(_UniformRandomNumberGenerator& __urng,
1997		   const param_type& __p)
1998	{
1999	  result_type __x = std::sin(this->_M_ud(__urng));
2000	  return (__x * (__p.b() - __p.a())
2001		  + __p.a() + __p.b()) / result_type(2);
2002	}
2003
2004      template<typename _ForwardIterator,
2005	       typename _UniformRandomNumberGenerator>
2006	void
2007	__generate(_ForwardIterator __f, _ForwardIterator __t,
2008		   _UniformRandomNumberGenerator& __urng)
2009	{ this->__generate(__f, __t, __urng, _M_param); }
2010
2011      template<typename _ForwardIterator,
2012	       typename _UniformRandomNumberGenerator>
2013	void
2014	__generate(_ForwardIterator __f, _ForwardIterator __t,
2015		   _UniformRandomNumberGenerator& __urng,
2016		   const param_type& __p)
2017	{ this->__generate_impl(__f, __t, __urng, __p); }
2018
2019      template<typename _UniformRandomNumberGenerator>
2020	void
2021	__generate(result_type* __f, result_type* __t,
2022		   _UniformRandomNumberGenerator& __urng,
2023		   const param_type& __p)
2024	{ this->__generate_impl(__f, __t, __urng, __p); }
2025
2026      /**
2027       * @brief Return true if two arcsine distributions have
2028       *        the same parameters and the sequences that would
2029       *        be generated are equal.
2030       */
2031      friend bool
2032      operator==(const arcsine_distribution& __d1,
2033		 const arcsine_distribution& __d2)
2034      { return (__d1._M_param == __d2._M_param
2035		&& __d1._M_ud == __d2._M_ud); }
2036
2037      /**
2038       * @brief Inserts a %arcsine_distribution random number distribution
2039       * @p __x into the output stream @p __os.
2040       *
2041       * @param __os An output stream.
2042       * @param __x  A %arcsine_distribution random number distribution.
2043       *
2044       * @returns The output stream with the state of @p __x inserted or in
2045       * an error state.
2046       */
2047      template<typename _RealType1, typename _CharT, typename _Traits>
2048	friend std::basic_ostream<_CharT, _Traits>&
2049	operator<<(std::basic_ostream<_CharT, _Traits>&,
2050		   const arcsine_distribution<_RealType1>&);
2051
2052      /**
2053       * @brief Extracts a %arcsine_distribution random number distribution
2054       * @p __x from the input stream @p __is.
2055       *
2056       * @param __is An input stream.
2057       * @param __x A %arcsine_distribution random number
2058       *            generator engine.
2059       *
2060       * @returns The input stream with @p __x extracted or in an error state.
2061       */
2062      template<typename _RealType1, typename _CharT, typename _Traits>
2063	friend std::basic_istream<_CharT, _Traits>&
2064	operator>>(std::basic_istream<_CharT, _Traits>&,
2065		   arcsine_distribution<_RealType1>&);
2066
2067    private:
2068      template<typename _ForwardIterator,
2069	       typename _UniformRandomNumberGenerator>
2070	void
2071	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2072			_UniformRandomNumberGenerator& __urng,
2073			const param_type& __p);
2074
2075      param_type _M_param;
2076
2077      std::uniform_real_distribution<result_type> _M_ud;
2078    };
2079
2080  /**
2081   * @brief Return true if two arcsine distributions are not equal.
2082   */
2083  template<typename _RealType>
2084    inline bool
2085    operator!=(const arcsine_distribution<_RealType>& __d1,
2086	       const arcsine_distribution<_RealType>& __d2)
2087    { return !(__d1 == __d2); }
2088
2089
2090  /**
2091   * @brief A Hoyt continuous distribution for random numbers.
2092   *
2093   * The formula for the Hoyt probability density function is
2094   * @f[
2095   *     p(x|q,\omega) = \frac{(1 + q^2)x}{q\omega}
2096   *                     \exp\left(-\frac{(1 + q^2)^2 x^2}{4 q^2 \omega}\right)
2097   *                       I_0\left(\frac{(1 - q^4) x^2}{4 q^2 \omega}\right)
2098   * @f]
2099   * where @f$I_0(z)@f$ is the modified Bessel function of the first kind
2100   * of order 0 and @f$0 < q < 1@f$.
2101   *
2102   * <table border=1 cellpadding=10 cellspacing=0>
2103   * <caption align=top>Distribution Statistics</caption>
2104   * <tr><td>Mean</td><td>@f$ \sqrt{\frac{2}{\pi}} \sqrt{\frac{\omega}{1 + q^2}}
2105   *                       E(1 - q^2) @f$</td></tr>
2106   * <tr><td>Variance</td><td>@f$ \omega \left(1 - \frac{2E^2(1 - q^2)}
2107   *                                      {\pi (1 + q^2)}\right) @f$</td></tr>
2108   * <tr><td>Range</td><td>@f$[0, \infty)@f$</td></tr>
2109   * </table>
2110   * where @f$E(x)@f$ is the elliptic function of the second kind.
2111   */
2112  template<typename _RealType = double>
2113    class
2114    hoyt_distribution
2115    {
2116      static_assert(std::is_floating_point<_RealType>::value,
2117		    "template argument not a floating point type");
2118
2119    public:
2120      /** The type of the range of the distribution. */
2121      typedef _RealType result_type;
2122      /** Parameter type. */
2123      struct param_type
2124      {
2125	typedef hoyt_distribution<result_type> distribution_type;
2126
2127	param_type(result_type __q = result_type(0.5L),
2128		   result_type __omega = result_type(1))
2129	: _M_q(__q), _M_omega(__omega)
2130	{
2131	  _GLIBCXX_DEBUG_ASSERT(_M_q > result_type(0));
2132	  _GLIBCXX_DEBUG_ASSERT(_M_q < result_type(1));
2133	}
2134
2135	result_type
2136	q() const
2137	{ return _M_q; }
2138
2139	result_type
2140	omega() const
2141	{ return _M_omega; }
2142
2143	friend bool
2144	operator==(const param_type& __p1, const param_type& __p2)
2145	{ return __p1._M_q == __p2._M_q
2146	      && __p1._M_omega == __p2._M_omega; }
2147
2148      private:
2149	void _M_initialize();
2150
2151	result_type _M_q;
2152	result_type _M_omega;
2153      };
2154
2155      /**
2156       * @brief Constructors.
2157       */
2158      explicit
2159      hoyt_distribution(result_type __q = result_type(0.5L),
2160			result_type __omega = result_type(1))
2161      : _M_param(__q, __omega),
2162	_M_ad(result_type(0.5L) * (result_type(1) + __q * __q),
2163	      result_type(0.5L) * (result_type(1) + __q * __q)
2164				/ (__q * __q)),
2165	_M_ed(result_type(1))
2166      { }
2167
2168      explicit
2169      hoyt_distribution(const param_type& __p)
2170      : _M_param(__p),
2171	_M_ad(result_type(0.5L) * (result_type(1) + __p.q() * __p.q()),
2172	      result_type(0.5L) * (result_type(1) + __p.q() * __p.q())
2173				/ (__p.q() * __p.q())),
2174	_M_ed(result_type(1))
2175      { }
2176
2177      /**
2178       * @brief Resets the distribution state.
2179       */
2180      void
2181      reset()
2182      {
2183	_M_ad.reset();
2184	_M_ed.reset();
2185      }
2186
2187      /**
2188       * @brief Return the parameters of the distribution.
2189       */
2190      result_type
2191      q() const
2192      { return _M_param.q(); }
2193
2194      result_type
2195      omega() const
2196      { return _M_param.omega(); }
2197
2198      /**
2199       * @brief Returns the parameter set of the distribution.
2200       */
2201      param_type
2202      param() const
2203      { return _M_param; }
2204
2205      /**
2206       * @brief Sets the parameter set of the distribution.
2207       * @param __param The new parameter set of the distribution.
2208       */
2209      void
2210      param(const param_type& __param)
2211      { _M_param = __param; }
2212
2213      /**
2214       * @brief Returns the greatest lower bound value of the distribution.
2215       */
2216      result_type
2217      min() const
2218      { return result_type(0); }
2219
2220      /**
2221       * @brief Returns the least upper bound value of the distribution.
2222       */
2223      result_type
2224      max() const
2225      { return std::numeric_limits<result_type>::max(); }
2226
2227      /**
2228       * @brief Generating functions.
2229       */
2230      template<typename _UniformRandomNumberGenerator>
2231	result_type
2232	operator()(_UniformRandomNumberGenerator& __urng);
2233
2234      template<typename _UniformRandomNumberGenerator>
2235	result_type
2236	operator()(_UniformRandomNumberGenerator& __urng,
2237		   const param_type& __p);
2238
2239      template<typename _ForwardIterator,
2240	       typename _UniformRandomNumberGenerator>
2241	void
2242	__generate(_ForwardIterator __f, _ForwardIterator __t,
2243		   _UniformRandomNumberGenerator& __urng)
2244	{ this->__generate(__f, __t, __urng, _M_param); }
2245
2246      template<typename _ForwardIterator,
2247	       typename _UniformRandomNumberGenerator>
2248	void
2249	__generate(_ForwardIterator __f, _ForwardIterator __t,
2250		   _UniformRandomNumberGenerator& __urng,
2251		   const param_type& __p)
2252	{ this->__generate_impl(__f, __t, __urng, __p); }
2253
2254      template<typename _UniformRandomNumberGenerator>
2255	void
2256	__generate(result_type* __f, result_type* __t,
2257		   _UniformRandomNumberGenerator& __urng,
2258		   const param_type& __p)
2259	{ this->__generate_impl(__f, __t, __urng, __p); }
2260
2261      /**
2262       * @brief Return true if two Hoyt distributions have
2263       *        the same parameters and the sequences that would
2264       *        be generated are equal.
2265       */
2266      friend bool
2267      operator==(const hoyt_distribution& __d1,
2268		 const hoyt_distribution& __d2)
2269      { return (__d1._M_param == __d2._M_param
2270		&& __d1._M_ad == __d2._M_ad
2271		&& __d1._M_ed == __d2._M_ed); }
2272
2273      /**
2274       * @brief Inserts a %hoyt_distribution random number distribution
2275       * @p __x into the output stream @p __os.
2276       *
2277       * @param __os An output stream.
2278       * @param __x  A %hoyt_distribution random number distribution.
2279       *
2280       * @returns The output stream with the state of @p __x inserted or in
2281       * an error state.
2282       */
2283      template<typename _RealType1, typename _CharT, typename _Traits>
2284	friend std::basic_ostream<_CharT, _Traits>&
2285	operator<<(std::basic_ostream<_CharT, _Traits>&,
2286		   const hoyt_distribution<_RealType1>&);
2287
2288      /**
2289       * @brief Extracts a %hoyt_distribution random number distribution
2290       * @p __x from the input stream @p __is.
2291       *
2292       * @param __is An input stream.
2293       * @param __x A %hoyt_distribution random number
2294       *            generator engine.
2295       *
2296       * @returns The input stream with @p __x extracted or in an error state.
2297       */
2298      template<typename _RealType1, typename _CharT, typename _Traits>
2299	friend std::basic_istream<_CharT, _Traits>&
2300	operator>>(std::basic_istream<_CharT, _Traits>&,
2301		   hoyt_distribution<_RealType1>&);
2302
2303    private:
2304      template<typename _ForwardIterator,
2305	       typename _UniformRandomNumberGenerator>
2306	void
2307	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2308			_UniformRandomNumberGenerator& __urng,
2309			const param_type& __p);
2310
2311      param_type _M_param;
2312
2313      __gnu_cxx::arcsine_distribution<result_type> _M_ad;
2314      std::exponential_distribution<result_type> _M_ed;
2315    };
2316
2317  /**
2318   * @brief Return true if two Hoyt distributions are not equal.
2319   */
2320  template<typename _RealType>
2321    inline bool
2322    operator!=(const hoyt_distribution<_RealType>& __d1,
2323	       const hoyt_distribution<_RealType>& __d2)
2324    { return !(__d1 == __d2); }
2325
2326
2327  /**
2328   * @brief A triangular distribution for random numbers.
2329   *
2330   * The formula for the triangular probability density function is
2331   * @f[
2332   *                  / 0                          for x < a
2333   *     p(x|a,b,c) = | \frac{2(x-a)}{(c-a)(b-a)}  for a <= x <= b
2334   *                  | \frac{2(c-x)}{(c-a)(c-b)}  for b < x <= c
2335   *                  \ 0                          for c < x
2336   * @f]
2337   *
2338   * <table border=1 cellpadding=10 cellspacing=0>
2339   * <caption align=top>Distribution Statistics</caption>
2340   * <tr><td>Mean</td><td>@f$ \frac{a+b+c}{2} @f$</td></tr>
2341   * <tr><td>Variance</td><td>@f$ \frac{a^2+b^2+c^2-ab-ac-bc}
2342   *                                   {18}@f$</td></tr>
2343   * <tr><td>Range</td><td>@f$[a, c]@f$</td></tr>
2344   * </table>
2345   */
2346  template<typename _RealType = double>
2347    class triangular_distribution
2348    {
2349      static_assert(std::is_floating_point<_RealType>::value,
2350		    "template argument not a floating point type");
2351
2352    public:
2353      /** The type of the range of the distribution. */
2354      typedef _RealType result_type;
2355      /** Parameter type. */
2356      struct param_type
2357      {
2358	friend class triangular_distribution<_RealType>;
2359
2360	explicit
2361	param_type(_RealType __a = _RealType(0),
2362		   _RealType __b = _RealType(0.5),
2363		   _RealType __c = _RealType(1))
2364	: _M_a(__a), _M_b(__b), _M_c(__c)
2365	{
2366	  _GLIBCXX_DEBUG_ASSERT(_M_a <= _M_b);
2367	  _GLIBCXX_DEBUG_ASSERT(_M_b <= _M_c);
2368	  _GLIBCXX_DEBUG_ASSERT(_M_a < _M_c);
2369
2370	  _M_r_ab = (_M_b - _M_a) / (_M_c - _M_a);
2371	  _M_f_ab_ac = (_M_b - _M_a) * (_M_c - _M_a);
2372	  _M_f_bc_ac = (_M_c - _M_b) * (_M_c - _M_a);
2373	}
2374
2375	_RealType
2376	a() const
2377	{ return _M_a; }
2378
2379	_RealType
2380	b() const
2381	{ return _M_b; }
2382
2383	_RealType
2384	c() const
2385	{ return _M_c; }
2386
2387	friend bool
2388	operator==(const param_type& __p1, const param_type& __p2)
2389	{ return (__p1._M_a == __p2._M_a && __p1._M_b == __p2._M_b
2390		  && __p1._M_c == __p2._M_c); }
2391
2392      private:
2393
2394	_RealType _M_a;
2395	_RealType _M_b;
2396	_RealType _M_c;
2397	_RealType _M_r_ab;
2398	_RealType _M_f_ab_ac;
2399	_RealType _M_f_bc_ac;
2400      };
2401
2402      /**
2403       * @brief Constructs a triangle distribution with parameters
2404       * @f$ a @f$, @f$ b @f$ and @f$ c @f$.
2405       */
2406      explicit
2407      triangular_distribution(result_type __a = result_type(0),
2408			      result_type __b = result_type(0.5),
2409			      result_type __c = result_type(1))
2410      : _M_param(__a, __b, __c)
2411      { }
2412
2413      explicit
2414      triangular_distribution(const param_type& __p)
2415      : _M_param(__p)
2416      { }
2417
2418      /**
2419       * @brief Resets the distribution state.
2420       */
2421      void
2422      reset()
2423      { }
2424
2425      /**
2426       * @brief Returns the @f$ a @f$ of the distribution.
2427       */
2428      result_type
2429      a() const
2430      { return _M_param.a(); }
2431
2432      /**
2433       * @brief Returns the @f$ b @f$ of the distribution.
2434       */
2435      result_type
2436      b() const
2437      { return _M_param.b(); }
2438
2439      /**
2440       * @brief Returns the @f$ c @f$ of the distribution.
2441       */
2442      result_type
2443      c() const
2444      { return _M_param.c(); }
2445
2446      /**
2447       * @brief Returns the parameter set of the distribution.
2448       */
2449      param_type
2450      param() const
2451      { return _M_param; }
2452
2453      /**
2454       * @brief Sets the parameter set of the distribution.
2455       * @param __param The new parameter set of the distribution.
2456       */
2457      void
2458      param(const param_type& __param)
2459      { _M_param = __param; }
2460
2461      /**
2462       * @brief Returns the greatest lower bound value of the distribution.
2463       */
2464      result_type
2465      min() const
2466      { return _M_param._M_a; }
2467
2468      /**
2469       * @brief Returns the least upper bound value of the distribution.
2470       */
2471      result_type
2472      max() const
2473      { return _M_param._M_c; }
2474
2475      /**
2476       * @brief Generating functions.
2477       */
2478      template<typename _UniformRandomNumberGenerator>
2479	result_type
2480	operator()(_UniformRandomNumberGenerator& __urng)
2481	{ return this->operator()(__urng, _M_param); }
2482
2483      template<typename _UniformRandomNumberGenerator>
2484	result_type
2485	operator()(_UniformRandomNumberGenerator& __urng,
2486		   const param_type& __p)
2487	{
2488	  std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2489	    __aurng(__urng);
2490	  result_type __rnd = __aurng();
2491	  if (__rnd <= __p._M_r_ab)
2492	    return __p.a() + std::sqrt(__rnd * __p._M_f_ab_ac);
2493	  else
2494	    return __p.c() - std::sqrt((result_type(1) - __rnd)
2495				       * __p._M_f_bc_ac);
2496	}
2497
2498      template<typename _ForwardIterator,
2499	       typename _UniformRandomNumberGenerator>
2500	void
2501	__generate(_ForwardIterator __f, _ForwardIterator __t,
2502		   _UniformRandomNumberGenerator& __urng)
2503	{ this->__generate(__f, __t, __urng, _M_param); }
2504
2505      template<typename _ForwardIterator,
2506	       typename _UniformRandomNumberGenerator>
2507	void
2508	__generate(_ForwardIterator __f, _ForwardIterator __t,
2509		   _UniformRandomNumberGenerator& __urng,
2510		   const param_type& __p)
2511	{ this->__generate_impl(__f, __t, __urng, __p); }
2512
2513      template<typename _UniformRandomNumberGenerator>
2514	void
2515	__generate(result_type* __f, result_type* __t,
2516		   _UniformRandomNumberGenerator& __urng,
2517		   const param_type& __p)
2518	{ this->__generate_impl(__f, __t, __urng, __p); }
2519
2520      /**
2521       * @brief Return true if two triangle distributions have the same
2522       *        parameters and the sequences that would be generated
2523       *        are equal.
2524       */
2525      friend bool
2526      operator==(const triangular_distribution& __d1,
2527		 const triangular_distribution& __d2)
2528      { return __d1._M_param == __d2._M_param; }
2529
2530      /**
2531       * @brief Inserts a %triangular_distribution random number distribution
2532       * @p __x into the output stream @p __os.
2533       *
2534       * @param __os An output stream.
2535       * @param __x  A %triangular_distribution random number distribution.
2536       *
2537       * @returns The output stream with the state of @p __x inserted or in
2538       * an error state.
2539       */
2540      template<typename _RealType1, typename _CharT, typename _Traits>
2541	friend std::basic_ostream<_CharT, _Traits>&
2542	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2543		   const __gnu_cxx::triangular_distribution<_RealType1>& __x);
2544
2545      /**
2546       * @brief Extracts a %triangular_distribution random number distribution
2547       * @p __x from the input stream @p __is.
2548       *
2549       * @param __is An input stream.
2550       * @param __x  A %triangular_distribution random number generator engine.
2551       *
2552       * @returns The input stream with @p __x extracted or in an error state.
2553       */
2554      template<typename _RealType1, typename _CharT, typename _Traits>
2555	friend std::basic_istream<_CharT, _Traits>&
2556	operator>>(std::basic_istream<_CharT, _Traits>& __is,
2557		   __gnu_cxx::triangular_distribution<_RealType1>& __x);
2558
2559    private:
2560      template<typename _ForwardIterator,
2561	       typename _UniformRandomNumberGenerator>
2562	void
2563	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2564			_UniformRandomNumberGenerator& __urng,
2565			const param_type& __p);
2566
2567      param_type _M_param;
2568    };
2569
2570  /**
2571   * @brief Return true if two triangle distributions are different.
2572   */
2573  template<typename _RealType>
2574    inline bool
2575    operator!=(const __gnu_cxx::triangular_distribution<_RealType>& __d1,
2576	       const __gnu_cxx::triangular_distribution<_RealType>& __d2)
2577   { return !(__d1 == __d2); }
2578
2579
2580  /**
2581   * @brief A von Mises distribution for random numbers.
2582   *
2583   * The formula for the von Mises probability density function is
2584   * @f[
2585   *     p(x|\mu,\kappa) = \frac{e^{\kappa \cos(x-\mu)}}
2586   *                            {2\pi I_0(\kappa)}
2587   * @f]
2588   *
2589   * The generating functions use the method according to:
2590   *
2591   * D. J. Best and N. I. Fisher, 1979. "Efficient Simulation of the
2592   * von Mises Distribution", Journal of the Royal Statistical Society.
2593   * Series C (Applied Statistics), Vol. 28, No. 2, pp. 152-157.
2594   *
2595   * <table border=1 cellpadding=10 cellspacing=0>
2596   * <caption align=top>Distribution Statistics</caption>
2597   * <tr><td>Mean</td><td>@f$ \mu @f$</td></tr>
2598   * <tr><td>Variance</td><td>@f$ 1-I_1(\kappa)/I_0(\kappa) @f$</td></tr>
2599   * <tr><td>Range</td><td>@f$[-\pi, \pi]@f$</td></tr>
2600   * </table>
2601   */
2602  template<typename _RealType = double>
2603    class von_mises_distribution
2604    {
2605      static_assert(std::is_floating_point<_RealType>::value,
2606		    "template argument not a floating point type");
2607
2608    public:
2609      /** The type of the range of the distribution. */
2610      typedef _RealType result_type;
2611      /** Parameter type. */
2612      struct param_type
2613      {
2614	friend class von_mises_distribution<_RealType>;
2615
2616	explicit
2617	param_type(_RealType __mu = _RealType(0),
2618		   _RealType __kappa = _RealType(1))
2619	: _M_mu(__mu), _M_kappa(__kappa)
2620	{
2621	  const _RealType __pi = __gnu_cxx::__math_constants<_RealType>::__pi;
2622	  _GLIBCXX_DEBUG_ASSERT(_M_mu >= -__pi && _M_mu <= __pi);
2623	  _GLIBCXX_DEBUG_ASSERT(_M_kappa >= _RealType(0));
2624
2625	  auto __tau = std::sqrt(_RealType(4) * _M_kappa * _M_kappa
2626				 + _RealType(1)) + _RealType(1);
2627	  auto __rho = ((__tau - std::sqrt(_RealType(2) * __tau))
2628			/ (_RealType(2) * _M_kappa));
2629	  _M_r = (_RealType(1) + __rho * __rho) / (_RealType(2) * __rho);
2630	}
2631
2632	_RealType
2633	mu() const
2634	{ return _M_mu; }
2635
2636	_RealType
2637	kappa() const
2638	{ return _M_kappa; }
2639
2640	friend bool
2641	operator==(const param_type& __p1, const param_type& __p2)
2642	{ return (__p1._M_mu == __p2._M_mu
2643		  && __p1._M_kappa == __p2._M_kappa); }
2644
2645      private:
2646	_RealType _M_mu;
2647	_RealType _M_kappa;
2648	_RealType _M_r;
2649      };
2650
2651      /**
2652       * @brief Constructs a von Mises distribution with parameters
2653       * @f$\mu@f$ and @f$\kappa@f$.
2654       */
2655      explicit
2656      von_mises_distribution(result_type __mu = result_type(0),
2657			     result_type __kappa = result_type(1))
2658	: _M_param(__mu, __kappa)
2659      { }
2660
2661      explicit
2662      von_mises_distribution(const param_type& __p)
2663      : _M_param(__p)
2664      { }
2665
2666      /**
2667       * @brief Resets the distribution state.
2668       */
2669      void
2670      reset()
2671      { }
2672
2673      /**
2674       * @brief Returns the @f$ \mu @f$ of the distribution.
2675       */
2676      result_type
2677      mu() const
2678      { return _M_param.mu(); }
2679
2680      /**
2681       * @brief Returns the @f$ \kappa @f$ of the distribution.
2682       */
2683      result_type
2684      kappa() const
2685      { return _M_param.kappa(); }
2686
2687      /**
2688       * @brief Returns the parameter set of the distribution.
2689       */
2690      param_type
2691      param() const
2692      { return _M_param; }
2693
2694      /**
2695       * @brief Sets the parameter set of the distribution.
2696       * @param __param The new parameter set of the distribution.
2697       */
2698      void
2699      param(const param_type& __param)
2700      { _M_param = __param; }
2701
2702      /**
2703       * @brief Returns the greatest lower bound value of the distribution.
2704       */
2705      result_type
2706      min() const
2707      {
2708	return -__gnu_cxx::__math_constants<result_type>::__pi;
2709      }
2710
2711      /**
2712       * @brief Returns the least upper bound value of the distribution.
2713       */
2714      result_type
2715      max() const
2716      {
2717	return __gnu_cxx::__math_constants<result_type>::__pi;
2718      }
2719
2720      /**
2721       * @brief Generating functions.
2722       */
2723      template<typename _UniformRandomNumberGenerator>
2724	result_type
2725	operator()(_UniformRandomNumberGenerator& __urng)
2726	{ return this->operator()(__urng, _M_param); }
2727
2728      template<typename _UniformRandomNumberGenerator>
2729	result_type
2730	operator()(_UniformRandomNumberGenerator& __urng,
2731		   const param_type& __p)
2732	{
2733	  const result_type __pi
2734	    = __gnu_cxx::__math_constants<result_type>::__pi;
2735	  std::__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
2736	    __aurng(__urng);
2737
2738	  result_type __f;
2739	  while (1)
2740	    {
2741	      result_type __rnd = std::cos(__pi * __aurng());
2742	      __f = (result_type(1) + __p._M_r * __rnd) / (__p._M_r + __rnd);
2743	      result_type __c = __p._M_kappa * (__p._M_r - __f);
2744
2745	      result_type __rnd2 = __aurng();
2746	      if (__c * (result_type(2) - __c) > __rnd2)
2747		break;
2748	      if (std::log(__c / __rnd2) >= __c - result_type(1))
2749		break;
2750	    }
2751
2752	  result_type __res = std::acos(__f);
2753#if _GLIBCXX_USE_C99_MATH_TR1
2754	  __res = std::copysign(__res, __aurng() - result_type(0.5));
2755#else
2756	  if (__aurng() < result_type(0.5))
2757	    __res = -__res;
2758#endif
2759	  __res += __p._M_mu;
2760	  if (__res > __pi)
2761	    __res -= result_type(2) * __pi;
2762	  else if (__res < -__pi)
2763	    __res += result_type(2) * __pi;
2764	  return __res;
2765	}
2766
2767      template<typename _ForwardIterator,
2768	       typename _UniformRandomNumberGenerator>
2769	void
2770	__generate(_ForwardIterator __f, _ForwardIterator __t,
2771		   _UniformRandomNumberGenerator& __urng)
2772	{ this->__generate(__f, __t, __urng, _M_param); }
2773
2774      template<typename _ForwardIterator,
2775	       typename _UniformRandomNumberGenerator>
2776	void
2777	__generate(_ForwardIterator __f, _ForwardIterator __t,
2778		   _UniformRandomNumberGenerator& __urng,
2779		   const param_type& __p)
2780	{ this->__generate_impl(__f, __t, __urng, __p); }
2781
2782      template<typename _UniformRandomNumberGenerator>
2783	void
2784	__generate(result_type* __f, result_type* __t,
2785		   _UniformRandomNumberGenerator& __urng,
2786		   const param_type& __p)
2787	{ this->__generate_impl(__f, __t, __urng, __p); }
2788
2789      /**
2790       * @brief Return true if two von Mises distributions have the same
2791       *        parameters and the sequences that would be generated
2792       *        are equal.
2793       */
2794      friend bool
2795      operator==(const von_mises_distribution& __d1,
2796		 const von_mises_distribution& __d2)
2797      { return __d1._M_param == __d2._M_param; }
2798
2799      /**
2800       * @brief Inserts a %von_mises_distribution random number distribution
2801       * @p __x into the output stream @p __os.
2802       *
2803       * @param __os An output stream.
2804       * @param __x  A %von_mises_distribution random number distribution.
2805       *
2806       * @returns The output stream with the state of @p __x inserted or in
2807       * an error state.
2808       */
2809      template<typename _RealType1, typename _CharT, typename _Traits>
2810	friend std::basic_ostream<_CharT, _Traits>&
2811	operator<<(std::basic_ostream<_CharT, _Traits>& __os,
2812		   const __gnu_cxx::von_mises_distribution<_RealType1>& __x);
2813
2814      /**
2815       * @brief Extracts a %von_mises_distribution random number distribution
2816       * @p __x from the input stream @p __is.
2817       *
2818       * @param __is An input stream.
2819       * @param __x  A %von_mises_distribution random number generator engine.
2820       *
2821       * @returns The input stream with @p __x extracted or in an error state.
2822       */
2823      template<typename _RealType1, typename _CharT, typename _Traits>
2824	friend std::basic_istream<_CharT, _Traits>&
2825	operator>>(std::basic_istream<_CharT, _Traits>& __is,
2826		   __gnu_cxx::von_mises_distribution<_RealType1>& __x);
2827
2828    private:
2829      template<typename _ForwardIterator,
2830	       typename _UniformRandomNumberGenerator>
2831	void
2832	__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
2833			_UniformRandomNumberGenerator& __urng,
2834			const param_type& __p);
2835
2836      param_type _M_param;
2837    };
2838
2839  /**
2840   * @brief Return true if two von Mises distributions are different.
2841   */
2842  template<typename _RealType>
2843    inline bool
2844    operator!=(const __gnu_cxx::von_mises_distribution<_RealType>& __d1,
2845	       const __gnu_cxx::von_mises_distribution<_RealType>& __d2)
2846   { return !(__d1 == __d2); }
2847
2848_GLIBCXX_END_NAMESPACE_VERSION
2849} // namespace __gnu_cxx
2850
2851#include "ext/opt_random.h"
2852#include "random.tcc"
2853
2854#endif // _GLIBCXX_USE_C99_STDINT_TR1
2855
2856#endif // C++11
2857
2858#endif // _EXT_RANDOM
2859