---

Rand.h

Go to the documentation of this file.

#ifndef RAND_H
#define RAND_H

#include <cmath>
#include <cstdlib>

const int BIGPRIME = 899999963;

class Rand {
        double c,cd,cm,u[97];
        int i97,j97 ;
        bool outputReady;
        double output;

public:
        Rand(int seed=1) {
                ranmarin(abs(seed%BIGPRIME));
                outputReady=false;
        }

/* Incompatible with Python?
        Rand() {
                ranmarin(1);
                outputReady=false;
        }
*/

        void seed(int seed) {
                ranmarin(abs(seed%BIGPRIME));
                outputReady=false;
        }

        double uniform()
        {
          double uni;
        
          uni=u[i97]-u[j97];
          if (uni<0.0) uni+=1.0;
          u[i97]=uni;
          if (--i97<0) i97=96;
          if (--j97<0) j97=96;
          c-=cd;
          if (c<0.0) c+=cm;
          uni-=c;
          if (uni<0.0) uni+=1.0;
          return(uni);
        }

        // Gives all numbers weighted equally, min and max inclusive
        int discrete(int min, int max)
        {
                int ret= (int)(floor(uniform()*(max+1-min)+min));
//              assert(ret<=max && ret>=min);
                return ret;
        }


        double gaussian(double sd)
        {
                double x,y,r,z;
                if(outputReady)
                {
                        outputReady = false;
                        return output*sd;
                }

                do
                {
                        x=(uniform()*2.0)-1.0;
                        y=(uniform()*2.0)-1.0;
                        r=x*x+y*y;
                } while (r>=1.0);

                z=sqrt(-2.0*log(r)/r);
                output = x*z;
                outputReady=true;
                return y*z*sd;
        }

private:
        void ranmarin(int ijkl)
        {
          int i,ii,j,jj,k,l,m ;
          double s,t ;

      int ij=ijkl/30082;
      int kl=ijkl-30082*ij;

          i=((ij/177)%177)+2 ;
                j=(ij%177)+2 ;
                k=((kl/169)%178)+1 ;
                l=kl%169 ;
                for (ii=0;ii<97;ii++) {
            s=0.0 ;
                        t=0.5 ;
                        for (jj=0;jj<24;jj++) {
              m=(((i*j)%179)*k)%179 ;
                                i=j;
                                j=k;
                                k=m;
                                l=(53*l+1)%169;
                                if (((l*m)%64)>=32) s+=t;
                                t*=0.5;
                        }
                        u[ii]=s;
                }
                c=362436.0/16777216.0;
                cd=7654321.0/16777216.0;
                cm=16777213.0/16777216.0;
                i97=96;
                j97=32;
        }
            
};      

#endif

---