pseudoRandom.cpp

Go to the documentation of this file.

#include <dlrCommon/exception.h>
#include <dlrLinearAlgebra/clapack.h>
#include <dlrPortability/timeUtilities.h>
#include <dlrRandom/pseudoRandom.h>

namespace dlr {

  namespace random {

    // The default constructor initializes the random number generator
    // with a seed derived from the system clock.
    PseudoRandom::
    PseudoRandom()
    {

      // Note that time resolution is not critical here.  The only
      // impact of low resolution is that two PseudorRandom instances
      // created at nearly the same time may wind up on the same
      // sequence.  If you're worried about that, you can use the
      // getCurrentSeed() member function to make sure it hasn't
      // happened.
      double currentTime = dlr::portability::getCurrentTime();
      Int32 seconds = static_cast<Int32>(currentTime);
      Int32 uSeconds = static_cast<Int32>(1000000.0 * (currentTime - seconds));
    
      // All seeds must be in the range [0, 4096).
      Int32 seed0 = static_cast<Int32>(seconds & 0x00000fff);
      Int32 seed1 = static_cast<Int32>((seconds & 0x00fff000) >> 12);
      Int32 seed2 = static_cast<Int32>((uSeconds & 0x00000fff));
      Int32 seed3 = static_cast<Int32>((uSeconds & 0x00fff000) >> 12);

      // seed3 must be odd.
      if(seed3 % 2 == 0) {
        seed3 -= 1;
      }

      // All seeds must be >= zero.
      if(seed3 < 0) {
        seed3 = -seed3;
      }

      // try {
      this->setLapackSeed(seed0, seed1, seed2, seed3);
      // } catch(const ValueException&) {
      // This won't happen because we've chosen the seed values wisely.
      // }
    }


    // This constructor sets the seed of the random number generator.
    PseudoRandom::
    PseudoRandom(Int64 seed)
    {
      this->setCurrentSeed(seed);
    }

  
    // This member function returns a double drawn from a Gaussian
    // distribution with the specified mean and standard deviation.
    double
    PseudoRandom::
    gaussian(double mu, double sigma)
    {
      double returnValue = this->normal();
      return (returnValue * sigma) + mu;
    }


    // This member function returns the current state of the random
    // number generator.
    Int64
    PseudoRandom::
    getCurrentSeed()
    {
      Int64 seed = 0;
      seed += static_cast<Int64>(m_seed[0]) << 35;
      seed += static_cast<Int64>(m_seed[1]) << 23;
      seed += static_cast<Int64>(m_seed[2]) << 11;
      seed += static_cast<Int64>((m_seed[3] & 0x0ffe) >> 1);
      return seed;
    }
  

    // This member function returns a double drawn from a Gaussian
    // distribution with equal to 0.0 and and standard deviation equal to
    // 1.0.
    double
    PseudoRandom::
    normal()
    {
      Int32 idist = 3;   // Tells lapack we want a normal distribution.
      double returnValue;   // The nuber we'll return.
      Int32 size = 1;    // Only one value needed.

      // Call the lapack routine
      dlarnv_(&idist, m_seed, &size, &returnValue);

      return returnValue;
    }


    //This member function sets the seed for the random number
    // generator.
    void
    PseudoRandom::
    setCurrentSeed(Int64 seed)
    {
      // All seeds must be in the range [0,4096).
      Int32 seed0 = static_cast<Int32>((seed & 0x00007ff800000000LL) >> 35);
      Int32 seed1 = static_cast<Int32>((seed & 0x00000007ff800000LL) >> 23);
      Int32 seed2 = static_cast<Int32>((seed & 0x00000000007ff800LL) >> 11);
    
      // Seed3 must be odd, 
      Int32 seed3 = static_cast<Int32>(((seed & 0x00000000000007ffLL) << 1) + 1);

      // try {
      this->setLapackSeed(seed0, seed1, seed2, seed3);
      // } catch(const ValueException&) {
      // This won't happen because we've chosen the seed values wisely.
      // }
    }

  
    // This member function returns a double, x, drawn from a uniform
    // distribution.
    double
    PseudoRandom::
    uniform(double lowerBound, double upperBound)
    {
      Int32 idist = 1;   // Tells lapack we want uniform [0,1)
      double returnValue;   // The number we'll return
      Int32 size = 1;    // Only one value needed.

      // Call the lapack routine
      dlarnv_(&idist, m_seed, &size, &returnValue);

      double range = upperBound - lowerBound;
      returnValue = returnValue * range + lowerBound;
      return returnValue;
    }


    // This member function returns a integer, x, drawn from a uniform
    // distribution.
    int
    PseudoRandom::
    uniformInt(int lowerBound, int upperBound)
    {
      return static_cast<int>(
        this->uniform(static_cast<double>(lowerBound),
                      static_cast<double>(upperBound)));
    }

    void
    PseudoRandom::
    setLapackSeed(Int32 seed0, Int32 seed1, Int32 seed2, Int32 seed3)
    {
      if(seed3 % 2 == 0) {
        DLR_THROW(ValueException,
                  "PseudoRandom::setLapackSeed(Int32, Int32, Int32, Int32)",
                  "Seed 3 must be odd.");
      }
      if(seed0 >= 4096 || seed1 >= 4096 || seed2 >= 4096 || seed3 >= 4096
         || seed0 < 0 || seed1 < 0 || seed2 < 0 || seed3 < 0 ) {
        DLR_THROW(ValueException,
                  "PseudoRandom::PseudoRandom(Int32, Int32, Int32, Int32)",
                  "All seeds must be >= 0 and < 4096.");
      }
      m_seed[0] = seed0;
      m_seed[1] = seed1;
      m_seed[2] = seed2;
      m_seed[3] = seed3;
    }

    
  } // namespace random
    
} // namespace dlr

---