UniformRandom.h

Go to the documentation of this file.

// -*- C++ -*-

#if !defined(__ads_UniformRandom_h__)
#define __ads_UniformRandom_h__

#include "../array/FixedArray.h"

#include <limits>
#include <functional>

#include <cassert>

BEGIN_NAMESPACE_ADS

//-----------------------------------------------------------------------------
// @{


template<typename T>
class UniformRandom {
};


class SubtractiveRNG : 
  public std::unary_function<unsigned int, unsigned int> {
private:

  unsigned int _table[55];
  std::size_t _index1;
  std::size_t _index2;

public:

  unsigned int 
  operator()(unsigned int limit) {
    _index1 = (_index1 + 1) % 55;
    _index2 = (_index2 + 1) % 55;
    _table[_index1] = _table[_index1] - _table[_index2];
    return _table[_index1] % limit;
  }

  void 
  initialize(unsigned int seed) {
    unsigned int k = 1;
    _table[54] = seed;
    size_t i;
    for (i = 0; i < 54; i++) {
      size_t ii = (21 * (i + 1) % 55) - 1;
      _table[ii] = k;
      k = seed - k;
      seed = _table[ii];
    }
    for (int __loop = 0; __loop < 4; __loop++) {
      for (i = 0; i < 55; i++) {
        _table[i] = _table[i] - _table[(1 + i + 30) % 55];
      }
    }
    _index1 = 0;
    _index2 = 31;
  }

  SubtractiveRNG(const unsigned int seed) :
    _index1(),
    _index2() { 
    initialize(seed);
  }

  SubtractiveRNG() :
    _index1(),
    _index2() { 
    initialize(161803398u); 
  }
};



class UniformRandomBase {
private:
  
  //
  // Member data.
  //

  mutable SubtractiveRNG _rng;

protected:

  //
  // Constructors, etc.
  //

  UniformRandomBase() :
    _rng()
  {}

  UniformRandomBase(const unsigned int n) :
    _rng(n)
  {}

  UniformRandomBase(const UniformRandomBase& x) :
    _rng(x._rng)
  {}

  UniformRandomBase&
  operator=( const UniformRandomBase& other) {
    if (&other != this) {
      _rng = other._rng;
    }
    return *this;
  }

  ~UniformRandomBase()
  {}

  //
  // Random number generation.
  //

  unsigned int
  random(const unsigned int n) const {
    return _rng(n);
  }

public:

  void
  initialize(const unsigned int n) {
    SubtractiveRNG tmp(n);
    _rng = tmp;
  }

};



template<typename T>
class UniformRandomInteger :
  public UniformRandomBase {
private:

  //
  // Private types.
  //

  typedef UniformRandomBase Base;

public:

  //
  // Public types.
  //

  typedef T Number;
  typedef void argument_type;
  typedef Number result_type;

private:

  //
  // Member data.
  //

  Number _min;
  Number _extent;

public:

  //--------------------------------------------------------------------------
  // \name Constructors, etc.

  UniformRandomInteger() :
    Base(),
    _min(0),
    _extent(0)
  {}

  UniformRandomInteger(const Number min, const Number max) :
    Base(),
    _min(min),
    _extent(max - min + 1) {
    assert(_extent > 0);
  }

  UniformRandomInteger(const Number min, const Number max,
                       const unsigned int n) :
    Base(n),
    _min(min),
    _extent(max - min + 1) {
    assert(_extent > 0);
  }

  UniformRandomInteger(const UniformRandomInteger& x) :
    Base(x),
    _min(x._min),
    _extent(x._extent)
  {}

  UniformRandomInteger&
  operator=(const UniformRandomInteger& other) {
    if (&other != this) {
      Base::operator=(other);
      _min = other._min;
      _extent = other._extent;
    }
    return *this;
  }
  
  ~UniformRandomInteger()
  {}

  //--------------------------------------------------------------------------
  // \name Random number generation.

  Number
  min() const {
    return _min;
  }
  
  Number
  max() const {
    return _min + _extent - 1;
  }
  
  result_type
  operator()() const {
    return _min + random(_extent);
  }

};



