// prandom.cp
// Pseudorandom number generation
// Šć1996 Stuart A. Kurtz
// All Rights Reserved

#include "prandom.h"

#include <iostream>
#include <stdlib.h>
#include <math.h>
#define nil NULL


// These routines are based on algorithms described in "Numerical Recipes in C,"
// 2nd edition, by Press, Teukolsky, Vetterling, and Flannery.  A very similar
// discussion can be found in the "Seminumerical Algorithms" volume of Knuth's
// "Art of Programming."
// 
// This code has been rendered in C++ largely for pedagogical purposes.
//
// It's a funny thing: the Numerical Recipes code is far more
// terse, but it involves significant duplication of code.  This
// code exhibits much better modularity, but is a lot longer.
// It seems overwhelming likely that the Numerical Recipes code
// is faster, while this code is easier to understand, easier
// to modify, and easier to experiment with. You take you're
// choice...

const long max_long = 2147483647;                                               // 2^31 - 1
const double scalef = (1.0)/(double(max_long));                 // 1/max_long

// ===================================================================
//      * puniform()
// ===================================================================
// Most of the work is done in the other classes.  This simply takes
// the "default" multiplicative congruental random number generator,
// which returns a result in the range from 1 to max_long, throws
// away the first few elements (because the default seed has very low
// Kolmogorov complexity...), shuffles, and then scales the resulting
// sequence to the range (0.0,1.0).
// 
// Note that we avoid the overhead of using the scale class by simply
// scaling the sequence ourself.  Likewise, the "shift" has been
// folded in to the shuffle class.
// 
// puniform is based on ran1 from Numerical Recipes.

double puniform(long seed)
{
        static shuffle seq(32,7,new mcs);
        
        if (seed != 0) {
                seed = labs(seed);
                seq.set_seed(seed);
                return seed;
        }
        
        while(1) {
                double result = scalef * seq.next();
                if (result < 1.0)       // avoid extreme values...
                        return result;
        }
}

// ===================================================================
//      * pnormal()
// ===================================================================
// This is remarkably similar to the code above, which is the reason
// for all the "overhead" classes.
// 
// pnormal is based on gasdev from Numerical Recipes, p. 289 ff.

double pnormal(long seed)
{
        static gauss_sequence seq(max_long,new shuffle(32,7,new mcs));
        
        if (seed != 0) {
                seed = labs(seed);
                seq.set_seed(seed);
                return seed;
        }
        
        return seq.next();
}


// ===================================================================
//      *       template<class T>
//              class sequence
// ===================================================================

// -------------------------------------------------------------------
//      *       sequence<T>::~sequence()
// -------------------------------------------------------------------

template<class T>
sequence<T>::~sequence()
{}

// ===================================================================
//      *       class mcs
// ===================================================================

// -------------------------------------------------------------------
//      *       mcs::mcs(long _seed, long _a, long _m)
// -------------------------------------------------------------------

const long good_a = 16807;
const long good_m = max_long;

mcs::mcs(long _seed, long _a, long _m) :
        state(labs(_seed)),
        a(_a),
        m(_m)
{
        if (a == 0)
                a = good_a;
        if (m == 0)
                m = good_m;
        if (state == 0)
                state = 1;              // not a great choice, perhaps...
        
        q = m/a;
        r = m%a;
}

// -------------------------------------------------------------------
//      *       mcs::~mcs()
// -------------------------------------------------------------------
// mcs inherits from a class with a virtual destructor, which it must
// override.

mcs::~mcs()
{}

// -------------------------------------------------------------------
//      *       mcs::next()
// -------------------------------------------------------------------

long mcs::next()
{
        const long q = 127773;
        const long r = 2836;
        
        // The following code amounts to
        //   state = (a * state) % m,
        // computed using Schrage's algorithm --
        // which avoids the possibility overflow,
        // cf. Numerical Recipes p. 278.

        long k = state/q;
        state = a * (state - k * q) - k * r;
        if (state < 0)
                state += m;

        return state;
}

