From 74c11a1da369a371925fbf4a7550837468bcda6e Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Mon, 10 Feb 2014 16:01:39 +0000 Subject: [PATCH] ranlxd random generator --- configure.ac | 10 + examples/Makefile.am | 7 +- examples/exRand.cpp | 42 +++ latan/Io.cpp | 11 + latan/Io.hpp | 6 +- latan/IoAsciiLexer.lpp | 1 + latan/IoAsciiParser.ypp | 56 +++- latan/IoObject.hpp | 7 +- latan/Makefile.am | 2 + latan/RandGen.cpp | 696 ++++++++++++++++++++++++++++++++++++++++ latan/RandGen.hpp | 101 ++++++ latan/includes.hpp | 2 + 12 files changed, 931 insertions(+), 10 deletions(-) create mode 100644 examples/exRand.cpp create mode 100644 latan/RandGen.cpp create mode 100644 latan/RandGen.hpp diff --git a/configure.ac b/configure.ac index 7431e52..8b5da33 100644 --- a/configure.ac +++ b/configure.ac @@ -25,6 +25,16 @@ m4_ifdef([AM_PROG_AR],[AM_PROG_AR]) LT_INIT +# Option to enable SSE random generator +AC_ARG_ENABLE([SSE], + [AS_HELP_STRING([--enable-SSE], + [compiles SSE version of ranlux random generator])], + [AC_DEFINE([HAVE_SSE], + [1], + [Define to 1 if your CPU support SSE instructions.])], + [] +) + # Get compilers informations AX_COMPILER_VENDOR AC_DEFINE_UNQUOTED([C_COMP_VENDOR],["$ax_cv_c_compiler_vendor"], diff --git a/examples/Makefile.am b/examples/Makefile.am index 32daf83..c34bd24 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -18,7 +18,8 @@ noinst_PROGRAMS = \ exCompiledDoubleFunction\ exMat \ exMathInterpreter \ - exPlot + exPlot \ + exRand exCompiledDoubleFunction_SOURCES = exCompiledDoubleFunction.cpp exCompiledDoubleFunction_CFLAGS = -g -O2 @@ -36,4 +37,8 @@ exPlot_SOURCES = exPlot.cpp exPlot_CFLAGS = -g -O2 exPlot_LDFLAGS = -L../latan/.libs -llatan +exRand_SOURCES = exRand.cpp +exRand_CFLAGS = -g -O2 +exRand_LDFLAGS = -L../latan/.libs -llatan + ACLOCAL_AMFLAGS = -I .buildutils/m4 diff --git a/examples/exRand.cpp b/examples/exRand.cpp new file mode 100644 index 0000000..2229aaa --- /dev/null +++ b/examples/exRand.cpp @@ -0,0 +1,42 @@ +#include +#include +#include + +using namespace std; +using namespace Latan; + +const int seqLength = 25; +const int saveStep = 9; +const string stateFileName = "exRand.seed"; + +int main(void) +{ + RandGen::State state; + RandGen gen[2]; + AsciiFile stateFile(stateFileName, File::Mode::write|File::Mode::read); + + cout << "- GENERATOR STATE I/O TESTS" << endl; + cout << "-- generating a " << seqLength << " steps random sequence..." + << endl; + for (int i = 0; i < seqLength; ++i) + { + if (i == saveStep) + { + state = gen[0].getState(); + stateFile.save(state, "exRand"); + cout << "generator state after step " << saveStep - 1 + << " saved in '" << stateFileName << "'" << endl; + } + cout << "step " << i << "\t: " << gen[0].uniform() <("exRand")); + cout << "-- generating a " << seqLength << " steps random sequence..." + << endl; + for (int i = 0; i < seqLength; ++i) + { + cout << "step " << i << "\t: " << gen[1].uniform() < #include #include +#include #include BEGIN_NAMESPACE @@ -63,7 +64,8 @@ public: unsigned int getMode(void) const; template const IoT & read(const std::string &name); - virtual void save(const DMat &m, const std::string &name) = 0; + virtual void save(const DMat &m, const std::string &name) = 0; + virtual void save(const RandGen::State &state, const std::string &name) = 0; // tests virtual bool isOpen(void) const = 0; // IO @@ -125,6 +127,7 @@ public: // public members std::stack dMatBuf; std::stack doubleBuf; + std::stack intBuf; private: // allocation/deallocation functions defined in IoAsciiLexer.lpp virtual void initScanner(void); @@ -138,6 +141,7 @@ public: virtual ~AsciiFile(void); // access virtual void save(const DMat &m, const std::string &name); + virtual void save(const RandGen::State &state, const std::string &name); // tests virtual bool isOpen(void) const; // IO diff --git a/latan/IoAsciiLexer.lpp b/latan/IoAsciiLexer.lpp index 8c87598..aeaa3de 100644 --- a/latan/IoAsciiLexer.lpp +++ b/latan/IoAsciiLexer.lpp @@ -75,6 +75,7 @@ BLANK [ \t] latan_begin {BEGIN(TYPE); RETTOK(OPEN);} latan_end {BEGIN(TYPE); RETTOK(CLOSE);} mat {BEGIN(INITIAL); RETTOK(MAT);} +rg_state {BEGIN(INITIAL); RETTOK(RG_STATE);} <*>\n {yylloc->last_column = 0;} <*>[ \t] <*>. {yylval->val_char = yytext[0]; RETTOK(ERR);} diff --git a/latan/IoAsciiParser.ypp b/latan/IoAsciiParser.ypp index 278f94d..cea7dc2 100644 --- a/latan/IoAsciiParser.ypp +++ b/latan/IoAsciiParser.ypp @@ -23,6 +23,8 @@ #include #include #include + #include + #include #include using namespace std; @@ -52,6 +54,7 @@ %token INT %token ID %token MAT +%token RG_STATE %token OPEN %{ @@ -74,26 +77,64 @@ datas: /* empty string */ - | datas mat + | datas data ; +data: + mat + | rg_state + ; + mat: OPEN MAT ID INT floats CLOSE MAT { const int nRow = state->doubleBuf.size()/$INT, nCol = $INT; - (*state->data)[$ID] = new DMat(nRow,nCol); - DMat &M = static_cast(*((*state->data)[$ID])); - int r,i,j; + (*state->data)[$ID] = new DMat(nRow, nCol); + DMat &m = static_cast(*((*state->data)[$ID])); + int r, i, j; r = 0; while (!state->doubleBuf.empty()) { j = r % nCol; i = (r - j)/nCol; - M(i,j) = state->doubleBuf.top(); + m(i, j) = state->doubleBuf.top(); state->doubleBuf.pop(); ++r; } + if (r != nRow*nCol) + { + LATAN_ERROR(Range, "matrix '" + *state->streamName + ":" + $ID + + "' has a wrong size"); + } + } + ; + +rg_state: + OPEN RG_STATE ID ints CLOSE RG_STATE + { + (*state->data)[$ID] = new RandGen::State; + RandGen::State &rgState = + static_cast(*((*state->data)[$ID])); + + for (int i = 0; i < RLXG_STATE_SIZE; ++i) + { + if (!state->intBuf.empty()) + { + rgState[i] = state->intBuf.top(); + state->intBuf.pop(); + } + else + { + LATAN_ERROR(Range, "random generator state '" + + *state->streamName + ":" + $ID + "' is too short"); + } + } + if (!state->intBuf.empty()) + { + LATAN_ERROR(Range, "random generator state '" + + *state->streamName + ":" + $ID + "' is too long"); + } } ; @@ -103,3 +144,8 @@ floats: | FLOAT {state->doubleBuf.push($1);} | INT {state->doubleBuf.push(static_cast($1));} ; + +ints: + INT ints {state->intBuf.push($1);} + | INT {state->intBuf.push($1);} + ; diff --git a/latan/IoObject.hpp b/latan/IoObject.hpp index 96bdd61..252a1a2 100644 --- a/latan/IoObject.hpp +++ b/latan/IoObject.hpp @@ -33,9 +33,10 @@ public: public: enum { - noType = 0, - dMat = 1, - sample = 2 + noType = 0, + dMat = 1, + sample = 2, + rgState = 3 }; }; public: diff --git a/latan/Makefile.am b/latan/Makefile.am index fd30603..6520e90 100644 --- a/latan/Makefile.am +++ b/latan/Makefile.am @@ -39,6 +39,7 @@ liblatan_la_SOURCES = \ MathParser.ypp \ MathLexer.lpp \ Plot.cpp \ + RandGen.cpp \ Sample.cpp \ ../config.h liblatan_ladir = $(pkgincludedir) @@ -53,6 +54,7 @@ liblatan_la_HEADERS = \ MathInterpreter.hpp \ ParserState.hpp \ Plot.hpp \ + RandGen.hpp \ Sample.hpp liblatan_la_CFLAGS = $(COM_CFLAGS) liblatan_la_CXXFLAGS = $(COM_CXXFLAGS) diff --git a/latan/RandGen.cpp b/latan/RandGen.cpp new file mode 100644 index 0000000..344ac65 --- /dev/null +++ b/latan/RandGen.cpp @@ -0,0 +1,696 @@ +/* + * RandGen.cpp, part of LatAnalyze 3 + * + * Copyright (C) 2013 - 2014 Antonin Portelli + * + * LatAnalyze 3 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. + * + * LatAnalyze 3 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 LatAnalyze 3. If not, see . + */ + +#include +#include + +#ifndef RLXD_LEVEL +#define RLXD_LEVEL 1 +#endif + +using namespace std; +using namespace Latan; + +/****************************************************************************** + * RandGen implementation * + ******************************************************************************/ +// State implentation ////////////////////////////////////////////////////////// +RandGen::State::State(void) +{} + +RandGen::State::~State(void) +{} + +unsigned int RandGen::State::getType(void) const +{ + return IoType::rgState; +} + +// RanLxd implementation /////////////////////////////////////////////////////// +RandGen::RanLxd::RanLxd(void) +: init(0) +{ + // avoid a warning in the SSE case + one_bit = 0.0; +} + +/****************************************************************************/ +/* Copyright (C) 2005 Martin Luescher (GPL) + * This software is distributed under the terms of the GNU General Public + * License (GPL) + * + * Random number generator "ranlxd". See the notes + * + * "User's guide for ranlxs and ranlxd v3.2" (December 2005) + * + * "Algorithms used in ranlux v3.0" (May 2001) + * + * for a detailed description + * + * The functions are + * + * void ranlxd(double r[],int n) + * Computes the next n double-precision random numbers and + * assigns them to the elements r[0],...,r[n-1] of the array r[] + * + * void rlxd_init(int level,int seed) + * Initialization of the generator + * + * int rlxd_size(void) + * Returns the number of integers required to save the state of + * the generator + * + * void rlxd_get(int state[]) + * Extracts the current state of the generator and stores the + * information in the array state[N] where N>=rlxd_size() + * + * void rlxd_reset(int state[]) + * Resets the generator to the state defined by the array state[N] + */ + +#ifdef HAVE_SSE + +#define STEP(pi,pj) \ +__asm__ __volatile__ ("movaps %4, %%xmm4 \n\t" \ + "movaps %%xmm2, %%xmm3 \n\t" \ + "subps %2, %%xmm4 \n\t" \ + "movaps %%xmm1, %%xmm5 \n\t" \ + "cmpps $0x6, %%xmm4, %%xmm2 \n\t" \ + "andps %%xmm2, %%xmm5 \n\t" \ + "subps %%xmm3, %%xmm4 \n\t" \ + "andps %%xmm0, %%xmm2 \n\t" \ + "addps %%xmm4, %%xmm5 \n\t" \ + "movaps %%xmm5, %0 \n\t" \ + "movaps %5, %%xmm6 \n\t" \ + "movaps %%xmm2, %%xmm3 \n\t" \ + "subps %3, %%xmm6 \n\t" \ + "movaps %%xmm1, %%xmm7 \n\t" \ + "cmpps $0x6, %%xmm6, %%xmm2 \n\t" \ + "andps %%xmm2, %%xmm7 \n\t" \ + "subps %%xmm3, %%xmm6 \n\t" \ + "andps %%xmm0, %%xmm2 \n\t" \ + "addps %%xmm6, %%xmm7 \n\t" \ + "movaps %%xmm7, %1" \ + : \ + "=m" ((*pi).c1), \ + "=m" ((*pi).c2) \ + : \ + "m" ((*pi).c1), \ + "m" ((*pi).c2), \ + "m" ((*pj).c1), \ + "m" ((*pj).c2) \ + : \ + "xmm2", "xmm3", "xmm4", "xmm5", "xmm6", "xmm7") + +void RandGen::RanLxd::error(int no) const +{ + switch(no) + { + case 1: + LATAN_ERROR(Range, "Bad choice of luxury level (should be 1 or 2)"); + break; + case 2: + LATAN_ERROR(Range, "Bad choice of seed (should be between 1 and 2^31-1)"); + break; + case 3: + LATAN_ERROR(Runtime, "Undefined state (ranlxd is not initialized)"); + break; + case 5: + LATAN_ERROR(Logic, "Unexpected input data"); + break; + } +} + +void RandGen::RanLxd::update(void) +{ + int k,kmax; + rlxd_dble_vec_t *pmin,*pmax,*pi,*pj; + + kmax=rlxd_pr; + pmin=&rlxd_x.vec[0]; + pmax=pmin+12; + pi=&rlxd_x.vec[ir]; + pj=&rlxd_x.vec[jr]; + + __asm__ __volatile__ ("movaps %0, %%xmm0 \n\t" + "movaps %1, %%xmm1 \n\t" + "movaps %2, %%xmm2" + : + : + "m" (one_bit_sse), + "m" (one_sse), + "m" (carry) + : + "xmm0", "xmm1", "xmm2"); + + for (k=0;k=12) + ir-=12; + if (jr>=12) + jr-=12; + is=8*ir; + is_old=is; +} + +void RandGen::RanLxd::define_constants(void) +{ + int k; + float b; + + one_sse.c1=1.0f; + one_sse.c2=1.0f; + one_sse.c3=1.0f; + one_sse.c4=1.0f; + + b=(float)(ldexp(1.0,-24)); + one_bit_sse.c1=b; + one_bit_sse.c2=b; + one_bit_sse.c3=b; + one_bit_sse.c4=b; + + for (k=0;k<96;k++) + { + next[k]=(k+1)%96; + if ((k%4)==3) + next[k]=(k+5)%96; + } +} + +void RandGen::RanLxd::rlxd_init(int level,int seed) +{ + int i,k,l; + int ibit,jbit,xbit[31]; + int ix,iy; + + define_constants(); + + if (level==1) + rlxd_pr=202; + else if (level==2) + rlxd_pr=397; + else + error(1); + + i=seed; + + for (k=0;k<31;k++) + { + xbit[k]=i%2; + i/=2; + } + + if ((seed<=0)||(i!=0)) + error(2); + + ibit=0; + jbit=18; + + for (i=0;i<4;i++) + { + for (k=0;k<24;k++) + { + ix=0; + + for (l=0;l<24;l++) + { + iy=xbit[ibit]; + ix=2*ix+iy; + + xbit[ibit]=(xbit[ibit]+xbit[jbit])%2; + ibit=(ibit+1)%31; + jbit=(jbit+1)%31; + } + + if ((k%4)!=i) + ix=16777215-ix; + + rlxd_x.num[4*k+i]=(float)(ldexp((double)(ix),-24)); + } + } + + carry.c1=0.0f; + carry.c2=0.0f; + carry.c3=0.0f; + carry.c4=0.0f; + + ir=0; + jr=7; + is=91; + is_old=0; + prm=rlxd_pr%12; + init=1; +} + +void RandGen::RanLxd::ranlxd(double r[],int n) +{ + int k; + + if (init==0) + rlxd_init(1,1); + + for (k=0;k=167777216)) + error(5); + + rlxd_x.num[k]=(float)(ldexp((double)(state[k+1]),-24)); + } + + if (((state[97]!=0)&&(state[97]!=1))|| + ((state[98]!=0)&&(state[98]!=1))|| + ((state[99]!=0)&&(state[99]!=1))|| + ((state[100]!=0)&&(state[100]!=1))) + error(5); + + carry.c1=(float)(ldexp((double)(state[97]),-24)); + carry.c2=(float)(ldexp((double)(state[98]),-24)); + carry.c3=(float)(ldexp((double)(state[99]),-24)); + carry.c4=(float)(ldexp((double)(state[100]),-24)); + + rlxd_pr=state[101]; + ir=state[102]; + jr=state[103]; + is=state[104]; + is_old=8*ir; + prm=rlxd_pr%12; + init=1; + + if (((rlxd_pr!=202)&&(rlxd_pr!=397))|| + (ir<0)||(ir>11)||(jr<0)||(jr>11)||(jr!=((ir+7)%12))|| + (is<0)||(is>91)) + error(5); +} + +#else + +#define BASE 0x1000000 +#define MASK 0xffffff + +#define STEP(pi,pj) \ +d=(*pj).c1.c1-(*pi).c1.c1-carry.c1; \ +(*pi).c2.c1+=(d<0); \ +d+=BASE; \ +(*pi).c1.c1=d&MASK; \ +d=(*pj).c1.c2-(*pi).c1.c2-carry.c2; \ +(*pi).c2.c2+=(d<0); \ +d+=BASE; \ +(*pi).c1.c2=d&MASK; \ +d=(*pj).c1.c3-(*pi).c1.c3-carry.c3; \ +(*pi).c2.c3+=(d<0); \ +d+=BASE; \ +(*pi).c1.c3=d&MASK; \ +d=(*pj).c1.c4-(*pi).c1.c4-carry.c4; \ +(*pi).c2.c4+=(d<0); \ +d+=BASE; \ +(*pi).c1.c4=d&MASK; \ +d=(*pj).c2.c1-(*pi).c2.c1; \ +carry.c1=(d<0); \ +d+=BASE; \ +(*pi).c2.c1=d&MASK; \ +d=(*pj).c2.c2-(*pi).c2.c2; \ +carry.c2=(d<0); \ +d+=BASE; \ +(*pi).c2.c2=d&MASK; \ +d=(*pj).c2.c3-(*pi).c2.c3; \ +carry.c3=(d<0); \ +d+=BASE; \ +(*pi).c2.c3=d&MASK; \ +d=(*pj).c2.c4-(*pi).c2.c4; \ +carry.c4=(d<0); \ +d+=BASE; \ +(*pi).c2.c4=d&MASK + +void RandGen::RanLxd::error(int no) const +{ + switch(no) + { + case 0: + LATAN_ERROR(System, "Arithmetic on this machine is not suitable for ranlxd"); + break; + case 1: + LATAN_ERROR(Range, "Bad choice of luxury level (should be 1 or 2)"); + break; + case 2: + LATAN_ERROR(Range, "Bad choice of seed (should be between 1 and 2^31-1)"); + break; + case 3: + LATAN_ERROR(Runtime, "Undefined state (ranlxd is not initialized)"); + case 4: + LATAN_ERROR(System, "Arithmetic on this machine is not suitable for ranlxd"); + break; + case 5: + LATAN_ERROR(Logic, "Unexpected input data"); + break; + } +} + +void RandGen::RanLxd::update(void) +{ + int k,kmax,d; + rlxd_dble_vec_t *pmin,*pmax,*pi,*pj; + + kmax=rlxd_pr; + pmin=&rlxd_x.vec[0]; + pmax=pmin+12; + pi=&rlxd_x.vec[ir]; + pj=&rlxd_x.vec[jr]; + + for (k=0;k=12) + ir-=12; + if (jr>=12) + jr-=12; + is=8*ir; + is_old=is; +} + +void RandGen::RanLxd::define_constants(void) +{ + int k; + + one_bit=ldexp(1.0,-24); + + for (k=0;k<96;k++) + { + next[k]=(k+1)%96; + if ((k%4)==3) + next[k]=(k+5)%96; + } +} + +void RandGen::RanLxd::rlxd_init(int level,int seed) +{ + int i,k,l; + int ibit,jbit,xbit[31]; + int ix,iy; + + if ((INT_MAX<2147483647)||(FLT_RADIX!=2)||(FLT_MANT_DIG<24)|| + (DBL_MANT_DIG<48)) + error(0); + + define_constants(); + + if (level==1) + rlxd_pr=202; + else if (level==2) + rlxd_pr=397; + else + error(1); + + i=seed; + + for (k=0;k<31;k++) + { + xbit[k]=i%2; + i/=2; + } + + if ((seed<=0)||(i!=0)) + error(2); + + ibit=0; + jbit=18; + + for (i=0;i<4;i++) + { + for (k=0;k<24;k++) + { + ix=0; + + for (l=0;l<24;l++) + { + iy=xbit[ibit]; + ix=2*ix+iy; + + xbit[ibit]=(xbit[ibit]+xbit[jbit])%2; + ibit=(ibit+1)%31; + jbit=(jbit+1)%31; + } + + if ((k%4)!=i) + ix=16777215-ix; + + rlxd_x.num[4*k+i]=ix; + } + } + + carry.c1=0; + carry.c2=0; + carry.c3=0; + carry.c4=0; + + ir=0; + jr=7; + is=91; + is_old=0; + prm=rlxd_pr%12; + init=1; +} + +void RandGen::RanLxd::ranlxd(double r[],int n) +{ + int k; + + if (init==0) + rlxd_init(1,1); + + for (k=0;k=167777216)) + error(5); + + rlxd_x.num[k]=state[k+1]; + } + + if (((state[97]!=0)&&(state[97]!=1))|| + ((state[98]!=0)&&(state[98]!=1))|| + ((state[99]!=0)&&(state[99]!=1))|| + ((state[100]!=0)&&(state[100]!=1))) + error(5); + + carry.c1=state[97]; + carry.c2=state[98]; + carry.c3=state[99]; + carry.c4=state[100]; + + rlxd_pr=state[101]; + ir=state[102]; + jr=state[103]; + is=state[104]; + is_old=8*ir; + prm=rlxd_pr%12; + init=1; + + if (((rlxd_pr!=202)&&(rlxd_pr!=397))|| + (ir<0)||(ir>11)||(jr<0)||(jr>11)||(jr!=((ir+7)%12))|| + (is<0)||(is>91)) + error(5); +} + +#endif + +// constructors //////////////////////////////////////////////////////////////// +RandGen::RandGen(void) +{ + generator_.rlxd_init(RLXD_LEVEL, (int)time(NULL)); +} + +RandGen::RandGen(const int seed) +{ + generator_.rlxd_init(RLXD_LEVEL, seed); +} + +RandGen::RandGen(const State &state) +{ + setState(state); +} + +// destructor ////////////////////////////////////////////////////////////////// +RandGen::~RandGen(void) +{} + +// state management //////////////////////////////////////////////////////////// +RandGen::State RandGen::getState(void) const +{ + State state; + + generator_.rlxd_get(state.data()); + + return state; +} + +void RandGen::setState(const State &state) +{ + generator_.rlxd_reset(state.data()); +} + +// generators ////////////////////////////////////////////////////////////////// +double RandGen::uniform(const double a, const double b) +{ + double rx; + + generator_.ranlxd(&rx, 1); + + return (b-a)*rx + a; +} + +double RandGen::discreteUniform(const unsigned int n) +{ + return ((unsigned int)(uniform()*(double)(n))); +} + +double RandGen::gaussian(const double mean, const double sigma) +{ + double rx, ry, sqNrm; + + do + { + rx = uniform(-1.0, 1.0); + ry = uniform(-1.0, 1.0); + sqNrm = rx*rx + ry*ry; + } while ((sqNrm > 1.0)||(sqNrm == 0.0)); + + return sigma*rx*sqrt(-2.0*log(sqNrm)/sqNrm) + mean; +} diff --git a/latan/RandGen.hpp b/latan/RandGen.hpp new file mode 100644 index 0000000..dad08ac --- /dev/null +++ b/latan/RandGen.hpp @@ -0,0 +1,101 @@ +/* + * RandGen.hpp, part of LatAnalyze 3 + * + * Copyright (C) 2013 - 2014 Antonin Portelli + * + * LatAnalyze 3 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. + * + * LatAnalyze 3 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 LatAnalyze 3. If not, see . + */ + +#ifndef Latan_RandGen_hpp_ +#define Latan_RandGen_hpp_ + +#include +#include + +#define RLXG_STATE_SIZE 105 + +BEGIN_NAMESPACE + +/****************************************************************************** + * Random generator class * + ******************************************************************************/ + +class RandGen +{ +public: + class State: public Eigen::Array, public IoObject + { + public: + // constructor + State(void); + // destructor + ~State(void); + // IO type + unsigned int getType(void) const; + }; +private: + // Martin Luescher's ranlxd generator interface + class RanLxd + { + private: + typedef struct + { + float c1,c2,c3,c4; + } rlxd_vec_t __attribute__ ((aligned (16))); + typedef struct + { + rlxd_vec_t c1,c2; + } rlxd_dble_vec_t __attribute__ ((aligned (16))); + public: + RanLxd(void); + void ranlxd(double r[],int n); + void rlxd_init(int level,int seed); + int rlxd_size(void) const; + void rlxd_get(int state[]) const; + void rlxd_reset(const int state[]); + private: + void error(int no) const; + void update(void); + void define_constants(void); + private: + int init, rlxd_pr, prm, ir, jr, is, is_old, next[96]; + rlxd_vec_t one_sse, one_bit_sse, carry; + double one_bit; + union + { + rlxd_dble_vec_t vec[12]; + float num[96]; + } rlxd_x __attribute__ ((aligned (16))); + }; +public: + // constructors + RandGen(void); + RandGen(const int seed); + RandGen(const State &state); + // destructor + virtual ~RandGen(void); + // state management + State getState(void) const; + void setState(const State &state); + // generators + double uniform(const double a = 0.0, const double b = 1.0); + double discreteUniform(const unsigned int n); + double gaussian(const double mean = 0.0, const double sigma = 1.0); +private: + RanLxd generator_; +}; + +END_NAMESPACE + +#endif // Latan_RandGen_hpp_ diff --git a/latan/includes.hpp b/latan/includes.hpp index eaba822..6e5d2ba 100644 --- a/latan/includes.hpp +++ b/latan/includes.hpp @@ -24,6 +24,8 @@ #include #include #include +#include +#include #include #include #include