| blob | c7321a996238c36dcf133d47076dd13a0fbc4a32 |
1 // Random number extensions -*- C++ -*-
3 // Copyright (C) 2012 Free Software Foundation, Inc.
4 //
5 // This file is part of the GNU ISO C++ Library. This library is free
6 // software; you can redistribute it and/or modify it under the
7 // terms of the GNU General Public License as published by the
8 // Free Software Foundation; either version 3, or (at your option)
9 // any later version.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU General Public License for more details.
16 // Under Section 7 of GPL version 3, you are granted additional
17 // permissions described in the GCC Runtime Library Exception, version
18 // 3.1, as published by the Free Software Foundation.
20 // You should have received a copy of the GNU General Public License and
21 // a copy of the GCC Runtime Library Exception along with this program;
22 // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
23 // <http://www.gnu.org/licenses/>.
25 /** @file ext/random
26 * This file is a GNU extension to the Standard C++ Library.
27 */
29 #ifndef _EXT_RANDOM
30 #define _EXT_RANDOM 1
32 #pragma GCC system_header
34 #include <random>
35 #include <array>
36 #ifdef __SSE2__
37 # include <x86intrin.h>
38 #endif
41 namespace __gnu_cxx _GLIBCXX_VISIBILITY(default)
42 {
43 _GLIBCXX_BEGIN_NAMESPACE_VERSION
45 /* Mersenne twister implementation optimized for vector operations.
46 *
47 * Reference: http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/
48 */
49 template<typename _UIntType, size_t __m,
50 size_t __pos1, size_t __sl1, size_t __sl2,
51 size_t __sr1, size_t __sr2,
52 uint32_t __msk1, uint32_t __msk2,
53 uint32_t __msk3, uint32_t __msk4,
54 uint32_t __parity1, uint32_t __parity2,
55 uint32_t __parity3, uint32_t __parity4>
56 class simd_fast_mersenne_twister_engine
57 {
58 static_assert(std::is_unsigned<_UIntType>::value, "template argument "
59 "substituting _UIntType not an unsigned integral type");
60 static_assert(__sr1 < 32, "first right shift too large");
61 static_assert(__sr2 < 16, "second right shift too large");
62 static_assert(__sl1 < 32, "first left shift too large");
63 static_assert(__sl2 < 16, "second left shift too large");
65 public:
66 typedef _UIntType result_type;
68 private:
69 static constexpr size_t m_w = sizeof(result_type) * 8;
70 static constexpr size_t _M_nstate = __m / 128 + 1;
71 static constexpr size_t _M_nstate32 = _M_nstate * 4;
73 static_assert(std::is_unsigned<_UIntType>::value, "template argument "
74 "substituting _UIntType not an unsigned integral type");
75 static_assert(__pos1 < _M_nstate, "POS1 not smaller than state size");
76 static_assert(16 % sizeof(_UIntType) == 0,
77 "UIntType size must divide 16");
79 public:
80 static constexpr size_t state_size = _M_nstate * (16
81 / sizeof(result_type));
82 static constexpr result_type default_seed = 5489u;
84 // constructors and member function
85 explicit
86 simd_fast_mersenne_twister_engine(result_type __sd = default_seed)
87 { seed(__sd); }
89 template<typename _Sseq, typename = typename
90 std::enable_if<!std::is_same<_Sseq,
91 simd_fast_mersenne_twister_engine>::value>
92 ::type>
93 explicit
94 simd_fast_mersenne_twister_engine(_Sseq& __q)
95 { seed(__q); }
97 void
98 seed(result_type __sd = default_seed);
100 template<typename _Sseq>
101 typename std::enable_if<std::is_class<_Sseq>::value>::type
102 seed(_Sseq& __q);
104 static constexpr result_type
105 min()
106 { return 0; };
108 static constexpr result_type
109 max()
110 { return std::numeric_limits<result_type>::max(); }
112 void
113 discard(unsigned long long __z);
115 result_type
116 operator()()
117 {
118 if (__builtin_expect(_M_pos >= state_size, 0))
119 _M_gen_rand();
121 return _M_stateT[_M_pos++];
122 }
124 template<typename _UIntType_2, size_t __m_2,
125 size_t __pos1_2, size_t __sl1_2, size_t __sl2_2,
126 size_t __sr1_2, size_t __sr2_2,
127 uint32_t __msk1_2, uint32_t __msk2_2,
128 uint32_t __msk3_2, uint32_t __msk4_2,
129 uint32_t __parity1_2, uint32_t __parity2_2,
130 uint32_t __parity3_2, uint32_t __parity4_2>
131 friend bool
132 operator==(const simd_fast_mersenne_twister_engine<_UIntType_2,
133 __m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2,
134 __msk1_2, __msk2_2, __msk3_2, __msk4_2,
135 __parity1_2, __parity2_2, __parity3_2, __parity4_2>& __lhs,
136 const simd_fast_mersenne_twister_engine<_UIntType_2,
137 __m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2,
138 __msk1_2, __msk2_2, __msk3_2, __msk4_2,
139 __parity1_2, __parity2_2, __parity3_2, __parity4_2>& __rhs);
141 template<typename _UIntType_2, size_t __m_2,
142 size_t __pos1_2, size_t __sl1_2, size_t __sl2_2,
143 size_t __sr1_2, size_t __sr2_2,
144 uint32_t __msk1_2, uint32_t __msk2_2,
145 uint32_t __msk3_2, uint32_t __msk4_2,
146 uint32_t __parity1_2, uint32_t __parity2_2,
147 uint32_t __parity3_2, uint32_t __parity4_2,
148 typename _CharT, typename _Traits>
149 friend std::basic_ostream<_CharT, _Traits>&
150 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
151 const __gnu_cxx::simd_fast_mersenne_twister_engine
152 <_UIntType_2,
153 __m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2,
154 __msk1_2, __msk2_2, __msk3_2, __msk4_2,
155 __parity1_2, __parity2_2, __parity3_2, __parity4_2>& __x);
157 template<typename _UIntType_2, size_t __m_2,
158 size_t __pos1_2, size_t __sl1_2, size_t __sl2_2,
159 size_t __sr1_2, size_t __sr2_2,
160 uint32_t __msk1_2, uint32_t __msk2_2,
161 uint32_t __msk3_2, uint32_t __msk4_2,
162 uint32_t __parity1_2, uint32_t __parity2_2,
163 uint32_t __parity3_2, uint32_t __parity4_2,
164 typename _CharT, typename _Traits>
165 friend std::basic_istream<_CharT, _Traits>&
166 operator>>(std::basic_istream<_CharT, _Traits>& __is,
167 __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType_2,
168 __m_2, __pos1_2, __sl1_2, __sl2_2, __sr1_2, __sr2_2,
169 __msk1_2, __msk2_2, __msk3_2, __msk4_2,
170 __parity1_2, __parity2_2, __parity3_2, __parity4_2>& __x);
172 private:
173 union
174 {
175 #ifdef __SSE2__
176 __m128i _M_state[_M_nstate];
177 #endif
178 uint32_t _M_state32[_M_nstate32];
179 result_type _M_stateT[state_size];
180 } __attribute__ ((__aligned__ (16)));
181 size_t _M_pos;
183 void _M_gen_rand(void);
184 void _M_period_certification();
185 };
188 template<typename _UIntType, size_t __m,
189 size_t __pos1, size_t __sl1, size_t __sl2,
190 size_t __sr1, size_t __sr2,
191 uint32_t __msk1, uint32_t __msk2,
192 uint32_t __msk3, uint32_t __msk4,
193 uint32_t __parity1, uint32_t __parity2,
194 uint32_t __parity3, uint32_t __parity4>
195 inline bool
196 operator!=(const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
197 __m, __pos1, __sl1, __sl2, __sr1, __sr2, __msk1, __msk2, __msk3,
198 __msk4, __parity1, __parity2, __parity3, __parity4>& __lhs,
199 const __gnu_cxx::simd_fast_mersenne_twister_engine<_UIntType,
200 __m, __pos1, __sl1, __sl2, __sr1, __sr2, __msk1, __msk2, __msk3,
201 __msk4, __parity1, __parity2, __parity3, __parity4>& __rhs)
202 { return !(__lhs == __rhs); }
205 /* Definitions for the SIMD-oriented Fast Mersenne Twister as defined
206 * in the C implementation by Daito and Matsumoto, as both a 32-bit
207 * and 64-bit version.
208 */
209 typedef simd_fast_mersenne_twister_engine<uint32_t, 607, 2,
210 15, 3, 13, 3,
211 0xfdff37ffU, 0xef7f3f7dU,
212 0xff777b7dU, 0x7ff7fb2fU,
213 0x00000001U, 0x00000000U,
214 0x00000000U, 0x5986f054U>
215 sfmt607;
217 typedef simd_fast_mersenne_twister_engine<uint64_t, 607, 2,
218 15, 3, 13, 3,
219 0xfdff37ffU, 0xef7f3f7dU,
220 0xff777b7dU, 0x7ff7fb2fU,
221 0x00000001U, 0x00000000U,
222 0x00000000U, 0x5986f054U>
223 sfmt607_64;
226 typedef simd_fast_mersenne_twister_engine<uint32_t, 1279, 7,
227 14, 3, 5, 1,
228 0xf7fefffdU, 0x7fefcfffU,
229 0xaff3ef3fU, 0xb5ffff7fU,
230 0x00000001U, 0x00000000U,
231 0x00000000U, 0x20000000U>
232 sfmt1279;
234 typedef simd_fast_mersenne_twister_engine<uint64_t, 1279, 7,
235 14, 3, 5, 1,
236 0xf7fefffdU, 0x7fefcfffU,
237 0xaff3ef3fU, 0xb5ffff7fU,
238 0x00000001U, 0x00000000U,
239 0x00000000U, 0x20000000U>
240 sfmt1279_64;
243 typedef simd_fast_mersenne_twister_engine<uint32_t, 2281, 12,
244 19, 1, 5, 1,
245 0xbff7ffbfU, 0xfdfffffeU,
246 0xf7ffef7fU, 0xf2f7cbbfU,
247 0x00000001U, 0x00000000U,
248 0x00000000U, 0x41dfa600U>
249 sfmt2281;
251 typedef simd_fast_mersenne_twister_engine<uint64_t, 2281, 12,
252 19, 1, 5, 1,
253 0xbff7ffbfU, 0xfdfffffeU,
254 0xf7ffef7fU, 0xf2f7cbbfU,
255 0x00000001U, 0x00000000U,
256 0x00000000U, 0x41dfa600U>
257 sfmt2281_64;
260 typedef simd_fast_mersenne_twister_engine<uint32_t, 4253, 17,
261 20, 1, 7, 1,
262 0x9f7bffffU, 0x9fffff5fU,
263 0x3efffffbU, 0xfffff7bbU,
264 0xa8000001U, 0xaf5390a3U,
265 0xb740b3f8U, 0x6c11486dU>
266 sfmt4253;
268 typedef simd_fast_mersenne_twister_engine<uint64_t, 4253, 17,
269 20, 1, 7, 1,
270 0x9f7bffffU, 0x9fffff5fU,
271 0x3efffffbU, 0xfffff7bbU,
272 0xa8000001U, 0xaf5390a3U,
273 0xb740b3f8U, 0x6c11486dU>
274 sfmt4253_64;
277 typedef simd_fast_mersenne_twister_engine<uint32_t, 11213, 68,
278 14, 3, 7, 3,
279 0xeffff7fbU, 0xffffffefU,
280 0xdfdfbfffU, 0x7fffdbfdU,
281 0x00000001U, 0x00000000U,
282 0xe8148000U, 0xd0c7afa3U>
283 sfmt11213;
285 typedef simd_fast_mersenne_twister_engine<uint64_t, 11213, 68,
286 14, 3, 7, 3,
287 0xeffff7fbU, 0xffffffefU,
288 0xdfdfbfffU, 0x7fffdbfdU,
289 0x00000001U, 0x00000000U,
290 0xe8148000U, 0xd0c7afa3U>
291 sfmt11213_64;
294 typedef simd_fast_mersenne_twister_engine<uint32_t, 19937, 122,
295 18, 1, 11, 1,
296 0xdfffffefU, 0xddfecb7fU,
297 0xbffaffffU, 0xbffffff6U,
298 0x00000001U, 0x00000000U,
299 0x00000000U, 0x13c9e684U>
300 sfmt19937;
302 typedef simd_fast_mersenne_twister_engine<uint64_t, 19937, 122,
303 18, 1, 11, 1,
304 0xdfffffefU, 0xddfecb7fU,
305 0xbffaffffU, 0xbffffff6U,
306 0x00000001U, 0x00000000U,
307 0x00000000U, 0x13c9e684U>
308 sfmt19937_64;
311 typedef simd_fast_mersenne_twister_engine<uint32_t, 44497, 330,
312 5, 3, 9, 3,
313 0xeffffffbU, 0xdfbebfffU,
314 0xbfbf7befU, 0x9ffd7bffU,
315 0x00000001U, 0x00000000U,
316 0xa3ac4000U, 0xecc1327aU>
317 sfmt44497;
319 typedef simd_fast_mersenne_twister_engine<uint64_t, 44497, 330,
320 5, 3, 9, 3,
321 0xeffffffbU, 0xdfbebfffU,
322 0xbfbf7befU, 0x9ffd7bffU,
323 0x00000001U, 0x00000000U,
324 0xa3ac4000U, 0xecc1327aU>
325 sfmt44497_64;
328 typedef simd_fast_mersenne_twister_engine<uint32_t, 86243, 366,
329 6, 7, 19, 1,
330 0xfdbffbffU, 0xbff7ff3fU,
331 0xfd77efffU, 0xbf9ff3ffU,
332 0x00000001U, 0x00000000U,
333 0x00000000U, 0xe9528d85U>
334 sfmt86243;
336 typedef simd_fast_mersenne_twister_engine<uint64_t, 86243, 366,
337 6, 7, 19, 1,
338 0xfdbffbffU, 0xbff7ff3fU,
339 0xfd77efffU, 0xbf9ff3ffU,
340 0x00000001U, 0x00000000U,
341 0x00000000U, 0xe9528d85U>
342 sfmt86243_64;
345 typedef simd_fast_mersenne_twister_engine<uint32_t, 132049, 110,
346 19, 1, 21, 1,
347 0xffffbb5fU, 0xfb6ebf95U,
348 0xfffefffaU, 0xcff77fffU,
349 0x00000001U, 0x00000000U,
350 0xcb520000U, 0xc7e91c7dU>
351 sfmt132049;
353 typedef simd_fast_mersenne_twister_engine<uint64_t, 132049, 110,
354 19, 1, 21, 1,
355 0xffffbb5fU, 0xfb6ebf95U,
356 0xfffefffaU, 0xcff77fffU,
357 0x00000001U, 0x00000000U,
358 0xcb520000U, 0xc7e91c7dU>
359 sfmt132049_64;
362 typedef simd_fast_mersenne_twister_engine<uint32_t, 216091, 627,
363 11, 3, 10, 1,
364 0xbff7bff7U, 0xbfffffffU,
365 0xbffffa7fU, 0xffddfbfbU,
366 0xf8000001U, 0x89e80709U,
367 0x3bd2b64bU, 0x0c64b1e4U>
368 sfmt216091;
370 typedef simd_fast_mersenne_twister_engine<uint64_t, 216091, 627,
371 11, 3, 10, 1,
372 0xbff7bff7U, 0xbfffffffU,
373 0xbffffa7fU, 0xffddfbfbU,
374 0xf8000001U, 0x89e80709U,
375 0x3bd2b64bU, 0x0c64b1e4U>
376 sfmt216091_64;
379 /**
380 * @brief A beta continuous distribution for random numbers.
381 *
382 * The formula for the beta probability density function is:
383 * @f[
384 * p(x|\alpha,\beta) = \frac{1}{B(\alpha,\beta)}
385 * x^{\alpha - 1} (1 - x)^{\beta - 1}
386 * @f]
387 */
388 template<typename _RealType = double>
389 class beta_distribution
390 {
391 static_assert(std::is_floating_point<_RealType>::value,
392 "template argument not a floating point type");
394 public:
395 /** The type of the range of the distribution. */
396 typedef _RealType result_type;
397 /** Parameter type. */
398 struct param_type
399 {
400 typedef beta_distribution<_RealType> distribution_type;
401 friend class beta_distribution<_RealType>;
403 explicit
404 param_type(_RealType __alpha_val = _RealType(1),
405 _RealType __beta_val = _RealType(1))
406 : _M_alpha(__alpha_val), _M_beta(__beta_val)
407 {
408 _GLIBCXX_DEBUG_ASSERT(_M_alpha > _RealType(0));
409 _GLIBCXX_DEBUG_ASSERT(_M_beta > _RealType(0));
410 }
412 _RealType
413 alpha() const
414 { return _M_alpha; }
416 _RealType
417 beta() const
418 { return _M_beta; }
420 friend bool
421 operator==(const param_type& __p1, const param_type& __p2)
422 { return (__p1._M_alpha == __p2._M_alpha
423 && __p1._M_beta == __p2._M_beta); }
425 private:
426 void
427 _M_initialize();
429 _RealType _M_alpha;
430 _RealType _M_beta;
431 };
433 public:
434 /**
435 * @brief Constructs a beta distribution with parameters
436 * @f$\alpha@f$ and @f$\beta@f$.
437 */
438 explicit
439 beta_distribution(_RealType __alpha_val = _RealType(1),
440 _RealType __beta_val = _RealType(1))
441 : _M_param(__alpha_val, __beta_val)
442 { }
444 explicit
445 beta_distribution(const param_type& __p)
446 : _M_param(__p)
447 { }
449 /**
450 * @brief Resets the distribution state.
451 */
452 void
453 reset()
454 { }
456 /**
457 * @brief Returns the @f$\alpha@f$ of the distribution.
458 */
459 _RealType
460 alpha() const
461 { return _M_param.alpha(); }
463 /**
464 * @brief Returns the @f$\beta@f$ of the distribution.
465 */
466 _RealType
467 beta() const
468 { return _M_param.beta(); }
470 /**
471 * @brief Returns the parameter set of the distribution.
472 */
473 param_type
474 param() const
475 { return _M_param; }
477 /**
478 * @brief Sets the parameter set of the distribution.
479 * @param __param The new parameter set of the distribution.
480 */
481 void
482 param(const param_type& __param)
483 { _M_param = __param; }
485 /**
486 * @brief Returns the greatest lower bound value of the distribution.
487 */
488 result_type
489 min() const
490 { return result_type(0); }
492 /**
493 * @brief Returns the least upper bound value of the distribution.
494 */
495 result_type
496 max() const
497 { return result_type(1); }
499 /**
500 * @brief Generating functions.
501 */
502 template<typename _UniformRandomNumberGenerator>
503 result_type
504 operator()(_UniformRandomNumberGenerator& __urng)
505 { return this->operator()(__urng, this->param()); }
507 template<typename _UniformRandomNumberGenerator>
508 result_type
509 operator()(_UniformRandomNumberGenerator& __urng,
510 const param_type& __p);
512 template<typename _ForwardIterator,
513 typename _UniformRandomNumberGenerator>
514 void
515 __generate(_ForwardIterator __f, _ForwardIterator __t,
516 _UniformRandomNumberGenerator& __urng)
517 { this->__generate(__f, __t, __urng, this->param()); }
519 template<typename _ForwardIterator,
520 typename _UniformRandomNumberGenerator>
521 void
522 __generate(_ForwardIterator __f, _ForwardIterator __t,
523 _UniformRandomNumberGenerator& __urng,
524 const param_type& __p)
525 { this->__generate_impl(__f, __t, __urng, __p); }
527 template<typename _UniformRandomNumberGenerator>
528 void
529 __generate(result_type* __f, result_type* __t,
530 _UniformRandomNumberGenerator& __urng,
531 const param_type& __p)
532 { this->__generate_impl(__f, __t, __urng, __p); }
534 /**
535 * @brief Inserts a %beta_distribution random number distribution
536 * @p __x into the output stream @p __os.
537 *
538 * @param __os An output stream.
539 * @param __x A %beta_distribution random number distribution.
540 *
541 * @returns The output stream with the state of @p __x inserted or in
542 * an error state.
543 */
544 template<typename _RealType1, typename _CharT, typename _Traits>
545 friend std::basic_ostream<_CharT, _Traits>&
546 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
547 const __gnu_cxx::beta_distribution<_RealType1>& __x);
549 /**
550 * @brief Extracts a %beta_distribution random number distribution
551 * @p __x from the input stream @p __is.
552 *
553 * @param __is An input stream.
554 * @param __x A %beta_distribution random number generator engine.
555 *
556 * @returns The input stream with @p __x extracted or in an error state.
557 */
558 template<typename _RealType1, typename _CharT, typename _Traits>
559 friend std::basic_istream<_CharT, _Traits>&
560 operator>>(std::basic_istream<_CharT, _Traits>& __is,
561 __gnu_cxx::beta_distribution<_RealType1>& __x);
563 private:
564 template<typename _ForwardIterator,
565 typename _UniformRandomNumberGenerator>
566 void
567 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
568 _UniformRandomNumberGenerator& __urng,
569 const param_type& __p);
571 param_type _M_param;
572 };
574 /**
575 * @brief Return true if two beta distributions have the same
576 * parameters and the sequences that would be generated
577 * are equal.
578 */
579 template<typename _RealType>
580 inline bool
581 operator==(const __gnu_cxx::beta_distribution<_RealType>& __d1,
582 const __gnu_cxx::beta_distribution<_RealType>& __d2)
583 { return __d1.param() == __d2.param(); }
585 /**
586 * @brief Return true if two beta distributions are different.
587 */
588 template<typename _RealType>
589 inline bool
590 operator!=(const __gnu_cxx::beta_distribution<_RealType>& __d1,
591 const __gnu_cxx::beta_distribution<_RealType>& __d2)
592 { return !(__d1 == __d2); }
595 /**
596 * @brief A multi-variate normal continuous distribution for random numbers.
597 *
598 * The formula for the normal probability density function is
599 * @f[
600 * p(\overrightarrow{x}|\overrightarrow{\mu },\Sigma) =
601 * \frac{1}{\sqrt{(2\pi )^k\det(\Sigma))}}
602 * e^{-\frac{1}{2}(\overrightarrow{x}-\overrightarrow{\mu})^\text{T}
603 * \Sigma ^{-1}(\overrightarrow{x}-\overrightarrow{\mu})}
604 * @f]
605 *
606 * where @f$\overrightarrow{x}@f$ and @f$\overrightarrow{\mu}@f$ are
607 * vectors of dimension @f$k@f$ and @f$\Sigma@f$ is the covariance
608 * matrix (which must be positive-definite).
609 */
610 template<std::size_t _Dimen, typename _RealType = double>
611 class normal_mv_distribution
612 {
613 static_assert(std::is_floating_point<_RealType>::value,
614 "template argument not a floating point type");
615 static_assert(_Dimen != 0, "dimension is zero");
617 public:
618 /** The type of the range of the distribution. */
619 typedef std::array<_RealType, _Dimen> result_type;
620 /** Parameter type. */
621 class param_type
622 {
623 static constexpr size_t _M_t_size = _Dimen * (_Dimen + 1) / 2;
625 public:
626 typedef normal_mv_distribution<_Dimen, _RealType> distribution_type;
627 friend class normal_mv_distribution<_Dimen, _RealType>;
629 param_type()
630 {
631 std::fill(_M_mean.begin(), _M_mean.end(), _RealType(0));
632 auto __it = _M_t.begin();
633 for (size_t __i = 0; __i < _Dimen; ++__i)
634 {
635 std::fill_n(__it, __i, _RealType(0));
636 __it += __i;
637 *__it++ = _RealType(1);
638 }
639 }
641 template<typename _ForwardIterator1, typename _ForwardIterator2>
642 param_type(_ForwardIterator1 __meanbegin,
643 _ForwardIterator1 __meanend,
644 _ForwardIterator2 __varcovbegin,
645 _ForwardIterator2 __varcovend)
646 {
647 __glibcxx_function_requires(_ForwardIteratorConcept<
648 _ForwardIterator1>)
649 __glibcxx_function_requires(_ForwardIteratorConcept<
650 _ForwardIterator2>)
651 _GLIBCXX_DEBUG_ASSERT(std::distance(__meanbegin, __meanend)
652 <= _Dimen);
653 const auto __dist = std::distance(__varcovbegin, __varcovend);
654 _GLIBCXX_DEBUG_ASSERT(__dist == _Dimen * _Dimen
655 || __dist == _Dimen * (_Dimen + 1) / 2
656 || __dist == _Dimen);
658 if (__dist == _Dimen * _Dimen)
659 _M_init_full(__meanbegin, __meanend, __varcovbegin, __varcovend);
660 else if (__dist == _Dimen * (_Dimen + 1) / 2)
661 _M_init_lower(__meanbegin, __meanend, __varcovbegin, __varcovend);
662 else
663 _M_init_diagonal(__meanbegin, __meanend,
664 __varcovbegin, __varcovend);
665 }
667 param_type(std::initializer_list<_RealType> __mean,
668 std::initializer_list<_RealType> __varcov)
669 {
670 _GLIBCXX_DEBUG_ASSERT(__mean.size() <= _Dimen);
671 _GLIBCXX_DEBUG_ASSERT(__varcov.size() == _Dimen * _Dimen
672 || __varcov.size() == _Dimen * (_Dimen + 1) / 2
673 || __varcov.size() == _Dimen);
675 if (__varcov.size() == _Dimen * _Dimen)
676 _M_init_full(__mean.begin(), __mean.end(),
677 __varcov.begin(), __varcov.end());
678 else if (__varcov.size() == _Dimen * (_Dimen + 1) / 2)
679 _M_init_lower(__mean.begin(), __mean.end(),
680 __varcov.begin(), __varcov.end());
681 else
682 _M_init_diagonal(__mean.begin(), __mean.end(),
683 __varcov.begin(), __varcov.end());
684 }
686 std::array<_RealType, _Dimen>
687 mean() const
688 { return _M_mean; }
690 std::array<_RealType, _M_t_size>
691 varcov() const
692 { return _M_t; }
694 friend bool
695 operator==(const param_type& __p1, const param_type& __p2)
696 { return __p1._M_mean == __p2._M_mean && __p1._M_t == __p2._M_t; }
698 private:
699 template <typename _InputIterator1, typename _InputIterator2>
700 void _M_init_full(_InputIterator1 __meanbegin,
701 _InputIterator1 __meanend,
702 _InputIterator2 __varcovbegin,
703 _InputIterator2 __varcovend);
704 template <typename _InputIterator1, typename _InputIterator2>
705 void _M_init_lower(_InputIterator1 __meanbegin,
706 _InputIterator1 __meanend,
707 _InputIterator2 __varcovbegin,
708 _InputIterator2 __varcovend);
709 template <typename _InputIterator1, typename _InputIterator2>
710 void _M_init_diagonal(_InputIterator1 __meanbegin,
711 _InputIterator1 __meanend,
712 _InputIterator2 __varbegin,
713 _InputIterator2 __varend);
715 std::array<_RealType, _Dimen> _M_mean;
716 std::array<_RealType, _M_t_size> _M_t;
717 };
719 public:
720 normal_mv_distribution()
721 : _M_param(), _M_nd()
722 { }
724 template<typename _ForwardIterator1, typename _ForwardIterator2>
725 normal_mv_distribution(_ForwardIterator1 __meanbegin,
726 _ForwardIterator1 __meanend,
727 _ForwardIterator2 __varcovbegin,
728 _ForwardIterator2 __varcovend)
729 : _M_param(__meanbegin, __meanend, __varcovbegin, __varcovend),
730 _M_nd()
731 { }
733 normal_mv_distribution(std::initializer_list<_RealType> __mean,
734 std::initializer_list<_RealType> __varcov)
735 : _M_param(__mean, __varcov), _M_nd()
736 { }
738 explicit
739 normal_mv_distribution(const param_type& __p)
740 : _M_param(__p), _M_nd()
741 { }
743 /**
744 * @brief Resets the distribution state.
745 */
746 void
747 reset()
748 { _M_nd.reset(); }
750 /**
751 * @brief Returns the mean of the distribution.
752 */
753 result_type
754 mean() const
755 { return _M_param.mean(); }
757 /**
758 * @brief Returns the compact form of the variance/covariance
759 * matrix of the distribution.
760 */
761 std::array<_RealType, _Dimen * (_Dimen + 1) / 2>
762 varcov() const
763 { return _M_param.varcov(); }
765 /**
766 * @brief Returns the parameter set of the distribution.
767 */
768 param_type
769 param() const
770 { return _M_param; }
772 /**
773 * @brief Sets the parameter set of the distribution.
774 * @param __param The new parameter set of the distribution.
775 */
776 void
777 param(const param_type& __param)
778 { _M_param = __param; }
780 /**
781 * @brief Returns the greatest lower bound value of the distribution.
782 */
783 result_type
784 min() const
785 { result_type __res;
786 __res.fill(std::numeric_limits<_RealType>::min());
787 return __res; }
789 /**
790 * @brief Returns the least upper bound value of the distribution.
791 */
792 result_type
793 max() const
794 { result_type __res;
795 __res.fill(std::numeric_limits<_RealType>::max());
796 return __res; }
798 /**
799 * @brief Generating functions.
800 */
801 template<typename _UniformRandomNumberGenerator>
802 result_type
803 operator()(_UniformRandomNumberGenerator& __urng)
804 { return this->operator()(__urng, this->param()); }
806 template<typename _UniformRandomNumberGenerator>
807 result_type
808 operator()(_UniformRandomNumberGenerator& __urng,
809 const param_type& __p);
811 template<typename _ForwardIterator,
812 typename _UniformRandomNumberGenerator>
813 void
814 __generate(_ForwardIterator __f, _ForwardIterator __t,
815 _UniformRandomNumberGenerator& __urng)
816 { return this->__generate_impl(__f, __t, __urng, this->param()); }
818 template<typename _ForwardIterator,
819 typename _UniformRandomNumberGenerator>
820 void
821 __generate(_ForwardIterator __f, _ForwardIterator __t,
822 _UniformRandomNumberGenerator& __urng,
823 const param_type& __p)
824 { return this->__generate_impl(__f, __t, __urng, __p); }
826 /**
827 * @brief Return true if two multi-variant normal distributions have
828 * the same parameters and the sequences that would
829 * be generated are equal.
830 */
831 template<size_t _Dimen1, typename _RealType1>
832 friend bool
833 operator==(const
834 __gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>&
835 __d1,
836 const
837 __gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>&
838 __d2);
840 /**
841 * @brief Inserts a %normal_mv_distribution random number distribution
842 * @p __x into the output stream @p __os.
843 *
844 * @param __os An output stream.
845 * @param __x A %normal_mv_distribution random number distribution.
846 *
847 * @returns The output stream with the state of @p __x inserted or in
848 * an error state.
849 */
850 template<size_t _Dimen1, typename _RealType1,
851 typename _CharT, typename _Traits>
852 friend std::basic_ostream<_CharT, _Traits>&
853 operator<<(std::basic_ostream<_CharT, _Traits>& __os,
854 const
855 __gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>&
856 __x);
858 /**
859 * @brief Extracts a %normal_mv_distribution random number distribution
860 * @p __x from the input stream @p __is.
861 *
862 * @param __is An input stream.
863 * @param __x A %normal_mv_distribution random number generator engine.
864 *
865 * @returns The input stream with @p __x extracted or in an error
866 * state.
867 */
868 template<size_t _Dimen1, typename _RealType1,
869 typename _CharT, typename _Traits>
870 friend std::basic_istream<_CharT, _Traits>&
871 operator>>(std::basic_istream<_CharT, _Traits>& __is,
872 __gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>&
873 __x);
875 private:
876 template<typename _ForwardIterator,
877 typename _UniformRandomNumberGenerator>
878 void
879 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
880 _UniformRandomNumberGenerator& __urng,
881 const param_type& __p);
883 param_type _M_param;
884 std::normal_distribution<_RealType> _M_nd;
885 };
887 /**
888 * @brief Return true if two multi-variate normal distributions are
889 * different.
890 */
891 template<size_t _Dimen, typename _RealType>
892 inline bool
893 operator!=(const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
894 __d1,
895 const __gnu_cxx::normal_mv_distribution<_Dimen, _RealType>&
896 __d2)
897 { return !(__d1 == __d2); }
900 /**
901 * @brief A Rice continuous distribution for random numbers.
902 *
903 * The formula for the Rice probability density function is
904 * @f[
905 * p(x|\nu,\sigma) = \frac{x}{\sigma^2}
906 * \exp\left(-\frac{x^2+\nu^2}{2\sigma^2}\right)
907 * I_0\left(\frac{x \nu}{\sigma^2}\right)
908 * @f]
909 * where @f$I_0(z)@f$ is the modified Bessel function of the first kind
910 * of order 0 and @f$\nu >= 0@f$ and @f$\sigma > 0@f$.
911 *
912 * <table border=1 cellpadding=10 cellspacing=0>
913 * <caption align=top>Distribution Statistics</caption>
914 * <tr><td>Mean</td><td>@f$\sqrt{\pi/2}L_{1/2}(-\nu^2/2\sigma^2)@f$</td></tr>
915 * <tr><td>Variance</td><td>@f$2\sigma^2 + \nu^2
916 * + (\pi\sigma^2/2)L^2_{1/2}(-\nu^2/2\sigma^2)@f$</td></tr>
917 * <tr><td>Range</td><td>@f$[0, \infty)@f$</td></tr>
918 * </table>
919 * where @f$L_{1/2}(x)@f$ is the Laguerre polynomial of order 1/2.
920 */
921 template<typename _RealType = double>
922 class
923 rice_distribution
924 {
925 static_assert(std::is_floating_point<_RealType>::value,
926 "template argument not a floating point type");
927 public:
928 /** The type of the range of the distribution. */
929 typedef _RealType result_type;
930 /** Parameter type. */
931 struct param_type
932 {
933 typedef rice_distribution<result_type> distribution_type;
935 param_type(result_type __nu_val = result_type(0),
936 result_type __sigma_val = result_type(1))
937 : _M_nu(__nu_val), _M_sigma(__sigma_val)
938 {
939 _GLIBCXX_DEBUG_ASSERT(_M_nu >= result_type(0));
940 _GLIBCXX_DEBUG_ASSERT(_M_sigma > result_type(0));
941 }
943 result_type
944 nu() const
945 { return _M_nu; }
947 result_type
948 sigma() const
949 { return _M_sigma; }
951 friend bool
952 operator==(const param_type& __p1, const param_type& __p2)
953 { return __p1._M_nu == __p2._M_nu
954 && __p1._M_sigma == __p2._M_sigma; }
956 private:
957 void _M_initialize();
959 result_type _M_nu;
960 result_type _M_sigma;
961 };
963 /**
964 * @brief Constructors.
965 */
966 explicit
967 rice_distribution(result_type __nu_val = result_type(0),
968 result_type __sigma_val = result_type(1))
969 : _M_param(__nu_val, __sigma_val),
970 _M_ndx(__nu_val, __sigma_val),
971 _M_ndy(result_type(0), __sigma_val)
972 { }
974 explicit
975 rice_distribution(const param_type& __p)
976 : _M_param(__p),
977 _M_ndx(__p.nu(), __p.sigma()),
978 _M_ndy(result_type(0), __p.sigma())
979 { }
981 /**
982 * @brief Resets the distribution state.
983 */
984 void
985 reset()
986 {
987 _M_ndx.reset();
988 _M_ndy.reset();
989 }
991 /**
992 * @brief Return the parameters of the distribution.
993 */
994 result_type
995 nu() const
996 { return _M_param.nu(); }
998 result_type
999 sigma() const
1000 { return _M_param.sigma(); }
1002 /**
1003 * @brief Returns the parameter set of the distribution.
1004 */
1005 param_type
1006 param() const
1007 { return _M_param; }
1009 /**
1010 * @brief Sets the parameter set of the distribution.
1011 * @param __param The new parameter set of the distribution.
1012 */
1013 void
1014 param(const param_type& __param)
1015 { _M_param = __param; }
1017 /**
1018 * @brief Returns the greatest lower bound value of the distribution.
1019 */
1020 result_type
1021 min() const
1022 { return result_type(0); }
1024 /**
1025 * @brief Returns the least upper bound value of the distribution.
1026 */
1027 result_type
1028 max() const
1029 { return std::numeric_limits<result_type>::max(); }
1031 /**
1032 * @brief Generating functions.
1033 */
1034 template<typename _UniformRandomNumberGenerator>
1035 result_type
1036 operator()(_UniformRandomNumberGenerator& __urng)
1037 {
1038 result_type __x = this->_M_ndx(__urng);
1039 result_type __y = this->_M_ndy(__urng);
1040 return std::hypot(__x, __y);
1041 }
1043 template<typename _UniformRandomNumberGenerator>
1044 result_type
1045 operator()(_UniformRandomNumberGenerator& __urng,
1046 const param_type& __p)
1047 {
1048 typename std::normal_distribution<result_type>::param_type
1049 __px(__p.nu(), __p.sigma()), __py(result_type(0), __p.sigma());
1050 result_type __x = this->_M_ndx(__px, __urng);
1051 result_type __y = this->_M_ndy(__py, __urng);
1052 return std::hypot(__x, __y);
1053 }
1055 template<typename _ForwardIterator,
1056 typename _UniformRandomNumberGenerator>
1057 void
1058 __generate(_ForwardIterator __f, _ForwardIterator __t,
1059 _UniformRandomNumberGenerator& __urng)
1060 { this->__generate(__f, __t, __urng, this->param()); }
1062 template<typename _ForwardIterator,
1063 typename _UniformRandomNumberGenerator>
1064 void
1065 __generate(_ForwardIterator __f, _ForwardIterator __t,
1066 _UniformRandomNumberGenerator& __urng,
1067 const param_type& __p)
1068 { this->__generate_impl(__f, __t, __urng, __p); }
1070 template<typename _UniformRandomNumberGenerator>
1071 void
1072 __generate(result_type* __f, result_type* __t,
1073 _UniformRandomNumberGenerator& __urng,
1074 const param_type& __p)
1075 { this->__generate_impl(__f, __t, __urng, __p); }
1077 /**
1078 * @brief Return true if two Rice distributions have
1079 * the same parameters and the sequences that would
1080 * be generated are equal.
1081 */
1082 friend bool
1083 operator==(const rice_distribution& __d1,
1084 const rice_distribution& __d2)
1085 { return (__d1.param() == __d2.param()
1086 && __d1._M_ndx == __d2._M_ndx
1087 && __d1._M_ndy == __d2._M_ndy); }
1089 /**
1090 * @brief Inserts a %rice_distribution random number distribution
1091 * @p __x into the output stream @p __os.
1092 *
1093 * @param __os An output stream.
1094 * @param __x A %rice_distribution random number distribution.
1095 *
1096 * @returns The output stream with the state of @p __x inserted or in
1097 * an error state.
1098 */
1099 template<typename _RealType1, typename _CharT, typename _Traits>
1100 friend std::basic_ostream<_CharT, _Traits>&
1101 operator<<(std::basic_ostream<_CharT, _Traits>&,
1102 const rice_distribution<_RealType1>&);
1104 /**
1105 * @brief Extracts a %rice_distribution random number distribution
1106 * @p __x from the input stream @p __is.
1107 *
1108 * @param __is An input stream.
1109 * @param __x A %rice_distribution random number
1110 * generator engine.
1111 *
1112 * @returns The input stream with @p __x extracted or in an error state.
1113 */
1114 template<typename _RealType1, typename _CharT, typename _Traits>
1115 friend std::basic_istream<_CharT, _Traits>&
1116 operator>>(std::basic_istream<_CharT, _Traits>&,
1117 rice_distribution<_RealType1>&);
1119 private:
1120 template<typename _ForwardIterator,
1121 typename _UniformRandomNumberGenerator>
1122 void
1123 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1124 _UniformRandomNumberGenerator& __urng,
1125 const param_type& __p);
1127 param_type _M_param;
1129 std::normal_distribution<result_type> _M_ndx;
1130 std::normal_distribution<result_type> _M_ndy;
1131 };
1133 /**
1134 * @brief Return true if two Rice distributions are not equal.
1135 */
1136 template<typename _RealType1>
1137 inline bool
1138 operator!=(const rice_distribution<_RealType1>& __d1,
1139 const rice_distribution<_RealType1>& __d2)
1140 { return !(__d1 == __d2); }
1143 /**
1144 * @brief A Nakagami continuous distribution for random numbers.
1145 *
1146 * The formula for the Nakagami probability density function is
1147 * @f[
1148 * p(x|\mu,\omega) = \frac{2\mu^\mu}{\Gamma(\mu)\omega^\mu}
1149 * x^{2\mu-1}e^{-\mu x / \omega}
1150 * @f]
1151 * where @f$\Gamma(z)@f$ is the gamma function and @f$\mu >= 0.5@f$
1152 * and @f$\omega > 0@f$.
1153 */
1154 template<typename _RealType = double>
1155 class
1156 nakagami_distribution
1157 {
1158 static_assert(std::is_floating_point<_RealType>::value,
1159 "template argument not a floating point type");
1161 public:
1162 /** The type of the range of the distribution. */
1163 typedef _RealType result_type;
1164 /** Parameter type. */
1165 struct param_type
1166 {
1167 typedef nakagami_distribution<result_type> distribution_type;
1169 param_type(result_type __mu_val = result_type(1),
1170 result_type __omega_val = result_type(1))
1171 : _M_mu(__mu_val), _M_omega(__omega_val)
1172 {
1173 _GLIBCXX_DEBUG_ASSERT(_M_mu >= result_type(0.5L));
1174 _GLIBCXX_DEBUG_ASSERT(_M_omega > result_type(0));
1175 }
1177 result_type
1178 mu() const
1179 { return _M_mu; }
1181 result_type
1182 omega() const
1183 { return _M_omega; }
1185 friend bool
1186 operator==(const param_type& __p1, const param_type& __p2)
1187 { return __p1._M_mu == __p2._M_mu
1188 && __p1._M_omega == __p2._M_omega; }
1190 private:
1191 void _M_initialize();
1193 result_type _M_mu;
1194 result_type _M_omega;
1195 };
1197 /**
1198 * @brief Constructors.
1199 */
1200 explicit
1201 nakagami_distribution(result_type __mu_val = result_type(1),
1202 result_type __omega_val = result_type(1))
1203 : _M_param(__mu_val, __omega_val),
1204 _M_gd(__mu_val, __omega_val / __mu_val)
1205 { }
1207 explicit
1208 nakagami_distribution(const param_type& __p)
1209 : _M_param(__p),
1210 _M_gd(__p.mu(), __p.omega() / __p.mu())
1211 { }
1213 /**
1214 * @brief Resets the distribution state.
1215 */
1216 void
1217 reset()
1218 { _M_gd.reset(); }
1220 /**
1221 * @brief Return the parameters of the distribution.
1222 */
1223 result_type
1224 mu() const
1225 { return _M_param.mu(); }
1227 result_type
1228 omega() const
1229 { return _M_param.omega(); }
1231 /**
1232 * @brief Returns the parameter set of the distribution.
1233 */
1234 param_type
1235 param() const
1236 { return _M_param; }
1238 /**
1239 * @brief Sets the parameter set of the distribution.
1240 * @param __param The new parameter set of the distribution.
1241 */
1242 void
1243 param(const param_type& __param)
1244 { _M_param = __param; }
1246 /**
1247 * @brief Returns the greatest lower bound value of the distribution.
1248 */
1249 result_type
1250 min() const
1251 { return result_type(0); }
1253 /**
1254 * @brief Returns the least upper bound value of the distribution.
1255 */
1256 result_type
1257 max() const
1258 { return std::numeric_limits<result_type>::max(); }
1260 /**
1261 * @brief Generating functions.
1262 */
1263 template<typename _UniformRandomNumberGenerator>
1264 result_type
1265 operator()(_UniformRandomNumberGenerator& __urng)
1266 { return std::sqrt(this->_M_gd(__urng)); }
1268 template<typename _UniformRandomNumberGenerator>
1269 result_type
1270 operator()(_UniformRandomNumberGenerator& __urng,
1271 const param_type& __p)
1272 {
1273 typename std::gamma_distribution<result_type>::param_type
1274 __pg(__p.mu(), __p.omega() / __p.mu());
1275 return std::sqrt(this->_M_gd(__pg, __urng));
1276 }
1278 template<typename _ForwardIterator,
1279 typename _UniformRandomNumberGenerator>
1280 void
1281 __generate(_ForwardIterator __f, _ForwardIterator __t,
1282 _UniformRandomNumberGenerator& __urng)
1283 { this->__generate(__f, __t, __urng, this->param()); }
1285 template<typename _ForwardIterator,
1286 typename _UniformRandomNumberGenerator>
1287 void
1288 __generate(_ForwardIterator __f, _ForwardIterator __t,
1289 _UniformRandomNumberGenerator& __urng,
1290 const param_type& __p)
1291 { this->__generate_impl(__f, __t, __urng, __p); }
1293 template<typename _UniformRandomNumberGenerator>
1294 void
1295 __generate(result_type* __f, result_type* __t,
1296 _UniformRandomNumberGenerator& __urng,
1297 const param_type& __p)
1298 { this->__generate_impl(__f, __t, __urng, __p); }
1300 /**
1301 * @brief Return true if two Nakagami distributions have
1302 * the same parameters and the sequences that would
1303 * be generated are equal.
1304 */
1305 friend bool
1306 operator==(const nakagami_distribution& __d1,
1307 const nakagami_distribution& __d2)
1308 { return (__d1.param() == __d2.param()
1309 && __d1._M_gd == __d2._M_gd); }
1311 /**
1312 * @brief Inserts a %nakagami_distribution random number distribution
1313 * @p __x into the output stream @p __os.
1314 *
1315 * @param __os An output stream.
1316 * @param __x A %nakagami_distribution random number distribution.
1317 *
1318 * @returns The output stream with the state of @p __x inserted or in
1319 * an error state.
1320 */
1321 template<typename _RealType1, typename _CharT, typename _Traits>
1322 friend std::basic_ostream<_CharT, _Traits>&
1323 operator<<(std::basic_ostream<_CharT, _Traits>&,
1324 const nakagami_distribution<_RealType1>&);
1326 /**
1327 * @brief Extracts a %nakagami_distribution random number distribution
1328 * @p __x from the input stream @p __is.
1329 *
1330 * @param __is An input stream.
1331 * @param __x A %nakagami_distribution random number
1332 * generator engine.
1333 *
1334 * @returns The input stream with @p __x extracted or in an error state.
1335 */
1336 template<typename _RealType1, typename _CharT, typename _Traits>
1337 friend std::basic_istream<_CharT, _Traits>&
1338 operator>>(std::basic_istream<_CharT, _Traits>&,
1339 nakagami_distribution<_RealType1>&);
1341 private:
1342 template<typename _ForwardIterator,
1343 typename _UniformRandomNumberGenerator>
1344 void
1345 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1346 _UniformRandomNumberGenerator& __urng,
1347 const param_type& __p);
1349 param_type _M_param;
1351 std::gamma_distribution<result_type> _M_gd;
1352 };
1354 /**
1355 * @brief Return true if two Nakagami distributions are not equal.
1356 */
1357 template<typename _RealType>
1358 inline bool
1359 operator!=(const nakagami_distribution<_RealType>& __d1,
1360 const nakagami_distribution<_RealType>& __d2)
1361 { return !(__d1 == __d2); }
1364 /**
1365 * @brief A Pareto continuous distribution for random numbers.
1366 *
1367 * The formula for the Pareto cumulative probability function is
1368 * @f[
1369 * P(x|\alpha,\mu) = 1 - \left(\frac{\mu}{x}\right)^\alpha
1370 * @f]
1371 * The formula for the Pareto probability density function is
1372 * @f[
1373 * p(x|\alpha,\mu) = \frac{\alpha + 1}{\mu}
1374 * \left(\frac{\mu}{x}\right)^{\alpha + 1}
1375 * @f]
1376 * where @f$x >= \mu@f$ and @f$\mu > 0@f$, @f$\alpha > 0@f$.
1377 *
1378 * <table border=1 cellpadding=10 cellspacing=0>
1379 * <caption align=top>Distribution Statistics</caption>
1380 * <tr><td>Mean</td><td>@f$\alpha \mu / (\alpha - 1)@f$
1381 * for @f$\alpha > 1@f$</td></tr>
1382 * <tr><td>Variance</td><td>@f$\alpha \mu^2 / [(\alpha - 1)^2(\alpha - 2)]@f$
1383 * for @f$\alpha > 2@f$</td></tr>
1384 * <tr><td>Range</td><td>@f$[\mu, \infty)@f$</td></tr>
1385 * </table>
1386 */
1387 template<typename _RealType = double>
1388 class
1389 pareto_distribution
1390 {
1391 static_assert(std::is_floating_point<_RealType>::value,
1392 "template argument not a floating point type");
1394 public:
1395 /** The type of the range of the distribution. */
1396 typedef _RealType result_type;
1397 /** Parameter type. */
1398 struct param_type
1399 {
1400 typedef pareto_distribution<result_type> distribution_type;
1402 param_type(result_type __alpha_val = result_type(1),
1403 result_type __mu_val = result_type(1))
1404 : _M_alpha(__alpha_val), _M_mu(__mu_val)
1405 {
1406 _GLIBCXX_DEBUG_ASSERT(_M_alpha > result_type(0));
1407 _GLIBCXX_DEBUG_ASSERT(_M_mu > result_type(0));
1408 }
1410 result_type
1411 alpha() const
1412 { return _M_alpha; }
1414 result_type
1415 mu() const
1416 { return _M_mu; }
1418 friend bool
1419 operator==(const param_type& __p1, const param_type& __p2)
1420 { return __p1._M_alpha == __p2._M_alpha && __p1._M_mu == __p2._M_mu; }
1422 private:
1423 void _M_initialize();
1425 result_type _M_alpha;
1426 result_type _M_mu;
1427 };
1429 /**
1430 * @brief Constructors.
1431 */
1432 explicit
1433 pareto_distribution(result_type __alpha_val = result_type(1),
1434 result_type __mu_val = result_type(1))
1435 : _M_param(__alpha_val, __mu_val),
1436 _M_ud()
1437 { }
1439 explicit
1440 pareto_distribution(const param_type& __p)
1441 : _M_param(__p),
1442 _M_ud()
1443 { }
1445 /**
1446 * @brief Resets the distribution state.
1447 */
1448 void
1449 reset()
1450 {
1451 _M_ud.reset();
1452 }
1454 /**
1455 * @brief Return the parameters of the distribution.
1456 */
1457 result_type
1458 alpha() const
1459 { return _M_param.alpha(); }
1461 result_type
1462 mu() const
1463 { return _M_param.mu(); }
1465 /**
1466 * @brief Returns the parameter set of the distribution.
1467 */
1468 param_type
1469 param() const
1470 { return _M_param; }
1472 /**
1473 * @brief Sets the parameter set of the distribution.
1474 * @param __param The new parameter set of the distribution.
1475 */
1476 void
1477 param(const param_type& __param)
1478 { _M_param = __param; }
1480 /**
1481 * @brief Returns the greatest lower bound value of the distribution.
1482 */
1483 result_type
1484 min() const
1485 { return this->mu(); }
1487 /**
1488 * @brief Returns the least upper bound value of the distribution.
1489 */
1490 result_type
1491 max() const
1492 { return std::numeric_limits<result_type>::max(); }
1494 /**
1495 * @brief Generating functions.
1496 */
1497 template<typename _UniformRandomNumberGenerator>
1498 result_type
1499 operator()(_UniformRandomNumberGenerator& __urng)
1500 {
1501 return this->mu() * std::pow(this->_M_ud(__urng),
1502 -result_type(1) / this->alpha());
1503 }
1505 template<typename _UniformRandomNumberGenerator>
1506 result_type
1507 operator()(_UniformRandomNumberGenerator& __urng,
1508 const param_type& __p)
1509 {
1510 return __p.mu() * std::pow(this->_M_ud(__urng),
1511 -result_type(1) / __p.alpha());
1512 }
1514 template<typename _ForwardIterator,
1515 typename _UniformRandomNumberGenerator>
1516 void
1517 __generate(_ForwardIterator __f, _ForwardIterator __t,
1518 _UniformRandomNumberGenerator& __urng)
1519 { this->__generate(__f, __t, __urng, this->param()); }
1521 template<typename _ForwardIterator,
1522 typename _UniformRandomNumberGenerator>
1523 void
1524 __generate(_ForwardIterator __f, _ForwardIterator __t,
1525 _UniformRandomNumberGenerator& __urng,
1526 const param_type& __p)
1527 { this->__generate_impl(__f, __t, __urng, __p); }
1529 template<typename _UniformRandomNumberGenerator>
1530 void
1531 __generate(result_type* __f, result_type* __t,
1532 _UniformRandomNumberGenerator& __urng,
1533 const param_type& __p)
1534 { this->__generate_impl(__f, __t, __urng, __p); }
1536 /**
1537 * @brief Return true if two Pareto distributions have
1538 * the same parameters and the sequences that would
1539 * be generated are equal.
1540 */
1541 friend bool
1542 operator==(const pareto_distribution& __d1,
1543 const pareto_distribution& __d2)
1544 { return (__d1.param() == __d2.param()
1545 && __d1._M_ud == __d2._M_ud); }
1547 /**
1548 * @brief Inserts a %pareto_distribution random number distribution
1549 * @p __x into the output stream @p __os.
1550 *
1551 * @param __os An output stream.
1552 * @param __x A %pareto_distribution random number distribution.
1553 *
1554 * @returns The output stream with the state of @p __x inserted or in
1555 * an error state.
1556 */
1557 template<typename _RealType1, typename _CharT, typename _Traits>
1558 friend std::basic_ostream<_CharT, _Traits>&
1559 operator<<(std::basic_ostream<_CharT, _Traits>&,
1560 const pareto_distribution<_RealType1>&);
1562 /**
1563 * @brief Extracts a %pareto_distribution random number distribution
1564 * @p __x from the input stream @p __is.
1565 *
1566 * @param __is An input stream.
1567 * @param __x A %pareto_distribution random number
1568 * generator engine.
1569 *
1570 * @returns The input stream with @p __x extracted or in an error state.
1571 */
1572 template<typename _RealType1, typename _CharT, typename _Traits>
1573 friend std::basic_istream<_CharT, _Traits>&
1574 operator>>(std::basic_istream<_CharT, _Traits>&,
1575 pareto_distribution<_RealType1>&);
1577 private:
1578 template<typename _ForwardIterator,
1579 typename _UniformRandomNumberGenerator>
1580 void
1581 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1582 _UniformRandomNumberGenerator& __urng,
1583 const param_type& __p);
1585 param_type _M_param;
1587 std::uniform_real_distribution<result_type> _M_ud;
1588 };
1590 /**
1591 * @brief Return true if two Pareto distributions are not equal.
1592 */
1593 template<typename _RealType>
1594 inline bool
1595 operator!=(const pareto_distribution<_RealType>& __d1,
1596 const pareto_distribution<_RealType>& __d2)
1597 { return !(__d1 == __d2); }
1600 /**
1601 * @brief A K continuous distribution for random numbers.
1602 *
1603 * The formula for the K probability density function is
1604 * @f[
1605 * p(x|\lambda, \mu, \nu) = \frac{2}{x}
1606 * \left(\frac{\lambda\nu x}{\mu}\right)^{\frac{\lambda + \nu}{2}}
1607 * \frac{1}{\Gamma(\lambda)\Gamma(\nu)}
1608 * K_{\nu - \lambda}\left(2\sqrt{\frac{\lambda\nu x}{\mu}}\right)
1609 * @f]
1610 * where @f$I_0(z)@f$ is the modified Bessel function of the second kind
1611 * of order @f$\nu - \lambda@f$ and @f$\lambda > 0@f$, @f$\mu > 0@f$
1612 * and @f$\nu > 0@f$.
1613 *
1614 * <table border=1 cellpadding=10 cellspacing=0>
1615 * <caption align=top>Distribution Statistics</caption>
1616 * <tr><td>Mean</td><td>@f$\mu@f$</td></tr>
1617 * <tr><td>Variance</td><td>@f$\mu^2\frac{\lambda + \nu + 1}{\lambda\nu}@f$</td></tr>
1618 * <tr><td>Range</td><td>@f$[0, \infty)@f$</td></tr>
1619 * </table>
1620 */
1621 template<typename _RealType = double>
1622 class
1623 k_distribution
1624 {
1625 static_assert(std::is_floating_point<_RealType>::value,
1626 "template argument not a floating point type");
1628 public:
1629 /** The type of the range of the distribution. */
1630 typedef _RealType result_type;
1631 /** Parameter type. */
1632 struct param_type
1633 {
1634 typedef k_distribution<result_type> distribution_type;
1636 param_type(result_type __lambda_val = result_type(1),
1637 result_type __mu_val = result_type(1),
1638 result_type __nu_val = result_type(1))
1639 : _M_lambda(__lambda_val), _M_mu(__mu_val), _M_nu(__nu_val)
1640 {
1641 _GLIBCXX_DEBUG_ASSERT(_M_lambda > result_type(0));
1642 _GLIBCXX_DEBUG_ASSERT(_M_mu > result_type(0));
1643 _GLIBCXX_DEBUG_ASSERT(_M_nu > result_type(0));
1644 }
1646 result_type
1647 lambda() const
1648 { return _M_lambda; }
1650 result_type
1651 mu() const
1652 { return _M_mu; }
1654 result_type
1655 nu() const
1656 { return _M_nu; }
1658 friend bool
1659 operator==(const param_type& __p1, const param_type& __p2)
1660 { return __p1._M_lambda == __p2._M_lambda
1661 && __p1._M_mu == __p2._M_mu
1662 && __p1._M_nu == __p2._M_nu; }
1664 private:
1665 void _M_initialize();
1667 result_type _M_lambda;
1668 result_type _M_mu;
1669 result_type _M_nu;
1670 };
1672 /**
1673 * @brief Constructors.
1674 */
1675 explicit
1676 k_distribution(result_type __lambda_val = result_type(1),
1677 result_type __mu_val = result_type(1),
1678 result_type __nu_val = result_type(1))
1679 : _M_param(__lambda_val, __mu_val, __nu_val),
1680 _M_gd1(__lambda_val, result_type(1) / __lambda_val),
1681 _M_gd2(__nu_val, __mu_val / __nu_val)
1682 { }
1684 explicit
1685 k_distribution(const param_type& __p)
1686 : _M_param(__p),
1687 _M_gd1(__p.lambda(), result_type(1) / __p.lambda()),
1688 _M_gd2(__p.nu(), __p.mu() / __p.nu())
1689 { }
1691 /**
1692 * @brief Resets the distribution state.
1693 */
1694 void
1695 reset()
1696 {
1697 _M_gd1.reset();
1698 _M_gd2.reset();
1699 }
1701 /**
1702 * @brief Return the parameters of the distribution.
1703 */
1704 result_type
1705 lambda() const
1706 { return _M_param.lambda(); }
1708 result_type
1709 mu() const
1710 { return _M_param.mu(); }
1712 result_type
1713 nu() const
1714 { return _M_param.nu(); }
1716 /**
1717 * @brief Returns the parameter set of the distribution.
1718 */
1719 param_type
1720 param() const
1721 { return _M_param; }
1723 /**
1724 * @brief Sets the parameter set of the distribution.
1725 * @param __param The new parameter set of the distribution.
1726 */
1727 void
1728 param(const param_type& __param)
1729 { _M_param = __param; }
1731 /**
1732 * @brief Returns the greatest lower bound value of the distribution.
1733 */
1734 result_type
1735 min() const
1736 { return result_type(0); }
1738 /**
1739 * @brief Returns the least upper bound value of the distribution.
1740 */
1741 result_type
1742 max() const
1743 { return std::numeric_limits<result_type>::max(); }
1745 /**
1746 * @brief Generating functions.
1747 */
1748 template<typename _UniformRandomNumberGenerator>
1749 result_type
1750 operator()(_UniformRandomNumberGenerator&);
1752 template<typename _UniformRandomNumberGenerator>
1753 result_type
1754 operator()(_UniformRandomNumberGenerator&, const param_type&);
1756 template<typename _ForwardIterator,
1757 typename _UniformRandomNumberGenerator>
1758 void
1759 __generate(_ForwardIterator __f, _ForwardIterator __t,
1760 _UniformRandomNumberGenerator& __urng)
1761 { this->__generate(__f, __t, __urng, this->param()); }
1763 template<typename _ForwardIterator,
1764 typename _UniformRandomNumberGenerator>
1765 void
1766 __generate(_ForwardIterator __f, _ForwardIterator __t,
1767 _UniformRandomNumberGenerator& __urng,
1768 const param_type& __p)
1769 { this->__generate_impl(__f, __t, __urng, __p); }
1771 template<typename _UniformRandomNumberGenerator>
1772 void
1773 __generate(result_type* __f, result_type* __t,
1774 _UniformRandomNumberGenerator& __urng,
1775 const param_type& __p)
1776 { this->__generate_impl(__f, __t, __urng, __p); }
1778 /**
1779 * @brief Return true if two K distributions have
1780 * the same parameters and the sequences that would
1781 * be generated are equal.
1782 */
1783 friend bool
1784 operator==(const k_distribution& __d1,
1785 const k_distribution& __d2)
1786 { return (__d1.param() == __d2.param()
1787 && __d1._M_gd1 == __d2._M_gd1
1788 && __d1._M_gd2 == __d2._M_gd2); }
1790 /**
1791 * @brief Inserts a %k_distribution random number distribution
1792 * @p __x into the output stream @p __os.
1793 *
1794 * @param __os An output stream.
1795 * @param __x A %k_distribution random number distribution.
1796 *
1797 * @returns The output stream with the state of @p __x inserted or in
1798 * an error state.
1799 */
1800 template<typename _RealType1, typename _CharT, typename _Traits>
1801 friend std::basic_ostream<_CharT, _Traits>&
1802 operator<<(std::basic_ostream<_CharT, _Traits>&,
1803 const k_distribution<_RealType1>&);
1805 /**
1806 * @brief Extracts a %k_distribution random number distribution
1807 * @p __x from the input stream @p __is.
1808 *
1809 * @param __is An input stream.
1810 * @param __x A %k_distribution random number
1811 * generator engine.
1812 *
1813 * @returns The input stream with @p __x extracted or in an error state.
1814 */
1815 template<typename _RealType1, typename _CharT, typename _Traits>
1816 friend std::basic_istream<_CharT, _Traits>&
1817 operator>>(std::basic_istream<_CharT, _Traits>&,
1818 k_distribution<_RealType1>&);
1820 private:
1821 template<typename _ForwardIterator,
1822 typename _UniformRandomNumberGenerator>
1823 void
1824 __generate_impl(_ForwardIterator __f, _ForwardIterator __t,
1825 _UniformRandomNumberGenerator& __urng,
1826 const param_type& __p);
1828 param_type _M_param;
1830 std::gamma_distribution<result_type> _M_gd1;
1831 std::gamma_distribution<result_type> _M_gd2;
1832 };
1834 /**
1835 * @brief Return true if two K distributions are not equal.
1836 */
1837 template<typename _RealType>
1838 inline bool
1839 operator!=(const k_distribution<_RealType>& __d1,
1840 const k_distribution<_RealType>& __d2)
1841 { return !(__d1 == __d2); }
1843 _GLIBCXX_END_NAMESPACE_VERSION
1844 } // namespace __gnu_cxx
1846 #include "opt_random.h"
1847 #include "random.tcc"
1849 #endif /* _EXT_RANDOM */