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