template<typename T>
class UniformRandomReal :
  public UniformRandomBase {
private:

  //
  // Private types.
  //

  typedef UniformRandomBase Base;

public:

  //
  // Public types.
  //

  typedef T Number;
  typedef void argument_type;
  typedef Number result_type;

private:

  //
  // Member data.
  //

  Number _min;
  Number _extent;

public:

  //--------------------------------------------------------------------------
  // \name Constructors, etc.

  UniformRandomReal() :
    Base(),
    _min(0),
    _extent(0)
  {}

  UniformRandomReal(const Number min, const Number max) :
    Base(),
    _min(min),
    _extent(max - min) {
    assert(_extent >= 0);
  }

  UniformRandomReal(const Number min, const Number max,
                     const unsigned int n) :
    Base(n),
    _min(min),
    _extent(max - min) {
    assert(_extent >= 0);
  }

  UniformRandomReal(const UniformRandomReal& x) :
    Base(x),
    _min(x._min),
    _extent(x._extent)
  {}

  UniformRandomReal&
  operator=(const UniformRandomReal& other) {
    if (&other != this) {
      Base::operator=(other);
      _min = other._min;
      _extent = other._extent;
    }
    return *this;
  }
  
  ~UniformRandomReal()
  {}

  //--------------------------------------------------------------------------
  // \name Random number generation.

  Number
  min() const {
    return _min;
  }
  
  Number
  max() const {
    return _min + _extent;
  }
  
  result_type
  operator()() const {
    return _min + _extent * 
      random(std::numeric_limits<unsigned int>::max()) /
      Number(std::numeric_limits<unsigned int>::max() - 1);
  }

};



#define UNIFORM_RANDOM(_T,_Base)\
template<>\
class UniformRandom<_T> :\
  public _Base<_T>\
{\
private:\
  typedef _Base<_T> Base;\
public:\
\
  typedef Base::Number Number;\
\
  typedef Base::argument_type argument_type;\
\
  typedef Base::result_type result_type;\
public:\
\
  UniformRandom() :\
    Base()\
  {}\
\
\
  UniformRandom(const Number min, const Number max) :\
    Base(min, max)\
  {}\
\
\
  UniformRandom(const Number min, const Number max,\
                 const unsigned int n) :\
    Base(min, max, n)\
  {}\
\
\
  UniformRandom(const UniformRandom& x) :\
    Base(x)\
  {}\
\
\
  UniformRandom&\
  operator=(const UniformRandom& x)\
  {\
    if (&x != this) {\
      Base::operator=(x);\
    }\
    return *this;\
  }\
\
\
  ~UniformRandom()\
  {}\
}


//
// The integer types.
//
UNIFORM_RANDOM(char,UniformRandomInteger);
UNIFORM_RANDOM(signed char,UniformRandomInteger);
UNIFORM_RANDOM(unsigned char,UniformRandomInteger);

// Conflicts with unsigned short.
// Uniform random number generator for wchar_t.
//UNIFORM_RANDOM(wchar_t,UniformRandomInteger);

UNIFORM_RANDOM(short,UniformRandomInteger);
UNIFORM_RANDOM(unsigned short,UniformRandomInteger);
UNIFORM_RANDOM(int,UniformRandomInteger);
UNIFORM_RANDOM(unsigned int,UniformRandomInteger);

//
// The real number types.
//
UNIFORM_RANDOM(float,UniformRandomReal);
UNIFORM_RANDOM(double,UniformRandomReal);
UNIFORM_RANDOM(long double,UniformRandomReal);





template<int N, typename T = double>
class UniformRandomPoint {
public:

  //
  // Public types.
  //

  typedef T Number;
  typedef void argument_type;
  typedef ads::FixedArray<N,Number> result_type;

private:

  //
  // Member data.
  //

  // Lower and upper bounds.
  result_type _min, _max;
  // A random number generator.
  UniformRandom<Number> _rng;

public:

  //--------------------------------------------------------------------------
  // \name Constructors, etc.

  UniformRandomPoint() :
    _min(),
    _max(),
    _rng(Number(0), Number(1))
  {}

  UniformRandomPoint(const result_type& min, const result_type& max) :
    _min(min),
    _max(max),
    _rng(Number(0), Number(1))
  {}

  UniformRandomPoint(const result_type& min, const result_type& max,
                     const unsigned int n) :
    _min(min),
    _max(max),
    _rng(Number(0), Number(1), n)
  {}

  UniformRandomPoint(const UniformRandomPoint& x) :
    _min(x._min),
    _max(x._max),
    _rng(x._rng)
  {}

  UniformRandomPoint&
  operator=(const UniformRandomPoint& other) {
    if (&other != this) {
      _min = other._min;
      _max = other._max;
      _rng = other._rng;
    }
    return *this;
  }
  
  ~UniformRandomPoint()
  {}

  //--------------------------------------------------------------------------
  // \name Random number generation.

  const result_type&
  min() const {
    return _min;
  }
  
  const result_type&
  max() const {
    return _max;
  }
  
  const result_type&
  operator()() const {
    static result_type x;
    for (int n = 0; n != N; ++n) {
      x[n] = _min[n] + _rng() * (_max[n] - _min[n]);
    }
    return x;
  }

};


// @}

END_NAMESPACE_ADS

#endif

---