Random Numbers class wrapper (based on Julienne Walker's implementation)

I'm wrapping Prelude's implementation of jsw_rand inside a class. So far this is what I have:

CRand.h

PHP Code:

---

#ifndef CLASS_RANDOM_GENERATOR_H
#define CLASS_RANDOM_GENERATOR_H

/*
  Random number Class

    - Wrapped around Julienne Walker's implementation of Mersenne Twister's found at
      http://eternallyconfuzzled.com/brain.html, as of July, 18, 2006
    - Slightly altered to accomodate the class implementation
*/

class CRand {
private:
    static unsigned int time_seed();

    static const int N = 624;
    static const int M = 397;
    static const unsigned long A  = 0x9908b0dfUL;
    static const unsigned long U  = 0x80000000UL;
    static const unsigned long L  = 0x7fffffffUL;

public:
    explicit CRand(unsigned long, unsigned long);

    // TODO: void reseed(unsigned long = time_seed());
    // TODO: unsigned long rethrow();
private:

    void seed(unsigned long = time_seed());
   // TODO: unsigned long generate();

    unsigned long low_;
    unsigned long high_;
    unsigned long res_; //holds random

    // internal state (per object)
    unsigned long x[N];
    int next;
};

#endif // CLASS_RANDOM_GENERATOR_H 


---

CRand.cpp

PHP Code:

---

#include <climits>
#include <ctime>
#include "crand.h"

CRand::CRand(unsigned long low, unsigned long high): low_(low), high_(high) {
    if(low_ > high_) {
        low_ = high;
        high_ = low;
    }

    if( low_ != high_ ) {
        seed();
        //TODO: res_ = generate();
    } else res_ = low_;
}

unsigned int CRand::time_seed() {
  std::time_t now = std::time(0);
  unsigned char *p = (unsigned char *)&now;
  unsigned seed = 0;
  std::size_t i;

  for (i = 0; i < sizeof now; i++)
    seed = seed * ( UCHAR_MAX + 2U ) + p[i];

  return seed;
}

/* Initialize internal state */
void CRand::seed(unsigned long s) {
  int i;

  x[0] = s & 0xffffffffUL;

  for ( i = 1; i < N; i++ ) {
    x[i] = ( 1812433253UL
      * ( x[i - 1] ^ ( x[i - 1] >> 30 ) ) + i );
    x[i] &= 0xffffffffUL;
  }
} 


---

I'm about to implement the actual jsw_rand() algorithm into the class. However, I'm afraid I couldn't fully understand the discussion about distribution found at http://eternallyconfuzzled.com/articles/rand.html.

My question is, will it be ok, after using Prelude's algorithm jsw_rand(), to use the modulos operator to bring the resulting long into the range defined by my data members low_ and high_? Or should I instead use the uniform deviation she discusses?

Basically... since I'm having an hard time interpreting the code inside jsw_rand() (and because math is not one of my best attributes) I'm at loss as to which it already provides a distribution safeguard(?). My guess is it doesn't and I should either call it repeatedly or provide a uniform deviation algorithm as she discusses.

Quote:

---

However, I'm afraid I couldn't fully understand the discussion about distribution found at http://eternallyconfuzzled.com/articles/rand.html.

---

Let's say you have a random number generator that produced integers in the range [0,2] with equal probability for each.
Suppose you wanted either 0 or 1 for your random number.

Using the modulo method, it is easy to see that you would end up with 0 twice as many times as 1, in the long run.
For example, if you produced 3000 random numbers, you would get around 1000 occurrences of 0, 1000 occurrences of 1, and 1000 occurrences of 2. But 2 is congruent to 0 under modulo 2, hence your final result would give you around 2000 occurrences of 0 and 1000 occurrences of 1.

Quote:

---

My question is, will it be ok, after using Prelude's algorithm jsw_rand(), to use the modulos operator to bring the resulting long into the range defined by my data members low_ and high_? Or should I instead use the uniform deviation she discusses?

---

Use her suggestions.

Ah! Excellent. I understand it now. And will do as suggested.

>My question is, will it be ok, after using Prelude's algorithm jsw_rand(), to use
>the modulos operator to bring the resulting long into the range defined by my
>data members low_ and high_? Or should I instead use the uniform deviation she discusses?
You can use the remainder operator to shrink the range. There's really nothing wrong with it if you don't mind a little munging of the distribution. If you really do need a uniform distribution (in a statistics program, for example), it's better to go with the methods that preserve it.

>I'm at loss as to which it already provides a distribution safeguard(?).
No, it doesn't. There's really no way to go about that since it's impossible to determine what range the user wants and the algorithm has to be tailored to the range. So the algorithm sets the range at something appropriately large and assumes the user can adjust it accordingly.

>My guess is it doesn't and I should either call it repeatedly or provide a uniform deviation algorithm as she discusses.
If I were wrapping it in a class, I'd provide several public interface functions that give the user a range of options. That way it's more than just a wrapper, it's a convenient interface. And because the time seed could be useful elsewhere, I'd break it out into a separate class of static seeding functions and make sure that the default seed is always predictable (for debugging). Something like this:

Code:

---

#include <vector>

namespace jsw {
  class seeder {
  public:
    typedef unsigned int seed_t;

    static seed_t sys_time();
    // Other seeding methods
  };

  class random {
  public:
    typedef unsigned long rand_t;
    typedef double frand_t;
    typedef std::vector<rand_t> range_t;
    typedef std::vector<frand_t> frange_t;

    static const rand_t rand_max = <max range of algo>;

    random ( seeder::seed_t init = 1 );
    void seed ( seeder::seed_t init = 1 );
    frand_t next();
    rand_t next ( rand_t low = 0, rand_t high = rand_max );
    frange_t next_range ( frange_t::size_type n );
    range_t next_range ( frange_t::size_type n,
      rand_t low = 0, rand_t high = rand_max );
  };
}

---

Thanks Prelude. I will seek to your suggestion for a second version of this class. Need to study yourcode further before I do though. :) Still struggle with C++ syntax here and there.

