STLdoc
STLdocumentation
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Macros | Functions
random.tcc File Reference

Macros

#define _GLIBCXX_TR1_RANDOM_TCC   1
 

Functions

namespace std _GLIBCXX_VISIBILITY (default)
 

Detailed Description

This is an internal header file, included by other library headers. Do not attempt to use it directly. {tr1/random}

Macro Definition Documentation

#define _GLIBCXX_TR1_RANDOM_TCC   1

Function Documentation

namespace std _GLIBCXX_VISIBILITY ( default  )

Seeds the LCR with integral value __x0, adjusted so that the ring identity is never a member of the convergence set.

Seeds the LCR engine with a value generated by __g.

Gets the next generated value in sequence.

A rejection algorithm when mean >= 12 and a simple method based upon the multiplication of uniform random variates otherwise. NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1 is defined.

Reference: Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).

A rejection algorithm when t * p >= 8 and a simple waiting time method - the second in the referenced book - otherwise. NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1 is defined.

Reference: Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, New York, 1986, Ch. X, Sect. 4 (+ Errata!).

Polar method due to Marsaglia.

Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, New York, 1986, Ch. V, Sect. 4.4.

Cheng's rejection algorithm GB for alpha >= 1 and a modification of Vaduva's rejection from Weibull algorithm due to Devroye for alpha < 1.

References: Cheng, R. C. The Generation of Gamma Random Variables with Non-integral Shape Parameter. Applied Statistics, 26, 71-75, 1977.

Vaduva, I. Computer Generation of Gamma Gandom Variables by Rejection and Composition Procedures. Math. Operationsforschung and Statistik, Series in Statistics, 8, 545-576, 1977.

Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, New York, 1986, Ch. IX, Sect. 3.4 (+ Errata!).

