src/gromacs/random/uniformintdistribution.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2015,2016, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35
  36 /*! \file
  37  * \brief The uniform integer distribution
  38  *
  39  * Portable version of the uniform integer that generates the same sequence
  40  * on all platforms. Since stdlibc++ and libc++ provide different sequences
  41  * we prefer this one so unit tests produce the same values on all platforms.
  42  *
  43  * \author Erik Lindahl <erik.lindahl@gmail.com>
  44  * \inpublicapi
  45  * \ingroup module_random
  46  */
  47
  48 #ifndef GMX_RANDOM_UNIFORMINTDISTRIBUTION_H
  49 #define GMX_RANDOM_UNIFORMINTDISTRIBUTION_H
  50
  51 #include <limits>
  52
  53 #include "gromacs/math/functions.h"
  54 #include "gromacs/utility/basedefinitions.h"
  55 #include "gromacs/utility/classhelpers.h"
  56 #include "gromacs/utility/gmxassert.h"
  57
  58 namespace gmx
  59 {
  60
  61 /*! \brief Uniform integer distribution
  62  *
  63  *  The C++ standard library does provide this distribution, but even
  64  *  though they all sample from the correct distribution different standard
  65  *  library implementations appear to return different sequences of numbers
  66  *  for the same random number generator. To make it easier to use GROMACS
  67  *  unit tests that depend on random numbers we have our own implementation.
  68  *
  69  * \tparam IntType Integer type, int by default.
  70  */
  71 template<class IntType = int>
  72 class UniformIntDistribution
  73 {
  74     public:
  75         /*! \brief Type of values returned */
  76         typedef IntType result_type;
  77
  78         /*! \brief Uniform int distribution parameters */
  79         class param_type
  80         {
  81             /*! \brief Lower end of range (inclusive) */
  82             result_type  a_;
  83             /*! \brief Upper end of range (inclusive) */
  84             result_type  b_;
  85
  86             public:
  87                 /*! \brief Reference back to the distribution class */
  88                 typedef UniformIntDistribution distribution_type;
  89
  90                 /*! \brief Construct parameter block
  91                  *
  92                  * \param a   Lower end of range (inclusive)
  93                  * \param b   Upper end of range (inclusive)
  94                  */
  95                 explicit param_type(result_type a = 0, result_type b = std::numeric_limits<result_type>::max())
  96                     : a_(a), b_(b)
  97                 {
  98                     GMX_RELEASE_ASSERT(a <= b, "The uniform integer distribution requires a<=b");
  99                 }
 100
 101                 /*! \brief Return lower range */
 102                 result_type a() const { return a_; }
 103                 /*! \brief Return upper range */
 104                 result_type b() const { return b_; }
 105
 106                 /*! \brief True if two parameter sets will return the same uniform int distribution.
 107                  *
 108                  * \param x  Instance to compare with.
 109                  */
 110                 bool
 111                 operator==(const param_type &x) const
 112                 {
 113                     // rangeBits is a function of a & b, so it does not have to be tested
 114                     return a_ == x.a_ && b_ == x.b_;
 115                 }
 116
 117                 /*! \brief True if two parameter sets will return different uniform int distributions
 118                  *
 119                  * \param x  Instance to compare with.
 120                  */
 121                 bool
 122                 operator!=(const param_type &x) const { return !operator==(x); }
 123         };
 124
 125     public:
 126
 127         /*! \brief Construct new distribution with given integer parameters.
 128          *
 129          * \param a   Lower end of range (inclusive)
 130          * \param b   Upper end of range (inclusive)
 131          */
 132         explicit UniformIntDistribution(result_type a = 0, result_type b = std::numeric_limits<result_type>::max())
 133             : param_(param_type(a, b)), savedRandomBits_(0), savedRandomBitsLeft_(0) {}
 134
 135         /*! \brief Construct new distribution from parameter class
 136          *
 137          * \param param  Parameter class as defined inside gmx::UniformIntDistribution.
 138          */
 139         explicit UniformIntDistribution(const param_type &param)
 140             : param_(param), savedRandomBits_(0), savedRandomBitsLeft_(0) {}
 141
 142         /*! \brief Flush all internal saved values  */
 143         void
 144         reset() { savedRandomBitsLeft_ = 0; }
 145
 146         /*! \brief Return values from uniform int distribution with internal parameters
 147          *
 148          * \tparam Rng  Uniform random engine class
 149          *
 150          * \param  g    Random engine
 151          */
 152         template<class Rng>
 153         result_type
 154         operator()(Rng &g) { return (*this)(g, param_); }
 155
 156         /*! \brief Return value from uniform int distribution with given parameters
 157          *
 158          * \tparam Rng   Uniform random engine class
 159          *
 160          * \param  g     Random engine
 161          * \param  param Parameters to use
 162          */
 163         template<class Rng>
 164         result_type
 165         operator()(Rng &g, const param_type &param)
 166         {
 167             static_assert(sizeof(typename Rng::result_type) >= sizeof(gmx_uint32_t),
 168                           "The random engine result_type should be 32 or 64 bits");
 169
 170             result_type  range = param.b() - param.a();
 171             unsigned int rangeBits;
 172             result_type  result;
 173
 174             if (range == 0)
 175             {
 176                 return param.a();
 177             }
 178             else if (range == std::numeric_limits<result_type>::max())
 179             {
 180                 rangeBits = std::numeric_limits<result_type>::digits; // Use all bits in type
 181             }
 182             else
 183             {
 184                 if (sizeof(result_type) == sizeof(gmx_uint32_t))
 185                 {
 186                     rangeBits = log2I(static_cast<gmx_uint32_t>(range));
 187                 }
 188                 else
 189                 {
 190                     rangeBits = log2I(range);
 191                 }
 192                 rangeBits += ((range >> rangeBits) > 0);
 193             }
 194
 195             do
 196             {
 197                 if (savedRandomBitsLeft_ < rangeBits)
 198                 {
 199                     savedRandomBits_     = static_cast<gmx_uint64_t>(g());
 200                     savedRandomBitsLeft_ = std::numeric_limits<typename Rng::result_type>::digits;
 201
 202                     if (sizeof(typename Rng::result_type) == sizeof(gmx_uint32_t))
 203                     {
 204                         savedRandomBits_    <<= std::numeric_limits<gmx_uint32_t>::digits;
 205                         savedRandomBits_     |= g();
 206                         savedRandomBitsLeft_ += std::numeric_limits<gmx_uint32_t>::digits;
 207                     }
 208                 }
 209                 result                   = savedRandomBits_;
 210                 savedRandomBits_       >>= rangeBits;
 211                 result                   = result - (savedRandomBits_ << rangeBits);
 212                 savedRandomBitsLeft_    -= rangeBits;
 213             }
 214             while (result > range);
 215
 216             return result + param.a();
 217         }
 218
 219         /*! \brief Return the lower range uniform int distribution */
 220         result_type
 221         a() const { return param_.a(); }
 222
 223         /*! \brief Return the upper range of the uniform int distribution */
 224         result_type
 225         b() const { return param_.b(); }
 226
 227         /*! \brief Return the full parameter class of the uniform int distribution */
 228         param_type param() const { return param_; }
 229
 230         /*! \brief Smallest value that can be returned from uniform int distribution */
 231         result_type
 232         min() const { return a(); }
 233
 234         /*! \brief Largest value that can be returned from uniform int distribution */
 235         result_type
 236         max() const { return b(); }
 237
 238         /*! \brief True if two uniform int distributions will produce the same values.
 239          *
 240          * \param  x     Instance to compare with.
 241          */
 242         bool
 243         operator==(const UniformIntDistribution &x) const
 244         { return param_ == x.param_; }
 245
 246         /*! \brief True if two uniform int distributions will produce different values.
 247          *
 248          * \param  x     Instance to compare with.
 249          */
 250         bool
 251         operator!=(const UniformIntDistribution &x) const
 252         { return !operator==(x); }
 253
 254     private:
 255         /*! \brief Internal value for parameters, can be overridden at generation time. */
 256         param_type      param_;
 257         /*! \brief Saved output from random engine, shifted tableBits right each time */
 258         gmx_uint64_t    savedRandomBits_;
 259         /*! \brief Number of valid bits remaining i savedRandomBits_ */
 260         unsigned int    savedRandomBitsLeft_;
 261
 262         GMX_DISALLOW_COPY_AND_ASSIGN(UniformIntDistribution);
 263 };
 264
 265 }      // namespace gmx
 266
 267 #endif // GMX_RANDOM_UNIFORMINTDISTRIBUTION_H