// -------------------------------------------------------------------
//      *       void mcs::set_seed(long seed)
// -------------------------------------------------------------------

void mcs::set_seed(long seed)
{
        state = seed;
}

// ===================================================================
//      *       template<class T, class U>
//              class filter_shuffle
// ===================================================================

// -------------------------------------------------------------------
//      *       filtered_sequence<T,U>::filtered_sequence(sequence<U> *_seq)
// -------------------------------------------------------------------

template<class T, class U>
filtered_sequence<T,U>::filtered_sequence(sequence<U> *_seq) :
        seq(_seq)
{
        if (seq == nil) {
                cerr << "Attempt to filter nil sequence." << endl;
                exit(-1);
        }
}

// -------------------------------------------------------------------
//      *       filtered_sequence<T,U>::~filtered_sequence()
// -------------------------------------------------------------------

template<class T, class U>
filtered_sequence<T,U>::~filtered_sequence()
{
        delete seq;
}

// -------------------------------------------------------------------
//      *       filtered_sequence<T,U>::void set_seed(long seed)
// -------------------------------------------------------------------

template<class T, class U>
void filtered_sequence<T,U>::set_seed(long seed)
{
        seq->set_seed(seed);
        reset();
}

// -------------------------------------------------------------------
//      *       void filtered_sequence<T,U>::reset()
// -------------------------------------------------------------------

template<class T, class U>
void filtered_sequence<T,U>::reset()
{}

// ===================================================================
//      *       template<class T>
//              shift<T>
// ===================================================================

// -------------------------------------------------------------------
//      *       shift<T>::shift(long delay,sequence<T> *seq)
// -------------------------------------------------------------------

template<class T>
shift<T>::shift(long delay,sequence<T> *seq) :
        filtered_sequence<T,T>(seq),
        delay(delay)
{
        reset();
}

// -------------------------------------------------------------------
//      *       shift<T>::~shift()
// -------------------------------------------------------------------
// shift<T> inherits from a class with a virtual destructor, which it
// must override.


template<class T>
shift<T>::~shift()
{}

// -------------------------------------------------------------------
//      *       T shift<T>::next()
// -------------------------------------------------------------------

template<class T>
T shift<T>::next()
{
        return seq->next();
}

// -------------------------------------------------------------------
//      *       void shift<T>::reset()
// -------------------------------------------------------------------

template<class T>
void shift<T>::reset()
{
        // throw away the first few elements of seq...
        for (int i = 0; i < delay; ++i)
                seq->next();
}

// ===================================================================
//      *       class shuffle
// ===================================================================

// -------------------------------------------------------------------
//      *       shuffle::shuffle(long tabsize,sequence<long> *seq)
// -------------------------------------------------------------------

shuffle::shuffle(long tabsize, long throw_away, sequence<long> *seq) :
        filtered_sequence<long,long>(seq),
        throw_away(throw_away),
        tabsize(tabsize),
        ndiv(1+(max_long-1)/tabsize)
{
        if (tabsize <= 0) {
                cerr << "Fatal error, nonpositive array size: " << tabsize << endl;
                exit(-1);
        }

        array = new long[tabsize];
        
        if (array == nil) {
                cerr << "Fatal error, could not allocate shuffle array." << endl;
                exit(-1);
        }
        
        reset();
}

// -------------------------------------------------------------------
//      *       shuffle::~shuffle
// -------------------------------------------------------------------

shuffle::~shuffle()
{
        delete array;
}

// -------------------------------------------------------------------
//      *       long shuffle::next()
// -------------------------------------------------------------------
// Here's where the real work takes place.  We use the last element
// to select an array entry.  Note that we use
//
//              long index = last / ndiv;
// 
// rather than the apparently more natural
// 
//              long index = last % tabsize;
// 
// The rationale is that last % tabsize depends only on the low-order
// bits of last, which often have an unexpectly short period.  The
// division effectively relies on the high-order bits, which do not
// usually have short periods.
// 
// Once the index is computed, we replace the index'th element of the
// array with a new value, returning the old value.

