1*e4b17023SJohn Marino // -*- C++ -*-
2*e4b17023SJohn Marino 
3*e4b17023SJohn Marino // Copyright (C) 2007, 2008, 2009 Free Software Foundation, Inc.
4*e4b17023SJohn Marino //
5*e4b17023SJohn Marino // This file is part of the GNU ISO C++ Library.  This library is free
6*e4b17023SJohn Marino // software; you can redistribute it and/or modify it under the terms
7*e4b17023SJohn Marino // of the GNU General Public License as published by the Free Software
8*e4b17023SJohn Marino // Foundation; either version 3, or (at your option) any later
9*e4b17023SJohn Marino // version.
10*e4b17023SJohn Marino 
11*e4b17023SJohn Marino // This library is distributed in the hope that it will be useful, but
12*e4b17023SJohn Marino // WITHOUT ANY WARRANTY; without even the implied warranty of
13*e4b17023SJohn Marino // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14*e4b17023SJohn Marino // General Public License for more details.
15*e4b17023SJohn Marino 
16*e4b17023SJohn Marino // Under Section 7 of GPL version 3, you are granted additional
17*e4b17023SJohn Marino // permissions described in the GCC Runtime Library Exception, version
18*e4b17023SJohn Marino // 3.1, as published by the Free Software Foundation.
19*e4b17023SJohn Marino 
20*e4b17023SJohn Marino // You should have received a copy of the GNU General Public License and
21*e4b17023SJohn Marino // a copy of the GCC Runtime Library Exception along with this program;
22*e4b17023SJohn Marino // see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
23*e4b17023SJohn Marino // <http://www.gnu.org/licenses/>.
24*e4b17023SJohn Marino 
25*e4b17023SJohn Marino /** @file parallel/random_number.h
26*e4b17023SJohn Marino  *  @brief Random number generator based on the Mersenne twister.
27*e4b17023SJohn Marino  *  This file is a GNU parallel extension to the Standard C++ Library.
28*e4b17023SJohn Marino  */
29*e4b17023SJohn Marino 
30*e4b17023SJohn Marino // Written by Johannes Singler.
31*e4b17023SJohn Marino 
32*e4b17023SJohn Marino #ifndef _GLIBCXX_PARALLEL_RANDOM_NUMBER_H
33*e4b17023SJohn Marino #define _GLIBCXX_PARALLEL_RANDOM_NUMBER_H 1
34*e4b17023SJohn Marino 
35*e4b17023SJohn Marino #include <parallel/types.h>
36*e4b17023SJohn Marino #include <tr1/random>
37*e4b17023SJohn Marino #include <limits>
38*e4b17023SJohn Marino 
39*e4b17023SJohn Marino namespace __gnu_parallel
40*e4b17023SJohn Marino {
41*e4b17023SJohn Marino   /** @brief Random number generator, based on the Mersenne twister. */
42*e4b17023SJohn Marino   class _RandomNumber
43*e4b17023SJohn Marino   {
44*e4b17023SJohn Marino   private:
45*e4b17023SJohn Marino     std::tr1::mt19937   _M_mt;
46*e4b17023SJohn Marino     uint64_t            _M_supremum;
47*e4b17023SJohn Marino     uint64_t            _M_rand_sup;
48*e4b17023SJohn Marino     double              _M_supremum_reciprocal;
49*e4b17023SJohn Marino     double              _M_rand_sup_reciprocal;
50*e4b17023SJohn Marino 
51*e4b17023SJohn Marino     // Assumed to be twice as long as the usual random number.
52*e4b17023SJohn Marino     uint64_t              __cache;
53*e4b17023SJohn Marino 
54*e4b17023SJohn Marino     // Bit results.
55*e4b17023SJohn Marino     int __bits_left;
56*e4b17023SJohn Marino 
57*e4b17023SJohn Marino     static uint32_t
__scale_down(uint64_t __x,uint64_t,double _M_supremum_reciprocal)58*e4b17023SJohn Marino     __scale_down(uint64_t __x,
59*e4b17023SJohn Marino #if _GLIBCXX_SCALE_DOWN_FPU
60*e4b17023SJohn Marino                uint64_t /*_M_supremum*/, double _M_supremum_reciprocal)
61*e4b17023SJohn Marino #else
62*e4b17023SJohn Marino                uint64_t _M_supremum, double /*_M_supremum_reciprocal*/)
63*e4b17023SJohn Marino #endif
64*e4b17023SJohn Marino         {
65*e4b17023SJohn Marino #if _GLIBCXX_SCALE_DOWN_FPU
66*e4b17023SJohn Marino           return uint32_t(__x * _M_supremum_reciprocal);
67*e4b17023SJohn Marino #else
68*e4b17023SJohn Marino           return static_cast<uint32_t>(__x % _M_supremum);
69*e4b17023SJohn Marino #endif
70*e4b17023SJohn Marino         }
71*e4b17023SJohn Marino 
72*e4b17023SJohn Marino   public:
73*e4b17023SJohn Marino     /** @brief Default constructor. Seed with 0. */
_RandomNumber()74*e4b17023SJohn Marino     _RandomNumber()
75*e4b17023SJohn Marino     : _M_mt(0), _M_supremum(0x100000000ULL),
76*e4b17023SJohn Marino       _M_rand_sup(1ULL << std::numeric_limits<uint32_t>::digits),
77*e4b17023SJohn Marino       _M_supremum_reciprocal(double(_M_supremum) / double(_M_rand_sup)),
78*e4b17023SJohn Marino       _M_rand_sup_reciprocal(1.0 / double(_M_rand_sup)),
79*e4b17023SJohn Marino       __cache(0), __bits_left(0) { }
80*e4b17023SJohn Marino 
81*e4b17023SJohn Marino     /** @brief Constructor.
82*e4b17023SJohn Marino      *  @param __seed Random __seed.
83*e4b17023SJohn Marino      *  @param _M_supremum Generate integer random numbers in the
84*e4b17023SJohn Marino      *                  interval @c [0,_M_supremum). */
85*e4b17023SJohn Marino     _RandomNumber(uint32_t __seed, uint64_t _M_supremum = 0x100000000ULL)
_M_mt(__seed)86*e4b17023SJohn Marino     : _M_mt(__seed), _M_supremum(_M_supremum),
87*e4b17023SJohn Marino       _M_rand_sup(1ULL << std::numeric_limits<uint32_t>::digits),
88*e4b17023SJohn Marino       _M_supremum_reciprocal(double(_M_supremum) / double(_M_rand_sup)),
89*e4b17023SJohn Marino       _M_rand_sup_reciprocal(1.0 / double(_M_rand_sup)),
90*e4b17023SJohn Marino       __cache(0), __bits_left(0) { }
91*e4b17023SJohn Marino 
92*e4b17023SJohn Marino     /** @brief Generate unsigned random 32-bit integer. */
93*e4b17023SJohn Marino     uint32_t
operator()94*e4b17023SJohn Marino     operator()()
95*e4b17023SJohn Marino     { return __scale_down(_M_mt(), _M_supremum, _M_supremum_reciprocal); }
96*e4b17023SJohn Marino 
97*e4b17023SJohn Marino     /** @brief Generate unsigned random 32-bit integer in the
98*e4b17023SJohn Marino         interval @c [0,local_supremum). */
99*e4b17023SJohn Marino     uint32_t
operator()100*e4b17023SJohn Marino     operator()(uint64_t local_supremum)
101*e4b17023SJohn Marino     {
102*e4b17023SJohn Marino       return __scale_down(_M_mt(), local_supremum,
103*e4b17023SJohn Marino                         double(local_supremum * _M_rand_sup_reciprocal));
104*e4b17023SJohn Marino     }
105*e4b17023SJohn Marino 
106*e4b17023SJohn Marino     /** @brief Generate a number of random bits, run-time parameter.
107*e4b17023SJohn Marino      *  @param __bits Number of bits to generate. */
108*e4b17023SJohn Marino     unsigned long
__genrand_bits(int __bits)109*e4b17023SJohn Marino     __genrand_bits(int __bits)
110*e4b17023SJohn Marino     {
111*e4b17023SJohn Marino       unsigned long __res = __cache & ((1 << __bits) - 1);
112*e4b17023SJohn Marino       __cache = __cache >> __bits;
113*e4b17023SJohn Marino       __bits_left -= __bits;
114*e4b17023SJohn Marino       if (__bits_left < 32)
115*e4b17023SJohn Marino         {
116*e4b17023SJohn Marino           __cache |= ((uint64_t(_M_mt())) << __bits_left);
117*e4b17023SJohn Marino           __bits_left += 32;
118*e4b17023SJohn Marino         }
119*e4b17023SJohn Marino       return __res;
120*e4b17023SJohn Marino     }
121*e4b17023SJohn Marino };
122*e4b17023SJohn Marino 
123*e4b17023SJohn Marino } // namespace __gnu_parallel
124*e4b17023SJohn Marino 
125*e4b17023SJohn Marino #endif /* _GLIBCXX_PARALLEL_RANDOM_NUMBER_H */
126