randomgenerator.h

Go to the documentation of this file.

/********************************************************************************************
**    iLand - an individual based forest landscape and disturbance model
**    http://iland-model.org
**    Copyright (C) 2009-  Werner Rammer, Rupert Seidl
**
**    This program is free software: you can redistribute it and/or modify
**    it under the terms of the GNU General Public License as published by
**    the Free Software Foundation, either version 3 of the License, or
**    (at your option) any later version.
**
**    This program is distributed in the hope that it will be useful,
**    but WITHOUT ANY WARRANTY; without even the implied warranty of
**    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
**    GNU General Public License for more details.
**
**    You should have received a copy of the GNU General Public License
**    along with this program.  If not, see <http://www.gnu.org/licenses/>.
********************************************************************************************/
 
#ifndef RANDOMGENERATOR_H
#define RANDOMGENERATOR_H
#include <cstdlib>
#include <math.h>
#include <time.h>
 
#define RANDOMGENERATORSIZE 500000
#define RANDOMGENERATORROTATIONS 0
// a new set of numbers is generated for every 5*500000 = 2.500.000 numbers
class RandomGenerator
{
public:
    enum ERandomGenerators { ergMersenneTwister, ergWellRNG512, ergXORShift96, ergFastRandom };
    RandomGenerator() { seed(0); setGeneratorType(ergMersenneTwister); }
    static void setGeneratorType(const ERandomGenerators gen) { mGeneratorType = gen; mRotationCount=RANDOMGENERATORROTATIONS+1;mIndex=0;mRefillCounter=0; }
    static void debugState(int &rIndex, int &rGeneration, int &rRefillCount) { rIndex = mIndex; rGeneration = mRotationCount; rRefillCount = mRefillCounter; }
    static int debugNRandomNumbers() { return mIndex + RANDOMGENERATORSIZE*mRotationCount + (RANDOMGENERATORROTATIONS+1)*RANDOMGENERATORSIZE*mRefillCounter; }
    static void checkGenerator() { if (mRotationCount>RANDOMGENERATORROTATIONS) { refill();  } }
    static void setup(const ERandomGenerators gen, const unsigned oneSeed) { setGeneratorType(gen); seed(oneSeed); checkGenerator(); }
    static void seed(const unsigned oneSeed);
    static inline double rand() { return next() * (1.0/4294967295.0); }
    static inline double rand(const double max_value) { return max_value * rand(); }
    static inline unsigned long randInt(){ return next(); }
    static inline unsigned long randInt(const int max_value) { return max_value>0?randInt()%max_value:0; }
    static inline double randNorm( const double mean, const double stddev );
 
private:
    static inline unsigned long next() { ++mIndex; if (mIndex>RANDOMGENERATORSIZE) { mRotationCount++; mIndex=0; checkGenerator(); }  return mBuffer[mIndex]; }
    static unsigned int mBuffer[RANDOMGENERATORSIZE+5];
    static int mIndex;
    static int mRotationCount;
    static int mRefillCounter;
    static ERandomGenerators mGeneratorType;
    static void refill();
};
 
 
 
inline double nrandom(const double& p1, const double& p2)
{
    return p1 + RandomGenerator::rand(p2-p1);
    //return p1 + (p2-p1)*(rand()/double(RAND_MAX));
}
inline double drandom()
{
    return RandomGenerator::rand();
    //return rand()/double(RAND_MAX);
}
inline int irandom(int from, int to)
{
    return from + RandomGenerator::randInt(to-from);
    //return from +  rand()%(to-from);
}
 
 
 
 
inline double RandomGenerator::randNorm( const double mean, const double stddev )
{
        // Return a real number from a normal (Gaussian) distribution with given
        // mean and standard deviation by polar form of Box-Muller transformation
        double x, y, r;
        do
        {
                x = 2.0 * rand() - 1.0;
                y = 2.0 * rand() - 1.0;
                r = x * x + y * y;
        }
        while ( r >= 1.0 || r == 0.0 );
        double s = sqrt( -2.0 * log(r) / r );
        return mean + x * s * stddev;
}
 
 
 
 
 
 
 
 
#endif // RANDOMGENERATOR_H