OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-stable/] [gcc-4.5.1/] [libstdc++-v3/] [include/] [tr1/] [random.tcc] - Blame information for rev 847

Go to most recent revision | Details | Compare with Previous | View Log

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

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.