long shuffle::next()
{
        long index = last / ndiv;
        last = array[index];
        array[index] = seq->next();
        
        return last;
}

// -------------------------------------------------------------------
//      *       long shuffle::reset()
// -------------------------------------------------------------------
// This simply loads both last and the shuffle array with the first
// few elements from the sequence.

void shuffle::reset()
{
        for (long i = 0; i < throw_away; ++i)
                seq->next();

        last = seq->next();
        
        for (long i = 0; i < tabsize; ++i)
                array[i] = seq->next();
}

// ===================================================================
//      *       template<class T, class U>
//              scale<T,U>
// ===================================================================

// -------------------------------------------------------------------
//      *       scale<T,U>::scale(T scale, sequence<U> *seq)
// -------------------------------------------------------------------

template<class T, class U>
scale<T,U>::scale(T scale, sequence<U> *seq) :
        filtered_sequence<T,U>(seq),
        s(scale)
{}

// -------------------------------------------------------------------
//      *       scale<T,U>::~scale()
// -------------------------------------------------------------------
// scale<T,U> inherits from a class with a virtual destructor, which it
// must override.

template<class T, class U>
scale<T,U>::~scale()
{}

// -------------------------------------------------------------------
//      *       T scale<T,U>::next()
// -------------------------------------------------------------------
// Here's the real work. The rest of this class -- in all its
// baroque glory -- exists to support this one little line.

template<class T, class U>
T scale<T,U>::next()
{
        return s * T(seq->next());
}

// ===================================================================
//      *       class gauss_sequence
// ===================================================================

// -------------------------------------------------------------------
//      *       gauss_sequence::gauss_sequence(sequence<double> *seq)
// -------------------------------------------------------------------

gauss_sequence::gauss_sequence(long max, sequence<long> *seq) :
        filtered_sequence<double,long>(seq),
        scalef(1.0/double(max))
{
        reset();
}

// -------------------------------------------------------------------
//      *       gauss_sequence::~gauss_sequence()
// -------------------------------------------------------------------
// gauss_sequence inherits from a class with a virtual destructor, which it
// must override.

gauss_sequence::~gauss_sequence()
{}

// -------------------------------------------------------------------
//      *       double gauss_sequence::next()
// -------------------------------------------------------------------
// This is the standard approach.  We assume that gauss_sequence is
// filtering a uniform sequence of longs in the range 0 .. max.  This
// sequence is scaled a double in the range (-1.0,1.0).  We generate
// reals in pairs -- essentially generating uniformly distributed points
// in the box (-1.0,1.0) x (-1.0,1.0).  This is done until we find a
// point in the unit ball, which is then subjected to the "usual"
// transformation, cf., Numerical Recipes in C++, p. 289 ff.
//
// This technique generates points in pairs.  One point is returned, the
// other cached for the next call.

double gauss_sequence::next()
{
        if (valid_cache) {
                valid_cache = false;
                return cached_value;
        }
        
        static double v1, v2, ss;
        
        do {
                v1 = 2.0 * scalef * double(seq->next()) - 1.0;
                v2 = 2.0 * scalef * double(seq->next()) - 1.0;
                ss = v1*v1 + v2*v2;
        } while (ss >= 1.0 || ss == 0.0);
        
        double scale = sqrt(-2.0 * log(ss)/ss);
        
        cached_value = v2 * scale;
        valid_cache = true;
        
        return v1 * scale;
}

// -------------------------------------------------------------------
//      *       void gauss_sequence::reset()
// -------------------------------------------------------------------
// Resetting the sequence requires that we start over -- thus any
// cached value must be invalidated.

void gauss_sequence::reset()
{
        valid_cache = false;
}