random_number.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2007, 2008 Free Software Foundation, Inc.
00004 //
00005 // This file is part of the GNU ISO C++ Library.  This library is free
00006 // software; you can redistribute it and/or modify it under the terms
00007 // of the GNU General Public License as published by the Free Software
00008 // Foundation; either version 2, or (at your option) any later
00009 // version.
00010 
00011 // This library is distributed in the hope that it will be useful, but
00012 // WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014 // General Public License for more details.
00015 
00016 // You should have received a copy of the GNU General Public License
00017 // along with this library; see the file COPYING.  If not, write to
00018 // the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
00019 // MA 02111-1307, USA.
00020 
00021 // As a special exception, you may use this file as part of a free
00022 // software library without restriction.  Specifically, if other files
00023 // instantiate templates or use macros or inline functions from this
00024 // file, or you compile this file and link it with other files to
00025 // produce an executable, this file does not by itself cause the
00026 // resulting executable to be covered by the GNU General Public
00027 // License.  This exception does not however invalidate any other
00028 // reasons why the executable file might be covered by the GNU General
00029 // Public License.
00030 
00031 /** @file parallel/random_number.h
00032  *  @brief Random number generator based on the Mersenne twister.
00033  *  This file is a GNU parallel extension to the Standard C++ Library.
00034  */
00035 
00036 // Written by Johannes Singler.
00037 
00038 #ifndef _GLIBCXX_PARALLEL_RANDOM_NUMBER_H
00039 #define _GLIBCXX_PARALLEL_RANDOM_NUMBER_H 1
00040 
00041 #include <parallel/types.h>
00042 #include <tr1/random>
00043 
00044 namespace __gnu_parallel
00045 {
00046   /** @brief Random number generator, based on the Mersenne twister. */
00047   class random_number
00048   {
00049   private:
00050     std::tr1::mt19937   mt;
00051     uint64      supremum;
00052     uint64      RAND_SUP;
00053     double      supremum_reciprocal;
00054     double      RAND_SUP_REC;
00055 
00056     // Assumed to be twice as long as the usual random number.
00057     uint64      cache;  
00058 
00059     // Bit results.
00060     int bits_left;
00061     
00062     static uint32
00063     scale_down(uint64 x,
00064 #if _GLIBCXX_SCALE_DOWN_FPU
00065            uint64 /*supremum*/, double supremum_reciprocal)
00066 #else
00067                uint64 supremum, double /*supremum_reciprocal*/)
00068 #endif
00069     {
00070 #if _GLIBCXX_SCALE_DOWN_FPU
00071       return uint32(x * supremum_reciprocal);
00072 #else
00073       return static_cast<uint32>(x % supremum);
00074 #endif
00075     }
00076 
00077   public:
00078     /** @brief Default constructor. Seed with 0. */
00079     random_number()
00080     : mt(0), supremum(0x100000000ULL),
00081       RAND_SUP(1ULL << (sizeof(uint32) * 8)),
00082       supremum_reciprocal(double(supremum) / double(RAND_SUP)),
00083       RAND_SUP_REC(1.0 / double(RAND_SUP)),
00084       cache(0), bits_left(0) { }
00085 
00086     /** @brief Constructor.
00087      *  @param seed Random seed.
00088      *  @param supremum Generate integer random numbers in the
00089      *                  interval @c [0,supremum). */
00090     random_number(uint32 seed, uint64 supremum = 0x100000000ULL)
00091     : mt(seed), supremum(supremum),
00092       RAND_SUP(1ULL << (sizeof(uint32) * 8)),
00093       supremum_reciprocal(double(supremum) / double(RAND_SUP)),
00094       RAND_SUP_REC(1.0 / double(RAND_SUP)),
00095       cache(0), bits_left(0) { }
00096 
00097     /** @brief Generate unsigned random 32-bit integer. */
00098     uint32
00099     operator()()
00100     { return scale_down(mt(), supremum, supremum_reciprocal); }
00101 
00102     /** @brief Generate unsigned random 32-bit integer in the
00103     interval @c [0,local_supremum). */
00104     uint32
00105     operator()(uint64 local_supremum)
00106     {
00107       return scale_down(mt(), local_supremum,
00108             double(local_supremum * RAND_SUP_REC));
00109     }
00110 
00111     /** @brief Generate a number of random bits, run-time parameter.
00112      *  @param bits Number of bits to generate. */
00113     unsigned long
00114     genrand_bits(int bits)
00115     {
00116       unsigned long res = cache & ((1 << bits) - 1);
00117       cache = cache >> bits;
00118       bits_left -= bits;
00119       if (bits_left < 32)
00120     {
00121       cache |= ((uint64(mt())) << bits_left);
00122       bits_left += 32;
00123     }
00124       return res;
00125     }
00126 };
00127 
00128 } // namespace __gnu_parallel
00129 
00130 #endif

Generated on Fri Jan 23 20:12:16 2009 for libstdc++ by  doxygen 1.5.6