random.cpp 2.67 KB
Newer Older
Stefan Weßel's avatar
ic  
Stefan Weßel committed
1
2
#include "random.h"

3
4
5
6
7
8
#ifdef MCL_RNG_MT
randomnumbergenerator :: randomnumbergenerator()
{
	ptr = new MTRand();
	my_seed=ptr->get_seed();
}
Stefan Weßel's avatar
ic  
Stefan Weßel committed
9

10
11
12
13
14
randomnumbergenerator :: randomnumbergenerator(luint seed)
{
	ptr = new MTRand(seed);
	my_seed=ptr->get_seed();
}
Stefan Weßel's avatar
ic  
Stefan Weßel committed
15

16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
void randomnumbergenerator :: write(odump& d)
{
	MTRand::uint32* randState;
	randState = new MTRand::uint32[ MTRand::SAVE ];
	ptr->save( randState );
	d.write(randState,MTRand::SAVE);
	d.write(my_seed);
	delete [] randState;
}

void randomnumbergenerator :: read(idump& d)
{
	MTRand::uint32* randState;
	randState = new MTRand::uint32[ MTRand::SAVE ];
	d.read(randState,MTRand::SAVE);
	d.read(my_seed);
	ptr->load( randState );
	delete [] randState;
}
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

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;
}

Stefan Weßel's avatar
ic  
Stefan Weßel committed
61
62
#endif

63
64
65
66
67
68
69
#ifdef  MCL_RNG_SPRNG_4
randomnumbergenerator :: randomnumbergenerator()
{
	int gtype = 4; //Multipl. Lagg. Fib. 
	ptr = SelectType(gtype);
	my_seed=make_sprng_seed();
	ptr->init_sprng(0, 1, my_seed, SPRNG_DEFAULT);
Stefan Weßel's avatar
ic  
Stefan Weßel committed
70
71
}

72
randomnumbergenerator :: randomnumbergenerator(luint seed)
Stefan Weßel's avatar
ic  
Stefan Weßel committed
73
{
74
75
76
77
	int gtype = 4; //Multipl. Lagg. Fib.
	ptr = SelectType(gtype);
	my_seed=seed;
	ptr->init_sprng(0, 1, my_seed, SPRNG_DEFAULT);
Stefan Weßel's avatar
ic  
Stefan Weßel committed
78
79
}

80
void randomnumbergenerator :: write(odump& d)
Stefan Weßel's avatar
ic  
Stefan Weßel committed
81
{
82
83
84
85
86
87
	char* buffer;
	int buffersize=ptr->pack_sprng(&buffer);
	d.write(buffersize);
	d.write(buffer,buffersize);
	d.write(my_seed);
	delete buffer;
Stefan Weßel's avatar
ic  
Stefan Weßel committed
88
89
}

90
91
92
93
94
95
96
97
98
99
100
101
102
103
void randomnumbergenerator :: read(idump& d)
{
	char* buffer;
	int buffersize;
	d.read(buffersize);
	buffer = new char[buffersize];
	d.read(buffer,buffersize);
	d.read(my_seed);
	ptr->unpack_sprng(buffer);
	delete buffer;
}
#endif

#ifdef MCL_RNG_BOOST
Stefan Weßel's avatar
ic  
Stefan Weßel committed
104
105
106
107
108
109
110
111
void randomnumbergenerator :: write(odump& d)
{
}

void randomnumbergenerator :: read(idump& d)
{
}

112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
randomnumbergenerator :: randomnumbergenerator()
{
	rng = new boost::mt19937;
	val = new boost::uniform_real<> (0.0, 1.0);
	
	my_seed = 5489;
}

randomnumbergenerator :: randomnumbergenerator(int seed)
{
	rng = new boost::mt19937;
	rng->seed(seed);
	val = new boost::uniform_real<> (0.0, 1.0);
	
	my_seed = seed;
}
#endif
// 
// #ifdef MCL_RNG_ACML
// randomnumbergenerator :: randomnumbergenerator(int seed)
// {
// 	gen = new acml_rand(seed, 0., 1.);
// }
//
// #endif