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