Meanwhile, I finished what can be considered a first version of the CRand class. A second one will come with Prelude's suggestions. I would like though a comment on the way I'm implementing the uniform deviation solution. You can find it under CRand::generate() towards the end of the function.

Code:

---

/*
  Random number Class

    - Wrapped around Julienne Walker's implementation of Mersenne Twister's found at
      http://eternallyconfuzzled.com/brain.html, as of July, 18, 2005
    - Slightly altered to accomodate the class implementation

    Constructor:
      - a pair of unsigned longs establishing an inclusive range
      - an enum establishing the type of distribution pretended

    Interface:
      Each instance of the CRand class represents an object that can be used to generate
      random integrals within the range provided by the constructor. When initialized, the
      object will seed itself and generate the first random. This number can be read
      through last_throw().
      Any further random number generation is done through rethrow() which both stores and
      returns the result. The method to get this value later on is always through
      last_throw(). Mind you, CRand only stores the last number generated.
      It is possible to reseed the object through seed_again(). Seeding defaults to
      time_seed(). So it can be used without any argument.
      The range can be changed and checked through the overloaded low() and high()
      functions. CRange automatically fixes the range if a low() is provided that is
      higher than high() or the opposite, a high() is providade that is lower than the
      current low() by swapping the values. This range checking is also done in the
      constructor.

    Usage Example:
        #include <cstdlib>
        #include <iostream>
        #include "crand.h"

        int main()
        {
            CRand test(0,3); // creates an object test. Random numbers are between 0 and 3
                            // and distribution type is NONE.

            std::cout << test.last_throw(); // returns the number generated when test was
                                            // initialized

            // run 20 times. Each time get a new random number between 0 and 3 using a
            // uniform deviation based distribution.
            for(int i = 0; i != 20; ++i) {
                std::cout << i + 1 << ": " ;
                std::cout << test.rethrow(CRand::Uniform_Deviation);
                std::cout << std::endl;
            }

            int value = test.last_throw() // get last random number generated.

            return EXIT_SUCCESS;
        }
*/

---

A reduced and less complicated version.

- Removed that silly idea of storing the last random generated and thus simplified the class usage considerably.

- Replaced high() and low() with max() and min(). Didn't know what I was thinking...

- Sizeof still reports 2508 bytes. But that is the price to pay for a Mersenne Twister, I guess.

- The more I work on this, the less I'm convinced a class exclusively meant to produce random numbers using only one method is a useful thing. :)

I do need someone's opinion on the way I'm implementing the uniform distribution method on CRand::generate(). I have this idea I'm totally missing the point. Using a Mersenne Twister and then use the result only to seed a rand() call makes little sense.

PHP Code:

---

unsigned long CRand::generate(Distribution dist) {
    double uniform_deviate;

    /* ... Mersenne Twister implemented here ...*/

    // last stages of Mersenne Twister
    y ^= (y >> 11);
    y ^= (y << 7) & 0x9d2c5680UL;
    y ^= (y << 15) & 0xefc60000UL;
    y ^= (y >> 18);

    /* set within range with chosen distribution */
    switch(dist) {
        case None:
            y = min_ + y % ( max_ + 1 - min_ );
            break;
        case Uniform_Deviation:
            srand(static_cast<int>(y));
            uniform_deviate = rand() * ( 1.0 / ( RAND_MAX + 1.0 ) );
            y = static_cast<unsigned long>( min_ + uniform_deviate * (max_ + 1 - min_) );
            break;
    }

    return y;
} 


---

Quote:

---

Originally Posted by Mario F.

- The more I work on this, the less I'm convinced a class exclusively meant to produce random numbers using only one method is a useful thing. :)

---

Which is why there's Boost Random, which has become part of TR1 and thus is an officially recommended extension to the C++ standard library.