Subversion Repositories Kolibri OS

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
5134 serge 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
// .
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 
36
#else
37
 
38
#include 
39
#include 
40
#include 
41
#ifdef __SSE2__
42
# include 
43
#endif
44
 
45
#ifdef _GLIBCXX_USE_C99_STDINT_TR1
46
 
47
namespace __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
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
98
	std::enable_if
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
109
	typename std::enable_if::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::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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
512
	result_type
513
	operator()(_UniformRandomNumberGenerator& __urng)
514
	{ return this->operator()(__urng, _M_param); }
515
 
516
      template
517
	result_type
518
	operator()(_UniformRandomNumberGenerator& __urng,
519
		   const param_type& __p);
520
 
521
      template
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
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
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
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
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
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
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
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
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 
708
	  void _M_init_full(_InputIterator1 __meanbegin,
709
			    _InputIterator1 __meanend,
710
			    _InputIterator2 __varcovbegin,
711
			    _InputIterator2 __varcovend);
712
	template 
713
	  void _M_init_lower(_InputIterator1 __meanbegin,
714
			     _InputIterator1 __meanend,
715
			     _InputIterator2 __varcovbegin,
716
			     _InputIterator2 __varcovend);
717
	template 
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
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
810
	result_type
811
	operator()(_UniformRandomNumberGenerator& __urng)
812
	{ return this->operator()(__urng, _M_param); }
813
 
814
      template
815
	result_type
816
	operator()(_UniformRandomNumberGenerator& __urng,
817
		   const param_type& __p);
818
 
819
      template
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
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
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
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
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
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
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
   * 
921
   * 
Distribution Statistics
922
   * 
Mean@f$\sqrt{\pi/2}L_{1/2}(-\nu^2/2\sigma^2)@f$
923
   * 
Variance@f$2\sigma^2 + \nu^2
924
   *                   + (\pi\sigma^2/2)L^2_{1/2}(-\nu^2/2\sigma^2)@f$
925
   * 
Range@f$[0, \infty)@f$
926
   * 
927
   * where @f$L_{1/2}(x)@f$ is the Laguerre polynomial of order 1/2.
928
   */
929
  template
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 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::max(); }
1038
 
1039
      /**
1040
       * @brief Generating functions.
1041
       */
1042
      template
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
1056
	result_type
1057
	operator()(_UniformRandomNumberGenerator& __urng,
1058
		   const param_type& __p)
1059
	{
1060
	  typename std::normal_distribution::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
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
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
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
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
1131
	friend std::basic_istream<_CharT, _Traits>&
1132
	operator>>(std::basic_istream<_CharT, _Traits>&,
1133
		   rice_distribution<_RealType1>&);
1134
 
1135
    private:
1136
      template
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 _M_ndx;
1146
      std::normal_distribution _M_ndy;
1147
    };
1148
 
1149
  /**
1150
   * @brief Return true if two Rice distributions are not equal.
1151
   */
1152
  template
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
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 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::max(); }
1275
 
1276
      /**
1277
       * @brief Generating functions.
1278
       */
1279
      template
1280
	result_type
1281
	operator()(_UniformRandomNumberGenerator& __urng)
1282
	{ return std::sqrt(this->_M_gd(__urng)); }
1283
 
1284
      template
1285
	result_type
1286
	operator()(_UniformRandomNumberGenerator& __urng,
1287
		   const param_type& __p)
1288
	{
1289
	  typename std::gamma_distribution::param_type
1290
	    __pg(__p.mu(), __p.omega() / __p.mu());
1291
	  return std::sqrt(this->_M_gd(__pg, __urng));
1292
	}
1293
 
1294
      template
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
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
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
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
1353
	friend std::basic_istream<_CharT, _Traits>&
1354
	operator>>(std::basic_istream<_CharT, _Traits>&,
1355
		   nakagami_distribution<_RealType1>&);
1356
 
1357
    private:
1358
      template
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 _M_gd;
1368
    };
1369
 
1370
  /**
1371
   * @brief Return true if two Nakagami distributions are not equal.
1372
   */
1373
  template
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
   * 
1395
   * 
Distribution Statistics
1396
   * 
Mean@f$\alpha \mu / (\alpha - 1)@f$
1397
   *              for @f$\alpha > 1@f$
1398
   * 
Variance@f$\alpha \mu^2 / [(\alpha - 1)^2(\alpha - 2)]@f$
1399
   *              for @f$\alpha > 2@f$
1400
   * 
Range@f$[\mu, \infty)@f$
1401
   * 
1402
   */
1403
  template
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 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::max(); }
1509
 
1510
      /**
1511
       * @brief Generating functions.
1512
       */
1513
      template
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
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
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
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
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
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
1589
	friend std::basic_istream<_CharT, _Traits>&
1590
	operator>>(std::basic_istream<_CharT, _Traits>&,
1591
		   pareto_distribution<_RealType1>&);
1592
 
1593
    private:
1594
      template
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 _M_ud;
1604
    };
1605
 
