Commit e9600998 authored by Lukas Weber's avatar Lukas Weber
Browse files

implement my own version of normal sampling

parent adfc54aa
...@@ -13,8 +13,6 @@ using namespace std; ...@@ -13,8 +13,6 @@ using namespace std;
class abstract_mc class abstract_mc
{ {
private: private:
randomnumbergenerator * rng = 0;
void param_init(string dir); void param_init(string dir);
void random_clear(); void random_clear();
void random_write(odump& d); void random_write(odump& d);
...@@ -26,6 +24,7 @@ class abstract_mc ...@@ -26,6 +24,7 @@ class abstract_mc
int therm = 0; int therm = 0;
protected: protected:
parser param; parser param;
randomnumbergenerator * rng = 0;
virtual void init() = 0; virtual void init() = 0;
virtual void write(odump &out) = 0; virtual void write(odump &out) = 0;
......
...@@ -32,6 +32,32 @@ void randomnumbergenerator :: read(idump& d) ...@@ -32,6 +32,32 @@ void randomnumbergenerator :: read(idump& d)
ptr->load( randState ); ptr->load( randState );
delete [] randState; delete [] randState;
} }
double randomnumbergenerator :: norm() { // adjusted from https://de.wikipedia.org/wiki/Polar-Methode
/* for now, I do not need the extra speed coming from this.
* If you add hidden state like this you have to change the dump format
* (breaking compatibility with your old runs)
if(norm_is_stored) {
norm_is_stored=false;
return norm_stored;
}
*/
double u, v, q, p;
do {
u = 2.0 * d() - 1;
v = 2.0 * d() - 1;
q = u * u + v * v;
} while (q <= 0.0 || q >= 1.0);
p = sqrt(-2 * log(q) / q);
/*
norm_stored = u * p;
norm_is_stored = true;
*/
return v * p;
}
#endif #endif
#ifdef MCL_RNG_SPRNG_4 #ifdef MCL_RNG_SPRNG_4
......
...@@ -21,6 +21,8 @@ public: ...@@ -21,6 +21,8 @@ public:
double d(double supp) {return ptr->randDblExc(supp);} //in (0,supp) double d(double supp) {return ptr->randDblExc(supp);} //in (0,supp)
int i(int bound) {return ptr->randInt(bound-1);} //in [0,bound) int i(int bound) {return ptr->randInt(bound-1);} //in [0,bound)
double norm(); // normal distribution, mean = 0, std = 1
private: private:
MTRand* ptr; MTRand* ptr;
luint my_seed; luint my_seed;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment