OSDN Git Service

2006-06-06 Paolo Carlini <pcarlini@suse.de>
[pf3gnuchains/gcc-fork.git] / libstdc++-v3 / include / tr1 / random.tcc
1 // random number generation (out of line) -*- C++ -*-
2
3 // Copyright (C) 2006 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 2, 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 // You should have received a copy of the GNU General Public License along
17 // with this library; see the file COPYING.  If not, write to the Free
18 // Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301,
19 // USA.
20
21 // As a special exception, you may use this file as part of a free software
22 // library without restriction.  Specifically, if other files instantiate
23 // templates or use macros or inline functions from this file, or you compile
24 // this file and link it with other files to produce an executable, this
25 // file does not by itself cause the resulting executable to be covered by
26 // the GNU General Public License.  This exception does not however
27 // invalidate any other reasons why the executable file might be covered by
28 // the GNU General Public License.
29
30 #include <limits>
31
32 namespace std
33 {
34 _GLIBCXX_BEGIN_NAMESPACE(tr1)
35
36   /*
37    * Implementation-space details.
38    */
39   namespace _Private
40   {
41     // General case for x = (ax + c) mod m -- use Schrage's algorithm to avoid
42     // integer overflow.
43     //
44     // Because a and c are compile-time integral constants the compiler kindly
45     // elides any unreachable paths.
46     //
47     // Preconditions:  a > 0, m > 0.
48     //
49     template<typename _Tp, _Tp __a, _Tp __c, _Tp __m, bool>
50       struct _Mod
51       {
52         static _Tp
53         __calc(_Tp __x)
54         {
55           if (__a == 1)
56             __x %= __m;
57           else
58             {
59               static const _Tp __q = __m / __a;
60               static const _Tp __r = __m % __a;
61               
62               _Tp __t1 = __a * (__x % __q);
63               _Tp __t2 = __r * (__x / __q);
64               if (__t1 >= __t2)
65                 __x = __t1 - __t2;
66               else
67                 __x = __m - __t2 + __t1;
68             }
69
70           if (__c != 0)
71             {
72               const _Tp __d = __m - __x;
73               if (__d > __c)
74                 __x += __c;
75               else
76                 __x = __c - __d;
77             }
78           return __x;
79         }
80       };
81
82     // Special case for m == 0 -- use unsigned integer overflow as modulo
83     // operator.
84     template<typename _Tp, _Tp __a, _Tp __c, _Tp __m>
85       struct _Mod<_Tp, __a, __c, __m, true>
86       {
87         static _Tp
88         __calc(_Tp __x)
89         { return __a * __x + __c; }
90       };
91
92     // Dispatch based on modulus value to prevent divide-by-zero compile-time
93     // errors when m == 0.
94     template<typename _Tp, _Tp __a, _Tp __c, _Tp __m>
95       inline _Tp
96       __mod(_Tp __x)
97       { return _Mod<_Tp, __a, __c, __m, __m == 0>::__calc(__x); }
98
99     // Like the above, for a == 1, c == 0, in terms of w.
100     template<typename _Tp, _Tp __w, bool>
101       struct _Mod_w
102       {
103         static _Tp
104         __calc(_Tp __x)
105         { return __x % (_Tp(1) << __w); }
106       };
107
108     template<typename _Tp, _Tp __w>
109       struct _Mod_w<_Tp, __w, true>
110       {
111         static _Tp
112         __calc(_Tp __x)
113         { return __x; }
114       };
115
116     template<typename _Tp, _Tp __w>
117       inline _Tp
118       __mod_w(_Tp __x)
119       { return _Mod_w<_Tp, __w,
120                       __w == std::numeric_limits<_Tp>::digits>::__calc(__x); }
121
122     // Selector to return the maximum value possible that will fit in 
123     // @p __w bits of @p _Tp.
124     template<typename _Tp, _Tp __w, bool>
125       struct _Max_w
126       {
127         static _Tp
128         __value()
129         { return (_Tp(1) << __w) - 1; }
130       };
131
132     template<typename _Tp, _Tp __w>
133       struct _Max_w<_Tp, __w, true>
134       {
135         static _Tp
136         __value()
137         { return std::numeric_limits<_Tp>::max(); }
138       };
139
140   } // namespace _Private
141
142
143   /**
144    * Constructs the LCR engine with integral seed @p __x0.
145    */
146   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
147     linear_congruential<_UIntType, __a, __c, __m>::
148     linear_congruential(unsigned long __x0)
149     { this->seed(__x0); }
150
151   /**
152    * Constructs the LCR engine with seed generated from @p __g.
153    */
154   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
155     template<class _Gen>
156       linear_congruential<_UIntType, __a, __c, __m>::
157       linear_congruential(_Gen& __g)
158       { this->seed(__g); }
159
160   /**
161    * Seeds the LCR with integral value @p __x0, adjusted so that the 
162    * ring identity is never a member of the convergence set.
163    */
164   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
165     void
166     linear_congruential<_UIntType, __a, __c, __m>::
167     seed(unsigned long __x0)
168     {
169       if ((_Private::__mod<_UIntType, 1, 0, __m>(__c) == 0)
170           && (_Private::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
171         _M_x = _Private::__mod<_UIntType, 1, 0, __m>(1);
172       else
173         _M_x = _Private::__mod<_UIntType, 1, 0, __m>(__x0);
174     }
175
176   /**
177    * Seeds the LCR engine with a value generated by @p __g.
178    */
179   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
180     template<class _Gen>
181       void
182       linear_congruential<_UIntType, __a, __c, __m>::
183       seed(_Gen& __g, false_type)
184       {
185         _UIntType __x0 = __g();
186         if ((_Private::__mod<_UIntType, 1, 0, __m>(__c) == 0)
187             && (_Private::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
188           _M_x = _Private::__mod<_UIntType, 1, 0, __m>(1);
189         else
190           _M_x = _Private::__mod<_UIntType, 1, 0, __m>(__x0);
191       }
192
193   /**
194    * Returns a value that is less than or equal to all values potentially
195    * returned by operator(). The return value of this function does not
196    * change during the lifetime of the object..
197    *
198    * The minumum depends on the @p __c parameter: if it is zero, the
199    * minimum generated must be > 0, otherwise 0 is allowed.
200    */
201   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
202     typename linear_congruential<_UIntType, __a, __c, __m>::result_type
203     linear_congruential<_UIntType, __a, __c, __m>::
204     min() const
205     { return (_Private::__mod<_UIntType, 1, 0, __m>(__c) == 0) ? 1 : 0; }
206
207   /**
208    * Gets the maximum possible value of the generated range.
209    *
210    * For a linear congruential generator, the maximum is always @p __m - 1.
211    */
212   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
213     typename linear_congruential<_UIntType, __a, __c, __m>::result_type
214     linear_congruential<_UIntType, __a, __c, __m>::
215     max() const
216     { return (__m == 0) ? std::numeric_limits<_UIntType>::max() : (__m - 1); }
217
218   /**
219    * Gets the next generated value in sequence.
220    */
221   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
222     typename linear_congruential<_UIntType, __a, __c, __m>::result_type
223     linear_congruential<_UIntType, __a, __c, __m>::
224     operator()()
225     {
226       _M_x = _Private::__mod<_UIntType, __a, __c, __m>(_M_x);
227       return _M_x;
228     }
229
230
231   template<class _UIntType, int __w, int __n, int __m, int __r,
232            _UIntType __a, int __u, int __s,
233            _UIntType __b, int __t, _UIntType __c, int __l>
234     void
235     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
236                      __b, __t, __c, __l>::
237     seed(unsigned long __value)
238     {
239       _M_x[0] = _Private::__mod_w<_UIntType, __w>(__value);
240
241       for (int __i = 1; __i < state_size; ++__i)
242         {
243           _UIntType __x = _M_x[__i - 1];
244           __x ^= __x >> (__w - 2);
245           __x *= 1812433253ul;
246           __x += __i;
247           _M_x[__i] = _Private::__mod_w<_UIntType, __w>(__x);     
248         }
249       _M_p = state_size;
250     }
251
252   template<class _UIntType, int __w, int __n, int __m, int __r,
253            _UIntType __a, int __u, int __s,
254            _UIntType __b, int __t, _UIntType __c, int __l>
255     template<class _Gen>
256       void
257       mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
258                        __b, __t, __c, __l>::
259       seed(_Gen& __gen, false_type)
260       {
261         for (int __i = 0; __i < state_size; ++__i)
262           _M_x[__i] = _Private::__mod_w<_UIntType, __w>(__gen());
263         _M_p = state_size;
264       }
265
266   template<class _UIntType, int __w, int __n, int __m, int __r,
267            _UIntType __a, int __u, int __s,
268            _UIntType __b, int __t, _UIntType __c, int __l>
269     typename
270     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
271                      __b, __t, __c, __l>::result_type
272     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
273                      __b, __t, __c, __l>::
274     max() const
275     {
276       using _Private::_Max_w;
277       using std::numeric_limits;
278       return _Max_w<_UIntType, __w,
279                     __w == numeric_limits<_UIntType>::digits>::__value();
280     }
281
282   template<class _UIntType, int __w, int __n, int __m, int __r,
283            _UIntType __a, int __u, int __s,
284            _UIntType __b, int __t, _UIntType __c, int __l>
285     typename
286     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
287                      __b, __t, __c, __l>::result_type
288     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
289                      __b, __t, __c, __l>::
290     operator()()
291     {
292       // Reload the vector - cost is O(n) amortized over n calls.
293       if (_M_p >= state_size)
294         {
295           const _UIntType __upper_mask = (~_UIntType()) << __r;
296           const _UIntType __lower_mask = ~__upper_mask;
297
298           for (int __k = 0; __k < (__n - __m); ++__k)
299             {
300               _UIntType __y = ((_M_x[__k] & __upper_mask)
301                                |(_M_x[__k + 1] & __lower_mask));
302               _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
303                            ^ ((__y & 0x01) ? __a : 0));
304             }
305
306           for (int __k = (__n - __m); __k < (__n - 1); ++__k)
307             {
308               _UIntType __y = ((_M_x[__k] & __upper_mask)
309                                | (_M_x[__k + 1] & __lower_mask));
310               _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
311                            ^ ((__y & 0x01) ? __a : 0));
312             }
313
314           _M_p = 0;
315         }
316
317       // Calculate o(x(i)).
318       result_type __z = _M_x[_M_p++];
319       __z ^= (__z >> __u);
320       __z ^= (__z << __s) & __b;
321       __z ^= (__z << __t) & __c;
322       __z ^= (__z >> __l);
323
324       return __z;
325     }
326
327
328   template<typename _IntType, _IntType __m, int __s, int __r>
329     void
330     subtract_with_carry<_IntType, __m, __s, __r>::
331     seed(unsigned long __value)
332     {
333       std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563>
334         __lcg(__value);
335
336       for (int __i = 0; __i < long_lag; ++__i)
337         _M_x[__i] = _Private::__mod<_IntType, 1, 0, modulus>(__lcg());
338
339       _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
340       _M_p = 0;
341     }
342
343   //
344   // This implementation differs from the tr1 spec because the tr1 spec refused
345   // to make any sense to me:  the exponent of the factor in the spec goes from
346   // 1 to (n-1), but it would only make sense to me if it went from 0 to (n-1).
347   //
348   // This algorithm is still problematic because it can overflow left right and
349   // center.
350   //
351   template<typename _IntType, _IntType __m, int __s, int __r>
352     template<class _Gen>
353     void
354     subtract_with_carry<_IntType, __m, __s, __r>::
355     seed(_Gen& __gen, false_type)
356     {
357       const int __n = (std::numeric_limits<_IntType>::digits + 31) / 32;
358       for (int __i = 0; __i < long_lag; ++__i)
359         {
360           _M_x[__i] = 0;
361           unsigned long __factor = 1;
362           for (int __j = 0; __j < __n; ++__j)
363             {
364               _M_x[__i] += __gen() * __factor;
365               __factor *= 0x80000000;
366             }
367           _M_x[__i] = _Private::__mod<_IntType, 1, 0, modulus>(_M_x[__i]);
368         }
369       _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
370       _M_p = 0;
371     }
372
373   template<typename _IntType, _IntType __m, int __s, int __r>
374     typename subtract_with_carry<_IntType, __m, __s, __r>::result_type
375     subtract_with_carry<_IntType, __m, __s, __r>::
376     operator()()
377     {
378       // Derive short lag index from current index.
379       int __ps = _M_p - short_lag;
380       if (__ps < 0)
381         __ps += long_lag;
382
383       // Calculate new x(i) without overflow or division.
384       _IntType __xi;
385       if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
386         {
387           __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
388           _M_carry = 0;
389         }
390       else
391         {
392           __xi = modulus - _M_x[_M_p] - _M_carry + _M_x[__ps];
393           _M_carry = 1;
394         }
395       _M_x[_M_p++] = __xi;
396
397       // Adjust current index to loop around in ring buffer.
398       if (_M_p >= long_lag)
399         _M_p = 0;
400
401       return __xi;
402     }
403
404
405   template<class _UniformRandomNumberGenerator, int __p, int __r>
406     typename discard_block<_UniformRandomNumberGenerator,
407                            __p, __r>::result_type
408     discard_block<_UniformRandomNumberGenerator, __p, __r>::
409     operator()()
410     {
411       if (_M_n >= used_block)
412         {
413           while (_M_n < block_size)
414             {
415               _M_b();
416               ++_M_n;
417             }
418           _M_n = 0;
419         }
420       ++_M_n;
421       return _M_b();
422     }
423
424 _GLIBCXX_END_NAMESPACE
425 }