libstdc++-v3/include/parallel/random_number.h

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