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