1606
  /**
1607
   * @brief Return true if two Pareto distributions are not equal.
1608
   */
1609
  template
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
   * 
1631
   * 
Distribution Statistics
1632
   * 
Mean@f$\mu@f$
1633
   * 
Variance@f$\mu^2\frac{\lambda + \nu + 1}{\lambda\nu}@f$
1634
   * 
Range@f$[0, \infty)@f$
1635
   * 
1636
   */
1637
  template
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 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::max(); }
1760
 
1761
      /**
1762
       * @brief Generating functions.
1763
       */
1764
      template
1765
	result_type
1766
	operator()(_UniformRandomNumberGenerator&);
1767
 
1768
      template
1769
	result_type
1770
	operator()(_UniformRandomNumberGenerator&, const param_type&);
1771
 
1772
      template
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
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
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
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
1832
	friend std::basic_istream<_CharT, _Traits>&
1833
	operator>>(std::basic_istream<_CharT, _Traits>&,
1834
		   k_distribution<_RealType1>&);
1835
 
1836
    private:
1837
      template
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 _M_gd1;
1847
      std::gamma_distribution _M_gd2;
1848
    };
1849
 
1850
  /**
1851
   * @brief Return true if two K distributions are not equal.
1852
   */
1853
  template
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
   * 
1870
   * 
Distribution Statistics
1871
   * 
Mean@f$ (a + b) / 2 @f$
1872
   * 
Variance@f$ (b - a)^2 / 8 @f$
1873
   * 
Range@f$[a, b]@f$
1874
   * 
1875
   */
1876
  template
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 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
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
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
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
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
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
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
2063
	friend std::basic_istream<_CharT, _Traits>&
2064
	operator>>(std::basic_istream<_CharT, _Traits>&,
2065
		   arcsine_distribution<_RealType1>&);
2066
 
2067
    private:
2068
      template
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 _M_ud;
2078
    };
2079
 
2080
  /**
2081
   * @brief Return true if two arcsine distributions are not equal.
2082
   */
2083
  template
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
   * 
2103
   * 
Distribution Statistics
2104
   * 
Mean@f$ \sqrt{\frac{2}{\pi}} \sqrt{\frac{\omega}{1 + q^2}}
2105
   *                       E(1 - q^2) @f$
2106
   * 
Variance@f$ \omega \left(1 - \frac{2E^2(1 - q^2)}
2107
   *                                      {\pi (1 + q^2)}\right) @f$
2108
   * 
Range@f$[0, \infty)@f$
2109
   * 
2110
   * where @f$E(x)@f$ is the elliptic function of the second kind.
2111
   */
2112
  template
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 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::max(); }
2226
 
2227
      /**
2228
       * @brief Generating functions.
2229
       */
2230
      template
2231
	result_type
2232
	operator()(_UniformRandomNumberGenerator& __urng);
2233
 
2234
      template
2235
	result_type
2236
	operator()(_UniformRandomNumberGenerator& __urng,
2237
		   const param_type& __p);
2238
 
2239
      template
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
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
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
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
2299
	friend std::basic_istream<_CharT, _Traits>&
2300
	operator>>(std::basic_istream<_CharT, _Traits>&,
2301
		   hoyt_distribution<_RealType1>&);
2302
 
2303
    private:
2304
      template
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 _M_ad;
2314
      std::exponential_distribution _M_ed;
2315
    };
2316
 
2317
  /**
2318
   * @brief Return true if two Hoyt distributions are not equal.
2319
   */
2320
  template
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
   * 
2339
   * 
Distribution Statistics
2340
   * 
Mean@f$ \frac{a+b+c}{2} @f$
2341
   * 
Variance@f$ \frac{a^2+b^2+c^2-ab-ac-bc}
2342
   *                                   {18}@f$
2343
   * 
Range@f$[a, c]@f$
2344
   * 
2345
   */
2346
  template
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
2479
	result_type
2480
	operator()(_UniformRandomNumberGenerator& __urng)
2481
	{ return this->operator()(__urng, _M_param); }
2482
 
2483
      template
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
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
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
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
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
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
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
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
   * 
2596
   * 
Distribution Statistics
2597
   * 
Mean@f$ \mu @f$
2598
   * 
Variance@f$ 1-I_1(\kappa)/I_0(\kappa) @f$
2599
   * 
Range@f$[-\pi, \pi]@f$
2600
   * 
2601
   */
2602
  template
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::__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::__pi;
2718
      }
2719
 
2720
      /**
2721
       * @brief Generating functions.
2722
       */
2723
      template
2724
	result_type
2725
	operator()(_UniformRandomNumberGenerator& __urng)
2726
	{ return this->operator()(__urng, _M_param); }
2727
 
2728
      template
2729
	result_type
2730
	operator()(_UniformRandomNumberGenerator& __urng,
2731
		   const param_type& __p)
2732
	{
2733
	  const result_type __pi
2734
	    = __gnu_cxx::__math_constants::__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
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
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
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
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
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
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
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