egs_rndm.h

Go to the documentation of this file.

/*
###############################################################################
#
#  EGSnrc egs++ random number generators headers
#  Copyright (C) 2015 National Research Council Canada
#
#  This file is part of EGSnrc.
#
#  EGSnrc is free software: you can redistribute it and/or modify it under
#  the terms of the GNU Affero General Public License as published by the
#  Free Software Foundation, either version 3 of the License, or (at your
#  option) any later version.
#
#  EGSnrc 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 Affero General Public License for
#  more details.
#
#  You should have received a copy of the GNU Affero General Public License
#  along with EGSnrc. If not, see <http://www.gnu.org/licenses/>.
#
###############################################################################
#
#  Author:          Iwan Kawrakow, 2005
#
#  Contributors:
#
###############################################################################
*/


#ifndef EGS_RANDOM_GENERATOR_
#define EGS_RANDOM_GENERATOR_

#include "egs_libconfig.h"
#include "egs_math.h"
#include "egs_functions.h"

class EGS_Input;

class EGS_EXPORT EGS_RandomGenerator {

public:

    EGS_RandomGenerator(int n=128);

    EGS_RandomGenerator(const EGS_RandomGenerator &r) : np(0) {
        copyBaseState(r);
    };

    virtual ~EGS_RandomGenerator() {
        delete [] rarray;
    };

    inline EGS_Float getUniform() {
        if (ip >= np) {
            fillArray(np,rarray);
            ip = 0;
        }
        return rarray[ip++];
    };

    EGS_I64 numbersGenerated() const {
        return count;
    };

    EGS_I64 numbersUsed() const {
        return ip<np ? count - np + ip : count;
    };

    inline void getAzimuth(EGS_Float &cphi, EGS_Float &sphi) {
#ifndef FAST_SINCOS
        register EGS_Float xphi,xphi2,yphi,yphi2,rhophi;
        do {
            xphi = 2*getUniform() - 1;
            xphi2 = xphi*xphi;
            yphi = getUniform();
            yphi2 = yphi*yphi;
            rhophi = xphi2 + yphi2;
        }
        while (rhophi > 1);
        cphi = (xphi2 - yphi2)/rhophi;
        sphi = 2*xphi*yphi/rhophi;
#else
        EGS_Float phi = 2*M_PI*getUniform();
        cphi = cos(phi);
        sphi = sin(phi);
#endif
    };

    inline EGS_Float getGaussian() {
        if (have_x) {
            have_x = false;
            return the_x;
        }
        EGS_Float r = sqrt(-2*log(1-getUniform()));
        EGS_Float cphi, sphi;
        getAzimuth(cphi,sphi);
        have_x = true;
        the_x = r*sphi;
        return r*cphi;
    };

    static EGS_RandomGenerator *createRNG(EGS_Input *inp, int sequence=0);

    static EGS_RandomGenerator *defaultRNG(int sequence=0);

    virtual void fillArray(int n, EGS_Float *array) = 0;


    bool storeState(ostream &data);
    bool setState(istream &data);
    bool addState(istream &data);
    void resetCounter() {
        count = 0;
    };

    virtual EGS_RandomGenerator *getCopy() = 0;

    virtual void setState(EGS_RandomGenerator *r) = 0;

    virtual void saveState() = 0;

    virtual void resetState() = 0;

    virtual int  rngSize() const = 0;

    virtual void describeRNG() const {};

protected:

    EGS_I64   count;  
    int       np;     
    int       ip;     
    EGS_Float *rarray;

    virtual bool storePrivateState(ostream &data) =  0;

    virtual bool setPrivateState(istream &data) =  0;

    void copyBaseState(const EGS_RandomGenerator &r);

    void allocate(int n);

    int  baseSize() const {
        return 2*sizeof(EGS_I32) + sizeof(EGS_I64) + sizeof(bool) +
               sizeof(EGS_Float) + np*sizeof(EGS_Float);
    };


private:

    bool      have_x; 
    EGS_Float the_x;  

};

#endif