35 {
36 namespace tr1
37 {
38  /*
39  * (Further) implementation-space details.
40  */
41  namespace __detail
42  {
43  _GLIBCXX_BEGIN_NAMESPACE_VERSION
44 
45  // General case for x = (ax + c) mod m -- use Schrage's algorithm to avoid
46  // integer overflow.
47  //
48  // Because a and c are compile-time integral constants the compiler kindly
49  // elides any unreachable paths.
50  //
51  // Preconditions: a > 0, m > 0.
52  //
53  template<typename _Tp, _Tp __a, _Tp __c, _Tp __m, bool>
54  struct _Mod
55  {
56  static _Tp
57  __calc(_Tp __x)
58  {
59  if (__a == 1)
60  __x %= __m;
61  else
62  {
63  static const _Tp __q = __m / __a;
64  static const _Tp __r = __m % __a;
65 
66  _Tp __t1 = __a * (__x % __q);
67  _Tp __t2 = __r * (__x / __q);
68  if (__t1 >= __t2)
69  __x = __t1 - __t2;
70  else
71  __x = __m - __t2 + __t1;
72  }
73 
74  if (__c != 0)
75  {
76  const _Tp __d = __m - __x;
77  if (__d > __c)
78  __x += __c;
79  else
80  __x = __c - __d;
81  }
82  return __x;
83  }
84  };
85 
86  // Special case for m == 0 -- use unsigned integer overflow as modulo
87  // operator.
88  template<typename _Tp, _Tp __a, _Tp __c, _Tp __m>
89  struct _Mod<_Tp, __a, __c, __m, true>
90  {
91  static _Tp
92  __calc(_Tp __x)
93  { return __a * __x + __c; }
94  };
95  _GLIBCXX_END_NAMESPACE_VERSION
96  } // namespace __detail
97 
98 _GLIBCXX_BEGIN_NAMESPACE_VERSION
99 
100  template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
101  const _UIntType
102  linear_congruential<_UIntType, __a, __c, __m>::multiplier;
103 
104  template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
105  const _UIntType
106  linear_congruential<_UIntType, __a, __c, __m>::increment;
107 
108  template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
109  const _UIntType
110  linear_congruential<_UIntType, __a, __c, __m>::modulus;
111 
116  template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
117  void
118  linear_congruential<_UIntType, __a, __c, __m>::
119  seed(unsigned long __x0)
120  {
121  if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0)
122  && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
123  _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1);
124  else
125  _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0);
126  }
127 
131  template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
132  template<class _Gen>
133  void
134  linear_congruential<_UIntType, __a, __c, __m>::
135  seed(_Gen& __g, false_type)
136  {
137  _UIntType __x0 = __g();
138  if ((__detail::__mod<_UIntType, 1, 0, __m>(__c) == 0)
139  && (__detail::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
140  _M_x = __detail::__mod<_UIntType, 1, 0, __m>(1);
141  else
142  _M_x = __detail::__mod<_UIntType, 1, 0, __m>(__x0);
143  }
144 
148  template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
149  typename linear_congruential<_UIntType, __a, __c, __m>::result_type
150  linear_congruential<_UIntType, __a, __c, __m>::
151  operator()()
152  {
153  _M_x = __detail::__mod<_UIntType, __a, __c, __m>(_M_x);
154  return _M_x;
155  }
156 
157  template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
158  typename _CharT, typename _Traits>
159  std::basic_ostream<_CharT, _Traits>&
160  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
161  const linear_congruential<_UIntType, __a, __c, __m>& __lcr)
162  {
163  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
164  typedef typename __ostream_type::ios_base __ios_base;
165 
166  const typename __ios_base::fmtflags __flags = __os.flags();
167  const _CharT __fill = __os.fill();
168  __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
169  __os.fill(__os.widen(' '));
170 
171  __os << __lcr._M_x;
172 
173  __os.flags(__flags);
174  __os.fill(__fill);
175  return __os;
176  }
177 
178  template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
179  typename _CharT, typename _Traits>
180  std::basic_istream<_CharT, _Traits>&
181  operator>>(std::basic_istream<_CharT, _Traits>& __is,
182  linear_congruential<_UIntType, __a, __c, __m>& __lcr)
183  {
184  typedef std::basic_istream<_CharT, _Traits> __istream_type;
185  typedef typename __istream_type::ios_base __ios_base;
186 
187  const typename __ios_base::fmtflags __flags = __is.flags();
188  __is.flags(__ios_base::dec);
189 
190  __is >> __lcr._M_x;
191 
192  __is.flags(__flags);
193  return __is;
194  }
195 
196 
197  template<class _UIntType, int __w, int __n, int __m, int __r,
198  _UIntType __a, int __u, int __s,
199  _UIntType __b, int __t, _UIntType __c, int __l>
200  const int
201  mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
202  __b, __t, __c, __l>::word_size;
203 
204  template<class _UIntType, int __w, int __n, int __m, int __r,
205  _UIntType __a, int __u, int __s,
206  _UIntType __b, int __t, _UIntType __c, int __l>
207  const int
208  mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
209  __b, __t, __c, __l>::state_size;
210 
211  template<class _UIntType, int __w, int __n, int __m, int __r,
212  _UIntType __a, int __u, int __s,
213  _UIntType __b, int __t, _UIntType __c, int __l>
214  const int
215  mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
216  __b, __t, __c, __l>::shift_size;
217 
218  template<class _UIntType, int __w, int __n, int __m, int __r,
219  _UIntType __a, int __u, int __s,
220  _UIntType __b, int __t, _UIntType __c, int __l>
221  const int
222  mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
223  __b, __t, __c, __l>::mask_bits;
224 
225  template<class _UIntType, int __w, int __n, int __m, int __r,
226  _UIntType __a, int __u, int __s,
227  _UIntType __b, int __t, _UIntType __c, int __l>
228  const _UIntType
229  mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
230  __b, __t, __c, __l>::parameter_a;
231 
232  template<class _UIntType, int __w, int __n, int __m, int __r,
233  _UIntType __a, int __u, int __s,
234  _UIntType __b, int __t, _UIntType __c, int __l>
235  const int
236  mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
237  __b, __t, __c, __l>::output_u;
238 
239  template<class _UIntType, int __w, int __n, int __m, int __r,
240  _UIntType __a, int __u, int __s,
241  _UIntType __b, int __t, _UIntType __c, int __l>
242  const int
243  mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
244  __b, __t, __c, __l>::output_s;
245 
246  template<class _UIntType, int __w, int __n, int __m, int __r,
247  _UIntType __a, int __u, int __s,
248  _UIntType __b, int __t, _UIntType __c, int __l>
249  const _UIntType
250  mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
251  __b, __t, __c, __l>::output_b;
252 
253  template<class _UIntType, int __w, int __n, int __m, int __r,
254  _UIntType __a, int __u, int __s,
255  _UIntType __b, int __t, _UIntType __c, int __l>
256  const int
257  mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
258  __b, __t, __c, __l>::output_t;
259 
260  template<class _UIntType, int __w, int __n, int __m, int __r,
261  _UIntType __a, int __u, int __s,
262  _UIntType __b, int __t, _UIntType __c, int __l>
263  const _UIntType
264  mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
265  __b, __t, __c, __l>::output_c;
266 
267  template<class _UIntType, int __w, int __n, int __m, int __r,
268  _UIntType __a, int __u, int __s,
269  _UIntType __b, int __t, _UIntType __c, int __l>
270  const int
271  mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
272  __b, __t, __c, __l>::output_l;
273 
274  template<class _UIntType, int __w, int __n, int __m, int __r,
275  _UIntType __a, int __u, int __s,
276  _UIntType __b, int __t, _UIntType __c, int __l>
277  void
278  mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
279  __b, __t, __c, __l>::
280  seed(unsigned long __value)
281  {
282  _M_x[0] = __detail::__mod<_UIntType, 1, 0,
283  __detail::_Shift<_UIntType, __w>::__value>(__value);
284 
285  for (int __i = 1; __i < state_size; ++__i)
286  {
287  _UIntType __x = _M_x[__i - 1];
288  __x ^= __x >> (__w - 2);
289  __x *= 1812433253ul;
290  __x += __i;
291  _M_x[__i] = __detail::__mod<_UIntType, 1, 0,
292  __detail::_Shift<_UIntType, __w>::__value>(__x);
293  }
294  _M_p = state_size;
295  }
296 
297  template<class _UIntType, int __w, int __n, int __m, int __r,
298  _UIntType __a, int __u, int __s,
299  _UIntType __b, int __t, _UIntType __c, int __l>
300  template<class _Gen>
301  void
302  mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
303  __b, __t, __c, __l>::
304  seed(_Gen& __gen, false_type)
305  {
306  for (int __i = 0; __i < state_size; ++__i)
307  _M_x[__i] = __detail::__mod<_UIntType, 1, 0,
308  __detail::_Shift<_UIntType, __w>::__value>(__gen());
309  _M_p = state_size;
310  }
311 
312  template<class _UIntType, int __w, int __n, int __m, int __r,
313  _UIntType __a, int __u, int __s,
314  _UIntType __b, int __t, _UIntType __c, int __l>
315  typename
316  mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
317  __b, __t, __c, __l>::result_type
318  mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
319  __b, __t, __c, __l>::
320  operator()()
321  {
322  // Reload the vector - cost is O(n) amortized over n calls.
323  if (_M_p >= state_size)
324  {
325  const _UIntType __upper_mask = (~_UIntType()) << __r;
326  const _UIntType __lower_mask = ~__upper_mask;
327 
328  for (int __k = 0; __k < (__n - __m); ++__k)
329  {
330  _UIntType __y = ((_M_x[__k] & __upper_mask)
331  | (_M_x[__k + 1] & __lower_mask));
332  _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
333  ^ ((__y & 0x01) ? __a : 0));
334  }
335 
336  for (int __k = (__n - __m); __k < (__n - 1); ++__k)
337  {
338  _UIntType __y = ((_M_x[__k] & __upper_mask)
339  | (_M_x[__k + 1] & __lower_mask));
340  _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
341  ^ ((__y & 0x01) ? __a : 0));
342  }
343 
344  _UIntType __y = ((_M_x[__n - 1] & __upper_mask)
345  | (_M_x[0] & __lower_mask));
346  _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
347  ^ ((__y & 0x01) ? __a : 0));
348  _M_p = 0;
349  }
350 
351  // Calculate o(x(i)).
352  result_type __z = _M_x[_M_p++];
353  __z ^= (__z >> __u);
354  __z ^= (__z << __s) & __b;
355  __z ^= (__z << __t) & __c;
356  __z ^= (__z >> __l);
357 
358  return __z;
359  }
360 
361  template<class _UIntType, int __w, int __n, int __m, int __r,
362  _UIntType __a, int __u, int __s, _UIntType __b, int __t,
363  _UIntType __c, int __l,
364  typename _CharT, typename _Traits>
365  std::basic_ostream<_CharT, _Traits>&
366  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
367  const mersenne_twister<_UIntType, __w, __n, __m,
368  __r, __a, __u, __s, __b, __t, __c, __l>& __x)
369  {
370  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
371  typedef typename __ostream_type::ios_base __ios_base;
372 
373  const typename __ios_base::fmtflags __flags = __os.flags();
374  const _CharT __fill = __os.fill();
375  const _CharT __space = __os.widen(' ');
376  __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
377  __os.fill(__space);
378 
379  for (int __i = 0; __i < __n - 1; ++__i)
380  __os << __x._M_x[__i] << __space;
381  __os << __x._M_x[__n - 1];
382 
383  __os.flags(__flags);
384  __os.fill(__fill);
385  return __os;
386  }
387 
388  template<class _UIntType, int __w, int __n, int __m, int __r,
389  _UIntType __a, int __u, int __s, _UIntType __b, int __t,
390  _UIntType __c, int __l,
391  typename _CharT, typename _Traits>
392  std::basic_istream<_CharT, _Traits>&
393  operator>>(std::basic_istream<_CharT, _Traits>& __is,
394  mersenne_twister<_UIntType, __w, __n, __m,
395  __r, __a, __u, __s, __b, __t, __c, __l>& __x)
396  {
397  typedef std::basic_istream<_CharT, _Traits> __istream_type;
398  typedef typename __istream_type::ios_base __ios_base;
399 
400  const typename __ios_base::fmtflags __flags = __is.flags();
401  __is.flags(__ios_base::dec | __ios_base::skipws);
402 
403  for (int __i = 0; __i < __n; ++__i)
404  __is >> __x._M_x[__i];
405 
406  __is.flags(__flags);
407  return __is;
408  }
409 
410 
411  template<typename _IntType, _IntType __m, int __s, int __r>
412  const _IntType
413  subtract_with_carry<_IntType, __m, __s, __r>::modulus;
414 
415  template<typename _IntType, _IntType __m, int __s, int __r>
416  const int
417  subtract_with_carry<_IntType, __m, __s, __r>::long_lag;
418 
419  template<typename _IntType, _IntType __m, int __s, int __r>
420  const int
421  subtract_with_carry<_IntType, __m, __s, __r>::short_lag;
422 
423  template<typename _IntType, _IntType __m, int __s, int __r>
424  void
425  subtract_with_carry<_IntType, __m, __s, __r>::
426  seed(unsigned long __value)
427  {
428  if (__value == 0)
429  __value = 19780503;
430 
431  std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563>
432  __lcg(__value);
433 
434  for (int __i = 0; __i < long_lag; ++__i)
435  _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__lcg());
436 
437  _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
438  _M_p = 0;
439  }
440 
441  template<typename _IntType, _IntType __m, int __s, int __r>
442  template<class _Gen>
443  void
444  subtract_with_carry<_IntType, __m, __s, __r>::
445  seed(_Gen& __gen, false_type)
446  {
447  const int __n = (std::numeric_limits<_UIntType>::digits + 31) / 32;
448 
449  for (int __i = 0; __i < long_lag; ++__i)
450  {
451  _UIntType __tmp = 0;
452  _UIntType __factor = 1;
453  for (int __j = 0; __j < __n; ++__j)
454  {
455  __tmp += __detail::__mod<__detail::_UInt32Type, 1, 0, 0>
456  (__gen()) * __factor;
457  __factor *= __detail::_Shift<_UIntType, 32>::__value;
458  }
459  _M_x[__i] = __detail::__mod<_UIntType, 1, 0, modulus>(__tmp);
460  }
461  _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
462  _M_p = 0;
463  }
464 
465  template<typename _IntType, _IntType __m, int __s, int __r>
466  typename subtract_with_carry<_IntType, __m, __s, __r>::result_type
467  subtract_with_carry<_IntType, __m, __s, __r>::
468  operator()()
469  {
470  // Derive short lag index from current index.
471  int __ps = _M_p - short_lag;
472  if (__ps < 0)
473  __ps += long_lag;
474 
475  // Calculate new x(i) without overflow or division.
476  // NB: Thanks to the requirements for _IntType, _M_x[_M_p] + _M_carry
477  // cannot overflow.
478  _UIntType __xi;
479  if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
480  {
481  __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
482  _M_carry = 0;
483  }
484  else
485  {
486  __xi = modulus - _M_x[_M_p] - _M_carry + _M_x[__ps];
487  _M_carry = 1;
488  }
489  _M_x[_M_p] = __xi;
490 
491  // Adjust current index to loop around in ring buffer.
492  if (++_M_p >= long_lag)
493  _M_p = 0;
494 
495  return __xi;
496  }
497 
498  template<typename _IntType, _IntType __m, int __s, int __r,
499  typename _CharT, typename _Traits>
500  std::basic_ostream<_CharT, _Traits>&
501  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
502  const subtract_with_carry<_IntType, __m, __s, __r>& __x)
503  {
504  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
505  typedef typename __ostream_type::ios_base __ios_base;
506 
507  const typename __ios_base::fmtflags __flags = __os.flags();
508  const _CharT __fill = __os.fill();
509  const _CharT __space = __os.widen(' ');
510  __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
511  __os.fill(__space);
512 
513  for (int __i = 0; __i < __r; ++__i)
514  __os << __x._M_x[__i] << __space;
515  __os << __x._M_carry;
516 
517  __os.flags(__flags);
518  __os.fill(__fill);
519  return __os;
520  }
521 
522  template<typename _IntType, _IntType __m, int __s, int __r,
523  typename _CharT, typename _Traits>
524  std::basic_istream<_CharT, _Traits>&
525  operator>>(std::basic_istream<_CharT, _Traits>& __is,
526  subtract_with_carry<_IntType, __m, __s, __r>& __x)
527  {
528  typedef std::basic_ostream<_CharT, _Traits> __istream_type;
529  typedef typename __istream_type::ios_base __ios_base;
530 
531  const typename __ios_base::fmtflags __flags = __is.flags();
532  __is.flags(__ios_base::dec | __ios_base::skipws);
533 
534  for (int __i = 0; __i < __r; ++__i)
535  __is >> __x._M_x[__i];
536  __is >> __x._M_carry;
537 
538  __is.flags(__flags);
539  return __is;
540  }
541 
542 
543  template<typename _RealType, int __w, int __s, int __r>
544  const int
545  subtract_with_carry_01<_RealType, __w, __s, __r>::word_size;
546 
547  template<typename _RealType, int __w, int __s, int __r>
548  const int
549  subtract_with_carry_01<_RealType, __w, __s, __r>::long_lag;
550 
551  template<typename _RealType, int __w, int __s, int __r>
552  const int
553  subtract_with_carry_01<_RealType, __w, __s, __r>::short_lag;
554 
555  template<typename _RealType, int __w, int __s, int __r>
556  void
557  subtract_with_carry_01<_RealType, __w, __s, __r>::
558  _M_initialize_npows()
559  {
560  for (int __j = 0; __j < __n; ++__j)
561 #if _GLIBCXX_USE_C99_MATH_TR1
562  _M_npows[__j] = std::tr1::ldexp(_RealType(1), -__w + __j * 32);
563 #else
564  _M_npows[__j] = std::pow(_RealType(2), -__w + __j * 32);
565 #endif
566  }
567 
568  template<typename _RealType, int __w, int __s, int __r>
569  void
570  subtract_with_carry_01<_RealType, __w, __s, __r>::
571  seed(unsigned long __value)
572  {
573  if (__value == 0)
574  __value = 19780503;
575 
576  // _GLIBCXX_RESOLVE_LIB_DEFECTS
577  // 512. Seeding subtract_with_carry_01 from a single unsigned long.
578  std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563>
579  __lcg(__value);
580 
581  this->seed(__lcg);
582  }
583 
584  template<typename _RealType, int __w, int __s, int __r>
585  template<class _Gen>
586  void
587  subtract_with_carry_01<_RealType, __w, __s, __r>::
588  seed(_Gen& __gen, false_type)
589  {
590  for (int __i = 0; __i < long_lag; ++__i)
591  {
592  for (int __j = 0; __j < __n - 1; ++__j)
593  _M_x[__i][__j] = __detail::__mod<_UInt32Type, 1, 0, 0>(__gen());
594  _M_x[__i][__n - 1] = __detail::__mod<_UInt32Type, 1, 0,
595  __detail::_Shift<_UInt32Type, __w % 32>::__value>(__gen());
596  }
597 
598  _M_carry = 1;
599  for (int __j = 0; __j < __n; ++__j)
600  if (_M_x[long_lag - 1][__j] != 0)
601  {
602  _M_carry = 0;
603  break;
604  }
605 
606  _M_p = 0;
607  }
608 
609  template<typename _RealType, int __w, int __s, int __r>
610  typename subtract_with_carry_01<_RealType, __w, __s, __r>::result_type
611  subtract_with_carry_01<_RealType, __w, __s, __r>::
612  operator()()
613  {
614  // Derive short lag index from current index.
615  int __ps = _M_p - short_lag;
616  if (__ps < 0)
617  __ps += long_lag;
618 
619  _UInt32Type __new_carry;
620  for (int __j = 0; __j < __n - 1; ++__j)
621  {
622  if (_M_x[__ps][__j] > _M_x[_M_p][__j]
623  || (_M_x[__ps][__j] == _M_x[_M_p][__j] && _M_carry == 0))
624  __new_carry = 0;
625  else
626  __new_carry = 1;
627 
628  _M_x[_M_p][__j] = _M_x[__ps][__j] - _M_x[_M_p][__j] - _M_carry;
629  _M_carry = __new_carry;
630  }
631 
632  if (_M_x[__ps][__n - 1] > _M_x[_M_p][__n - 1]
633  || (_M_x[__ps][__n - 1] == _M_x[_M_p][__n - 1] && _M_carry == 0))
634  __new_carry = 0;
635  else
636  __new_carry = 1;
637 
638  _M_x[_M_p][__n - 1] = __detail::__mod<_UInt32Type, 1, 0,
639  __detail::_Shift<_UInt32Type, __w % 32>::__value>
640  (_M_x[__ps][__n - 1] - _M_x[_M_p][__n - 1] - _M_carry);
641  _M_carry = __new_carry;
642 
643  result_type __ret = 0.0;
644  for (int __j = 0; __j < __n; ++__j)
645  __ret += _M_x[_M_p][__j] * _M_npows[__j];
646 
647  // Adjust current index to loop around in ring buffer.
648  if (++_M_p >= long_lag)
649  _M_p = 0;
650 
651  return __ret;
652  }
653 
654  template<typename _RealType, int __w, int __s, int __r,
655  typename _CharT, typename _Traits>
656  std::basic_ostream<_CharT, _Traits>&
657  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
658  const subtract_with_carry_01<_RealType, __w, __s, __r>& __x)
659  {
660  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
661  typedef typename __ostream_type::ios_base __ios_base;
662 
663  const typename __ios_base::fmtflags __flags = __os.flags();
664  const _CharT __fill = __os.fill();
665  const _CharT __space = __os.widen(' ');
666  __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
667  __os.fill(__space);
668 
669  for (int __i = 0; __i < __r; ++__i)
670  for (int __j = 0; __j < __x.__n; ++__j)
671  __os << __x._M_x[__i][__j] << __space;
672  __os << __x._M_carry;
673 
674  __os.flags(__flags);
675  __os.fill(__fill);
676  return __os;
677  }
678 
679  template<typename _RealType, int __w, int __s, int __r,
680  typename _CharT, typename _Traits>
681  std::basic_istream<_CharT, _Traits>&
682  operator>>(std::basic_istream<_CharT, _Traits>& __is,
683  subtract_with_carry_01<_RealType, __w, __s, __r>& __x)
684  {
685  typedef std::basic_istream<_CharT, _Traits> __istream_type;
686  typedef typename __istream_type::ios_base __ios_base;
687 
688  const typename __ios_base::fmtflags __flags = __is.flags();
689  __is.flags(__ios_base::dec | __ios_base::skipws);
690 
691  for (int __i = 0; __i < __r; ++__i)
692  for (int __j = 0; __j < __x.__n; ++__j)
693  __is >> __x._M_x[__i][__j];
694  __is >> __x._M_carry;
695 
696  __is.flags(__flags);
697  return __is;
698  }
699 
700  template<class _UniformRandomNumberGenerator, int __p, int __r>
701  const int
702  discard_block<_UniformRandomNumberGenerator, __p, __r>::block_size;
703 
704  template<class _UniformRandomNumberGenerator, int __p, int __r>
705  const int
706  discard_block<_UniformRandomNumberGenerator, __p, __r>::used_block;
707 
708  template<class _UniformRandomNumberGenerator, int __p, int __r>
709  typename discard_block<_UniformRandomNumberGenerator,
710  __p, __r>::result_type
711  discard_block<_UniformRandomNumberGenerator, __p, __r>::
712  operator()()
713  {
714  if (_M_n >= used_block)
715  {
716  while (_M_n < block_size)
717  {
718  _M_b();
719  ++_M_n;
720  }
721  _M_n = 0;
722  }
723  ++_M_n;
724  return _M_b();
725  }
726 
727  template<class _UniformRandomNumberGenerator, int __p, int __r,
728  typename _CharT, typename _Traits>
729  std::basic_ostream<_CharT, _Traits>&
730  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
731  const discard_block<_UniformRandomNumberGenerator,
732  __p, __r>& __x)
733  {
734  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
735  typedef typename __ostream_type::ios_base __ios_base;
736 
737  const typename __ios_base::fmtflags __flags = __os.flags();
738  const _CharT __fill = __os.fill();
739  const _CharT __space = __os.widen(' ');
740  __os.flags(__ios_base::dec | __ios_base::fixed
741  | __ios_base::left);
742  __os.fill(__space);
743 
744  __os << __x._M_b << __space << __x._M_n;
745 
746  __os.flags(__flags);
747  __os.fill(__fill);
748  return __os;
749  }
750 
751  template<class _UniformRandomNumberGenerator, int __p, int __r,
752  typename _CharT, typename _Traits>
753  std::basic_istream<_CharT, _Traits>&
754  operator>>(std::basic_istream<_CharT, _Traits>& __is,
755  discard_block<_UniformRandomNumberGenerator, __p, __r>& __x)
756  {
757  typedef std::basic_istream<_CharT, _Traits> __istream_type;
758  typedef typename __istream_type::ios_base __ios_base;
759 
760  const typename __ios_base::fmtflags __flags = __is.flags();
761  __is.flags(__ios_base::dec | __ios_base::skipws);
762 
763  __is >> __x._M_b >> __x._M_n;
764 
765  __is.flags(__flags);
766  return __is;
767  }
768 
769 
770  template<class _UniformRandomNumberGenerator1, int __s1,
771  class _UniformRandomNumberGenerator2, int __s2>
772  const int
773  xor_combine<_UniformRandomNumberGenerator1, __s1,
774  _UniformRandomNumberGenerator2, __s2>::shift1;
775 
776  template<class _UniformRandomNumberGenerator1, int __s1,
777  class _UniformRandomNumberGenerator2, int __s2>
778  const int
779  xor_combine<_UniformRandomNumberGenerator1, __s1,
780  _UniformRandomNumberGenerator2, __s2>::shift2;
781 
782  template<class _UniformRandomNumberGenerator1, int __s1,
783  class _UniformRandomNumberGenerator2, int __s2>
784  void
785  xor_combine<_UniformRandomNumberGenerator1, __s1,
786  _UniformRandomNumberGenerator2, __s2>::
787  _M_initialize_max()
788  {
789  const int __w = std::numeric_limits<result_type>::digits;
790 
791  const result_type __m1 =
792  std::min(result_type(_M_b1.max() - _M_b1.min()),
793  __detail::_Shift<result_type, __w - __s1>::__value - 1);
794 
795  const result_type __m2 =
796  std::min(result_type(_M_b2.max() - _M_b2.min()),
797  __detail::_Shift<result_type, __w - __s2>::__value - 1);
798 
799  // NB: In TR1 s1 is not required to be >= s2.
800  if (__s1 < __s2)
801  _M_max = _M_initialize_max_aux(__m2, __m1, __s2 - __s1) << __s1;
802  else
803  _M_max = _M_initialize_max_aux(__m1, __m2, __s1 - __s2) << __s2;
804  }
805 
806  template<class _UniformRandomNumberGenerator1, int __s1,
807  class _UniformRandomNumberGenerator2, int __s2>
808  typename xor_combine<_UniformRandomNumberGenerator1, __s1,
809  _UniformRandomNumberGenerator2, __s2>::result_type
810  xor_combine<_UniformRandomNumberGenerator1, __s1,
811  _UniformRandomNumberGenerator2, __s2>::
812  _M_initialize_max_aux(result_type __a, result_type __b, int __d)
813  {
814  const result_type __two2d = result_type(1) << __d;
815  const result_type __c = __a * __two2d;
816 
817  if (__a == 0 || __b < __two2d)
818  return __c + __b;
819 
820  const result_type __t = std::max(__c, __b);
821  const result_type __u = std::min(__c, __b);
822 
823  result_type __ub = __u;
824  result_type __p;
825  for (__p = 0; __ub != 1; __ub >>= 1)
826  ++__p;
827 
828  const result_type __two2p = result_type(1) << __p;
829  const result_type __k = __t / __two2p;
830 
831  if (__k & 1)
832  return (__k + 1) * __two2p - 1;
833 
834  if (__c >= __b)
835  return (__k + 1) * __two2p + _M_initialize_max_aux((__t % __two2p)
836  / __two2d,
837  __u % __two2p, __d);
838  else
839  return (__k + 1) * __two2p + _M_initialize_max_aux((__u % __two2p)
840  / __two2d,
841  __t % __two2p, __d);
842  }
843 
844  template<class _UniformRandomNumberGenerator1, int __s1,
845  class _UniformRandomNumberGenerator2, int __s2,
846  typename _CharT, typename _Traits>
847  std::basic_ostream<_CharT, _Traits>&
848  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
849  const xor_combine<_UniformRandomNumberGenerator1, __s1,
850  _UniformRandomNumberGenerator2, __s2>& __x)
851  {
852  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
853  typedef typename __ostream_type::ios_base __ios_base;
854 
855  const typename __ios_base::fmtflags __flags = __os.flags();
856  const _CharT __fill = __os.fill();
857  const _CharT __space = __os.widen(' ');
858  __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
859  __os.fill(__space);
860 
861  __os << __x.base1() << __space << __x.base2();
862 
863  __os.flags(__flags);
864  __os.fill(__fill);
865  return __os;
866  }
867 
868  template<class _UniformRandomNumberGenerator1, int __s1,
869  class _UniformRandomNumberGenerator2, int __s2,
870  typename _CharT, typename _Traits>
871  std::basic_istream<_CharT, _Traits>&
872  operator>>(std::basic_istream<_CharT, _Traits>& __is,
873  xor_combine<_UniformRandomNumberGenerator1, __s1,
874  _UniformRandomNumberGenerator2, __s2>& __x)
875  {
876  typedef std::basic_istream<_CharT, _Traits> __istream_type;
877  typedef typename __istream_type::ios_base __ios_base;
878 
879  const typename __ios_base::fmtflags __flags = __is.flags();
880  __is.flags(__ios_base::skipws);
881 
882  __is >> __x._M_b1 >> __x._M_b2;
883 
884  __is.flags(__flags);
885  return __is;
886  }
887 
888 
889  template<typename _IntType>
890  template<typename _UniformRandomNumberGenerator>
891  typename uniform_int<_IntType>::result_type
892  uniform_int<_IntType>::
893  _M_call(_UniformRandomNumberGenerator& __urng,
894  result_type __min, result_type __max, true_type)
895  {
896  // XXX Must be fixed to work well for *arbitrary* __urng.max(),
897  // __urng.min(), __max, __min. Currently works fine only in the
898  // most common case __urng.max() - __urng.min() >= __max - __min,
899  // with __urng.max() > __urng.min() >= 0.
900  typedef typename __gnu_cxx::__add_unsigned<typename
901  _UniformRandomNumberGenerator::result_type>::__type __urntype;
902  typedef typename __gnu_cxx::__add_unsigned<result_type>::__type
903  __utype;
904  typedef typename __gnu_cxx::__conditional_type<(sizeof(__urntype)
905  > sizeof(__utype)),
906  __urntype, __utype>::__type __uctype;
907 
908  result_type __ret;
909 
910  const __urntype __urnmin = __urng.min();
911  const __urntype __urnmax = __urng.max();
912  const __urntype __urnrange = __urnmax - __urnmin;
913  const __uctype __urange = __max - __min;
914  const __uctype __udenom = (__urnrange <= __urange
915  ? 1 : __urnrange / (__urange + 1));
916  do
917  __ret = (__urntype(__urng()) - __urnmin) / __udenom;
918  while (__ret > __max - __min);
919 
920  return __ret + __min;
921  }
922 
923  template<typename _IntType, typename _CharT, typename _Traits>
924  std::basic_ostream<_CharT, _Traits>&
925  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
926  const uniform_int<_IntType>& __x)
927  {
928  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
929  typedef typename __ostream_type::ios_base __ios_base;
930 
931  const typename __ios_base::fmtflags __flags = __os.flags();
932  const _CharT __fill = __os.fill();
933  const _CharT __space = __os.widen(' ');
934  __os.flags(__ios_base::scientific | __ios_base::left);
935  __os.fill(__space);
936 
937  __os << __x.min() << __space << __x.max();
938 
939  __os.flags(__flags);
940  __os.fill(__fill);
941  return __os;
942  }
943 
944  template<typename _IntType, typename _CharT, typename _Traits>
945  std::basic_istream<_CharT, _Traits>&
946  operator>>(std::basic_istream<_CharT, _Traits>& __is,
947  uniform_int<_IntType>& __x)
948  {
949  typedef std::basic_istream<_CharT, _Traits> __istream_type;
950  typedef typename __istream_type::ios_base __ios_base;
951 
952  const typename __ios_base::fmtflags __flags = __is.flags();
953  __is.flags(__ios_base::dec | __ios_base::skipws);
954 
955  __is >> __x._M_min >> __x._M_max;
956 
957  __is.flags(__flags);
958  return __is;
959  }
960 
961 
962  template<typename _CharT, typename _Traits>
963  std::basic_ostream<_CharT, _Traits>&
964  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
965  const bernoulli_distribution& __x)
966  {
967  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
968  typedef typename __ostream_type::ios_base __ios_base;
969 
970  const typename __ios_base::fmtflags __flags = __os.flags();
971  const _CharT __fill = __os.fill();
972  const std::streamsize __precision = __os.precision();
973  __os.flags(__ios_base::scientific | __ios_base::left);
974  __os.fill(__os.widen(' '));
975  __os.precision(__gnu_cxx::__numeric_traits<double>::__max_digits10);
976 
977  __os << __x.p();
978 
979  __os.flags(__flags);
980  __os.fill(__fill);
981  __os.precision(__precision);
982  return __os;
983  }
984 
985 
986  template<typename _IntType, typename _RealType>
987  template<class _UniformRandomNumberGenerator>
988  typename geometric_distribution<_IntType, _RealType>::result_type
989  geometric_distribution<_IntType, _RealType>::
990  operator()(_UniformRandomNumberGenerator& __urng)
991  {
992  // About the epsilon thing see this thread:
993  // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
994  const _RealType __naf =
995  (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
996  // The largest _RealType convertible to _IntType.
997  const _RealType __thr =
998  std::numeric_limits<_IntType>::max() + __naf;
999 
1000  _RealType __cand;
1001  do
1002  __cand = std::ceil(std::log(__urng()) / _M_log_p);
1003  while (__cand >= __thr);
1004 
1005  return result_type(__cand + __naf);
1006  }
1007 
1008  template<typename _IntType, typename _RealType,
1009  typename _CharT, typename _Traits>
1010  std::basic_ostream<_CharT, _Traits>&
1011  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1012  const geometric_distribution<_IntType, _RealType>& __x)
1013  {
1014  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1015  typedef typename __ostream_type::ios_base __ios_base;
1016 
1017  const typename __ios_base::fmtflags __flags = __os.flags();
1018  const _CharT __fill = __os.fill();
1019  const std::streamsize __precision = __os.precision();
1020  __os.flags(__ios_base::scientific | __ios_base::left);
1021  __os.fill(__os.widen(' '));
1022  __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1023 
1024  __os << __x.p();
1025 
1026  __os.flags(__flags);
1027  __os.fill(__fill);
1028  __os.precision(__precision);
1029  return __os;
1030  }
1031 
1032 
1033  template<typename _IntType, typename _RealType>
1034  void
1035  poisson_distribution<_IntType, _RealType>::
1036  _M_initialize()
1037  {
1038 #if _GLIBCXX_USE_C99_MATH_TR1
1039  if (_M_mean >= 12)
1040  {
1041  const _RealType __m = std::floor(_M_mean);
1042  _M_lm_thr = std::log(_M_mean);
1043  _M_lfm = std::tr1::lgamma(__m + 1);
1044  _M_sm = std::sqrt(__m);
1045 
1046  const _RealType __pi_4 = 0.7853981633974483096156608458198757L;
1047  const _RealType __dx = std::sqrt(2 * __m * std::log(32 * __m
1048  / __pi_4));
1049  _M_d = std::tr1::round(std::max(_RealType(6),
1050  std::min(__m, __dx)));
1051  const _RealType __cx = 2 * __m + _M_d;
1052  _M_scx = std::sqrt(__cx / 2);
1053  _M_1cx = 1 / __cx;
1054 
1055  _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
1056  _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2)) / _M_d;
1057  }
1058  else
1059 #endif
1060  _M_lm_thr = std::exp(-_M_mean);
1061  }
1062 
1073  template<typename _IntType, typename _RealType>
1074  template<class _UniformRandomNumberGenerator>
1075  typename poisson_distribution<_IntType, _RealType>::result_type
1076  poisson_distribution<_IntType, _RealType>::
1077  operator()(_UniformRandomNumberGenerator& __urng)
1078  {
1079 #if _GLIBCXX_USE_C99_MATH_TR1
1080  if (_M_mean >= 12)
1081  {
1082  _RealType __x;
1083 
1084  // See comments above...
1085  const _RealType __naf =
1086  (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
1087  const _RealType __thr =
1088  std::numeric_limits<_IntType>::max() + __naf;
1089 
1090  const _RealType __m = std::floor(_M_mean);
1091  // sqrt(pi / 2)
1092  const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
1093  const _RealType __c1 = _M_sm * __spi_2;
1094  const _RealType __c2 = _M_c2b + __c1;
1095  const _RealType __c3 = __c2 + 1;
1096  const _RealType __c4 = __c3 + 1;
1097  // e^(1 / 78)
1098  const _RealType __e178 = 1.0129030479320018583185514777512983L;
1099  const _RealType __c5 = __c4 + __e178;
1100  const _RealType __c = _M_cb + __c5;
1101  const _RealType __2cx = 2 * (2 * __m + _M_d);
1102 
1103  bool __reject = true;
1104  do
1105  {
1106  const _RealType __u = __c * __urng();
1107  const _RealType __e = -std::log(__urng());
1108 
1109  _RealType __w = 0.0;
1110 
1111  if (__u <= __c1)
1112  {
1113  const _RealType __n = _M_nd(__urng);
1114  const _RealType __y = -std::abs(__n) * _M_sm - 1;
1115  __x = std::floor(__y);
1116  __w = -__n * __n / 2;
1117  if (__x < -__m)
1118  continue;
1119  }
1120  else if (__u <= __c2)
1121  {
1122  const _RealType __n = _M_nd(__urng);
1123  const _RealType __y = 1 + std::abs(__n) * _M_scx;
1124  __x = std::ceil(__y);
1125  __w = __y * (2 - __y) * _M_1cx;
1126  if (__x > _M_d)
1127  continue;
1128  }
1129  else if (__u <= __c3)
1130  // NB: This case not in the book, nor in the Errata,
1131  // but should be ok...
1132  __x = -1;
1133  else if (__u <= __c4)
1134  __x = 0;
1135  else if (__u <= __c5)
1136  __x = 1;
1137  else
1138  {
1139  const _RealType __v = -std::log(__urng());
1140  const _RealType __y = _M_d + __v * __2cx / _M_d;
1141  __x = std::ceil(__y);
1142  __w = -_M_d * _M_1cx * (1 + __y / 2);
1143  }
1144 
1145  __reject = (__w - __e - __x * _M_lm_thr
1146  > _M_lfm - std::tr1::lgamma(__x + __m + 1));
1147 
1148  __reject |= __x + __m >= __thr;
1149 
1150  } while (__reject);
1151 
1152  return result_type(__x + __m + __naf);
1153  }
1154  else
1155 #endif
1156  {
1157  _IntType __x = 0;
1158  _RealType __prod = 1.0;
1159 
1160  do
1161  {
1162  __prod *= __urng();
1163  __x += 1;
1164  }
1165  while (__prod > _M_lm_thr);
1166 
1167  return __x - 1;
1168  }
1169  }
1170 
1171  template<typename _IntType, typename _RealType,
1172  typename _CharT, typename _Traits>
1173  std::basic_ostream<_CharT, _Traits>&
1174  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1175  const poisson_distribution<_IntType, _RealType>& __x)
1176  {
1177  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1178  typedef typename __ostream_type::ios_base __ios_base;
1179 
1180  const typename __ios_base::fmtflags __flags = __os.flags();
1181  const _CharT __fill = __os.fill();
1182  const std::streamsize __precision = __os.precision();
1183  const _CharT __space = __os.widen(' ');
1184  __os.flags(__ios_base::scientific | __ios_base::left);
1185  __os.fill(__space);
1186  __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1187 
1188  __os << __x.mean() << __space << __x._M_nd;
1189 
1190  __os.flags(__flags);
1191  __os.fill(__fill);
1192  __os.precision(__precision);
1193  return __os;
1194  }
1195 
1196  template<typename _IntType, typename _RealType,
1197  typename _CharT, typename _Traits>
1198  std::basic_istream<_CharT, _Traits>&
1199  operator>>(std::basic_istream<_CharT, _Traits>& __is,
1200  poisson_distribution<_IntType, _RealType>& __x)
1201  {
1202  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1203  typedef typename __istream_type::ios_base __ios_base;
1204 
1205  const typename __ios_base::fmtflags __flags = __is.flags();
1206  __is.flags(__ios_base::skipws);
1207 
1208  __is >> __x._M_mean >> __x._M_nd;
1209  __x._M_initialize();
1210 
1211  __is.flags(__flags);
1212  return __is;
1213  }
1214 
1215 
1216  template<typename _IntType, typename _RealType>
1217  void
1218  binomial_distribution<_IntType, _RealType>::
1219  _M_initialize()
1220  {
1221  const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1222 
1223  _M_easy = true;
1224 
1225 #if _GLIBCXX_USE_C99_MATH_TR1
1226  if (_M_t * __p12 >= 8)
1227  {
1228  _M_easy = false;
1229  const _RealType __np = std::floor(_M_t * __p12);
1230  const _RealType __pa = __np / _M_t;
1231  const _RealType __1p = 1 - __pa;
1232 
1233  const _RealType __pi_4 = 0.7853981633974483096156608458198757L;
1234  const _RealType __d1x =
1235  std::sqrt(__np * __1p * std::log(32 * __np
1236  / (81 * __pi_4 * __1p)));
1237  _M_d1 = std::tr1::round(std::max(_RealType(1), __d1x));
1238  const _RealType __d2x =
1239  std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
1240  / (__pi_4 * __pa)));
1241  _M_d2 = std::tr1::round(std::max(_RealType(1), __d2x));
1242 
1243  // sqrt(pi / 2)
1244  const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
1245  _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
1246  _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
1247  _M_c = 2 * _M_d1 / __np;
1248  _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
1249  const _RealType __a12 = _M_a1 + _M_s2 * __spi_2;
1250  const _RealType __s1s = _M_s1 * _M_s1;
1251  _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
1252  * 2 * __s1s / _M_d1
1253  * std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
1254  const _RealType __s2s = _M_s2 * _M_s2;
1255  _M_s = (_M_a123 + 2 * __s2s / _M_d2
1256  * std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
1257  _M_lf = (std::tr1::lgamma(__np + 1)
1258  + std::tr1::lgamma(_M_t - __np + 1));
1259  _M_lp1p = std::log(__pa / __1p);
1260 
1261  _M_q = -std::log(1 - (__p12 - __pa) / __1p);
1262  }
1263  else
1264 #endif
1265  _M_q = -std::log(1 - __p12);
1266  }
1267 
1268  template<typename _IntType, typename _RealType>
1269  template<class _UniformRandomNumberGenerator>
1270  typename binomial_distribution<_IntType, _RealType>::result_type
1271  binomial_distribution<_IntType, _RealType>::
1272  _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t)
1273  {
1274  _IntType __x = 0;
1275  _RealType __sum = 0;
1276 
1277  do
1278  {
1279  const _RealType __e = -std::log(__urng());
1280  __sum += __e / (__t - __x);
1281  __x += 1;
1282  }
1283  while (__sum <= _M_q);
1284 
1285  return __x - 1;
1286  }
1287 
1298  template<typename _IntType, typename _RealType>
1299  template<class _UniformRandomNumberGenerator>
1300  typename binomial_distribution<_IntType, _RealType>::result_type
1301  binomial_distribution<_IntType, _RealType>::
1302  operator()(_UniformRandomNumberGenerator& __urng)
1303  {
1304  result_type __ret;
1305  const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
1306 
1307 #if _GLIBCXX_USE_C99_MATH_TR1
1308  if (!_M_easy)
1309  {
1310  _RealType __x;
1311 
1312  // See comments above...
1313  const _RealType __naf =
1314  (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
1315  const _RealType __thr =
1316  std::numeric_limits<_IntType>::max() + __naf;
1317 
1318  const _RealType __np = std::floor(_M_t * __p12);
1319  const _RealType __pa = __np / _M_t;
1320 
1321  // sqrt(pi / 2)
1322  const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
1323  const _RealType __a1 = _M_a1;
1324  const _RealType __a12 = __a1 + _M_s2 * __spi_2;
1325  const _RealType __a123 = _M_a123;
1326  const _RealType __s1s = _M_s1 * _M_s1;
1327  const _RealType __s2s = _M_s2 * _M_s2;
1328 
1329  bool __reject;
1330  do
1331  {
1332  const _RealType __u = _M_s * __urng();
1333 
1334  _RealType __v;
1335 
1336  if (__u <= __a1)
1337  {
1338  const _RealType __n = _M_nd(__urng);
1339  const _RealType __y = _M_s1 * std::abs(__n);
1340  __reject = __y >= _M_d1;
1341  if (!__reject)
1342  {
1343  const _RealType __e = -std::log(__urng());
1344  __x = std::floor(__y);
1345  __v = -__e - __n * __n / 2 + _M_c;
1346  }
1347  }
1348  else if (__u <= __a12)
1349  {
1350  const _RealType __n = _M_nd(__urng);
1351  const _RealType __y = _M_s2 * std::abs(__n);
1352  __reject = __y >= _M_d2;
1353  if (!__reject)
1354  {
1355  const _RealType __e = -std::log(__urng());
1356  __x = std::floor(-__y);
1357  __v = -__e - __n * __n / 2;
1358  }
1359  }
1360  else if (__u <= __a123)
1361  {
1362  const _RealType __e1 = -std::log(__urng());
1363  const _RealType __e2 = -std::log(__urng());
1364 
1365  const _RealType __y = _M_d1 + 2 * __s1s * __e1 / _M_d1;
1366  __x = std::floor(__y);
1367  __v = (-__e2 + _M_d1 * (1 / (_M_t - __np)
1368  -__y / (2 * __s1s)));
1369  __reject = false;
1370  }
1371  else
1372  {
1373  const _RealType __e1 = -std::log(__urng());
1374  const _RealType __e2 = -std::log(__urng());
1375 
1376  const _RealType __y = _M_d2 + 2 * __s2s * __e1 / _M_d2;
1377  __x = std::floor(-__y);
1378  __v = -__e2 - _M_d2 * __y / (2 * __s2s);
1379  __reject = false;
1380  }
1381 
1382  __reject = __reject || __x < -__np || __x > _M_t - __np;
1383  if (!__reject)
1384  {
1385  const _RealType __lfx =
1386  std::tr1::lgamma(__np + __x + 1)
1387  + std::tr1::lgamma(_M_t - (__np + __x) + 1);
1388  __reject = __v > _M_lf - __lfx + __x * _M_lp1p;
1389  }
1390 
1391  __reject |= __x + __np >= __thr;
1392  }
1393  while (__reject);
1394 
1395  __x += __np + __naf;
1396 
1397  const _IntType __z = _M_waiting(__urng, _M_t - _IntType(__x));
1398  __ret = _IntType(__x) + __z;
1399  }
1400  else
1401 #endif
1402  __ret = _M_waiting(__urng, _M_t);
1403 
1404  if (__p12 != _M_p)
1405  __ret = _M_t - __ret;
1406  return __ret;
1407  }
1408 
1409  template<typename _IntType, typename _RealType,
1410  typename _CharT, typename _Traits>
1411  std::basic_ostream<_CharT, _Traits>&
1412  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1413  const binomial_distribution<_IntType, _RealType>& __x)
1414  {
1415  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1416  typedef typename __ostream_type::ios_base __ios_base;
1417 
1418  const typename __ios_base::fmtflags __flags = __os.flags();
1419  const _CharT __fill = __os.fill();
1420  const std::streamsize __precision = __os.precision();
1421  const _CharT __space = __os.widen(' ');
1422  __os.flags(__ios_base::scientific | __ios_base::left);
1423  __os.fill(__space);
1424  __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1425 
1426  __os << __x.t() << __space << __x.p()
1427  << __space << __x._M_nd;
1428 
1429  __os.flags(__flags);
1430  __os.fill(__fill);
1431  __os.precision(__precision);
1432  return __os;
1433  }
1434 
1435  template<typename _IntType, typename _RealType,
1436  typename _CharT, typename _Traits>
1437  std::basic_istream<_CharT, _Traits>&
1438  operator>>(std::basic_istream<_CharT, _Traits>& __is,
1439  binomial_distribution<_IntType, _RealType>& __x)
1440  {
1441  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1442  typedef typename __istream_type::ios_base __ios_base;
1443 
1444  const typename __ios_base::fmtflags __flags = __is.flags();
1445  __is.flags(__ios_base::dec | __ios_base::skipws);
1446 
1447  __is >> __x._M_t >> __x._M_p >> __x._M_nd;
1448  __x._M_initialize();
1449 
1450  __is.flags(__flags);
1451  return __is;
1452  }
1453 
1454 
1455  template<typename _RealType, typename _CharT, typename _Traits>
1456  std::basic_ostream<_CharT, _Traits>&
1457  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1458  const uniform_real<_RealType>& __x)
1459  {
1460  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1461  typedef typename __ostream_type::ios_base __ios_base;
1462 
1463  const typename __ios_base::fmtflags __flags = __os.flags();
1464  const _CharT __fill = __os.fill();
1465  const std::streamsize __precision = __os.precision();
1466  const _CharT __space = __os.widen(' ');
1467  __os.flags(__ios_base::scientific | __ios_base::left);
1468  __os.fill(__space);
1469  __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1470 
1471  __os << __x.min() << __space << __x.max();
1472 
1473  __os.flags(__flags);
1474  __os.fill(__fill);
1475  __os.precision(__precision);
1476  return __os;
1477  }
1478 
1479  template<typename _RealType, typename _CharT, typename _Traits>
1480  std::basic_istream<_CharT, _Traits>&
1481  operator>>(std::basic_istream<_CharT, _Traits>& __is,
1482  uniform_real<_RealType>& __x)
1483  {
1484  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1485  typedef typename __istream_type::ios_base __ios_base;
1486 
1487  const typename __ios_base::fmtflags __flags = __is.flags();
1488  __is.flags(__ios_base::skipws);
1489 
1490  __is >> __x._M_min >> __x._M_max;
1491 
1492  __is.flags(__flags);
1493  return __is;
1494  }
1495 
1496 
1497  template<typename _RealType, typename _CharT, typename _Traits>
1498  std::basic_ostream<_CharT, _Traits>&
1499  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1500  const exponential_distribution<_RealType>& __x)
1501  {
1502  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1503  typedef typename __ostream_type::ios_base __ios_base;
1504 
1505  const typename __ios_base::fmtflags __flags = __os.flags();
1506  const _CharT __fill = __os.fill();
1507  const std::streamsize __precision = __os.precision();
1508  __os.flags(__ios_base::scientific | __ios_base::left);
1509  __os.fill(__os.widen(' '));
1510  __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1511 
1512  __os << __x.lambda();
1513 
1514  __os.flags(__flags);
1515  __os.fill(__fill);
1516  __os.precision(__precision);
1517  return __os;
1518  }
1519 
1520 
1527  template<typename _RealType>
1528  template<class _UniformRandomNumberGenerator>
1529  typename normal_distribution<_RealType>::result_type
1530  normal_distribution<_RealType>::
1531  operator()(_UniformRandomNumberGenerator& __urng)
1532  {
1533  result_type __ret;
1534 
1535  if (_M_saved_available)
1536  {
1537  _M_saved_available = false;
1538  __ret = _M_saved;
1539  }
1540  else
1541  {
1542  result_type __x, __y, __r2;
1543  do
1544  {
1545  __x = result_type(2.0) * __urng() - 1.0;
1546  __y = result_type(2.0) * __urng() - 1.0;
1547  __r2 = __x * __x + __y * __y;
1548  }
1549  while (__r2 > 1.0 || __r2 == 0.0);
1550 
1551  const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
1552  _M_saved = __x * __mult;
1553  _M_saved_available = true;
1554  __ret = __y * __mult;
1555  }
1556 
1557  __ret = __ret * _M_sigma + _M_mean;
1558  return __ret;
1559  }
1560 
1561  template<typename _RealType, typename _CharT, typename _Traits>
1562  std::basic_ostream<_CharT, _Traits>&
1563  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1564  const normal_distribution<_RealType>& __x)
1565  {
1566  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1567  typedef typename __ostream_type::ios_base __ios_base;
1568 
1569  const typename __ios_base::fmtflags __flags = __os.flags();
1570  const _CharT __fill = __os.fill();
1571  const std::streamsize __precision = __os.precision();
1572  const _CharT __space = __os.widen(' ');
1573  __os.flags(__ios_base::scientific | __ios_base::left);
1574  __os.fill(__space);
1575  __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1576 
1577  __os << __x._M_saved_available << __space
1578  << __x.mean() << __space
1579  << __x.sigma();
1580  if (__x._M_saved_available)
1581  __os << __space << __x._M_saved;
1582 
1583  __os.flags(__flags);
1584  __os.fill(__fill);
1585  __os.precision(__precision);
1586  return __os;
1587  }
1588 
1589  template<typename _RealType, typename _CharT, typename _Traits>
1590  std::basic_istream<_CharT, _Traits>&
1591  operator>>(std::basic_istream<_CharT, _Traits>& __is,
1592  normal_distribution<_RealType>& __x)
1593  {
1594  typedef std::basic_istream<_CharT, _Traits> __istream_type;
1595  typedef typename __istream_type::ios_base __ios_base;
1596 
1597  const typename __ios_base::fmtflags __flags = __is.flags();
1598  __is.flags(__ios_base::dec | __ios_base::skipws);
1599 
1600  __is >> __x._M_saved_available >> __x._M_mean
1601  >> __x._M_sigma;
1602  if (__x._M_saved_available)
1603  __is >> __x._M_saved;
1604 
1605  __is.flags(__flags);
1606  return __is;
1607  }
1608 
1609 
1610  template<typename _RealType>
1611  void
1612  gamma_distribution<_RealType>::
1613  _M_initialize()
1614  {
1615  if (_M_alpha >= 1)
1616  _M_l_d = std::sqrt(2 * _M_alpha - 1);
1617  else
1618  _M_l_d = (std::pow(_M_alpha, _M_alpha / (1 - _M_alpha))
1619  * (1 - _M_alpha));
1620  }
1621 
1638  template<typename _RealType>
1639  template<class _UniformRandomNumberGenerator>
1640  typename gamma_distribution<_RealType>::result_type
1641  gamma_distribution<_RealType>::
1642  operator()(_UniformRandomNumberGenerator& __urng)
1643  {
1644  result_type __x;
1645 
1646  bool __reject;
1647  if (_M_alpha >= 1)
1648  {
1649  // alpha - log(4)
1650  const result_type __b = _M_alpha
1651  - result_type(1.3862943611198906188344642429163531L);
1652  const result_type __c = _M_alpha + _M_l_d;
1653  const result_type __1l = 1 / _M_l_d;
1654 
1655  // 1 + log(9 / 2)
1656  const result_type __k = 2.5040773967762740733732583523868748L;
1657 
1658  do
1659  {
1660  const result_type __u = __urng();
1661  const result_type __v = __urng();
1662 
1663  const result_type __y = __1l * std::log(__v / (1 - __v));
1664  __x = _M_alpha * std::exp(__y);
1665 
1666  const result_type __z = __u * __v * __v;
1667  const result_type __r = __b + __c * __y - __x;
1668 
1669  __reject = __r < result_type(4.5) * __z - __k;
1670  if (__reject)
1671  __reject = __r < std::log(__z);
1672  }
1673  while (__reject);
1674  }
1675  else
1676  {
1677  const result_type __c = 1 / _M_alpha;
1678 
1679  do
1680  {
1681  const result_type __z = -std::log(__urng());
1682  const result_type __e = -std::log(__urng());
1683 
1684  __x = std::pow(__z, __c);
1685 
1686  __reject = __z + __e < _M_l_d + __x;
1687  }
1688  while (__reject);
1689  }
1690 
1691  return __x;
1692  }
1693 
1694  template<typename _RealType, typename _CharT, typename _Traits>
1695  std::basic_ostream<_CharT, _Traits>&
1696  operator<<(std::basic_ostream<_CharT, _Traits>& __os,
1697  const gamma_distribution<_RealType>& __x)
1698  {
1699  typedef std::basic_ostream<_CharT, _Traits> __ostream_type;
1700  typedef typename __ostream_type::ios_base __ios_base;
1701 
1702  const typename __ios_base::fmtflags __flags = __os.flags();
1703  const _CharT __fill = __os.fill();
1704  const std::streamsize __precision = __os.precision();
1705  __os.flags(__ios_base::scientific | __ios_base::left);
1706  __os.fill(__os.widen(' '));
1707  __os.precision(__gnu_cxx::__numeric_traits<_RealType>::__max_digits10);
1708 
1709  __os << __x.alpha();
1710 
1711  __os.flags(__flags);
1712  __os.fill(__fill);
1713  __os.precision(__precision);
1714  return __os;
1715  }
1716 
1717 _GLIBCXX_END_NAMESPACE_VERSION
1718 }
1719 }
true
#define true
Definition: stdbool.h:34
__gnu_parallel::min
const _Tp & min(const _Tp &__a, const _Tp &__b)
Equivalent to std::min.
Definition: base.h:144
__gnu_debug::operator>>
std::basic_istream< _CharT, _Traits > & operator>>(std::basic_istream< _CharT, _Traits > &__is, basic_string< _CharT, _Traits, _Allocator > &__str)
Definition: string:1122
__gnu_parallel::max
const _Tp & max(const _Tp &__a, const _Tp &__b)
Equivalent to std::max.
Definition: base.h:150
__gnu_pbds::detail::true_type
std::tr1::integral_constant< int, 1 > true_type
Definition: type_utils.hpp:70
__gnu_pbds::detail::false_type
std::tr1::integral_constant< int, 0 > false_type
Definition: type_utils.hpp:71
Generated on Sat Jul 22 2017 14:41:18 for STLdoc by   doxygen 1.8.6