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