45 #define Table0size 200
46 #define Table1size 250
47 #define Table2size 200
48 #define Table3size 250
49 #define Table4size 1000
50 #define TableSize (Table0size+Table1size+Table2size+Table3size+Table4size)
58 #define Table0step (2.0E-13)
59 #define Table1step (4.0E-11)
60 #define Table2step (1.0E-8)
61 #define Table3step (2.0E-6)
62 #define Table4step (5.0E-4)
70 #define Table0offset 0
71 #define Table1offset (2*(Table0size))
72 #define Table2offset (2*(Table0size + Table1size))
73 #define Table3offset (2*(Table0size + Table1size + Table2size))
74 #define Table4offset (2*(Table0size + Table1size + Table2size + Table3size))
85 #include "CLHEP/Random/gaussTables.cdat"
94 std::cout <<
"r = " << r <<
"\n";
101 std::cout <<
"r = " << r <<
"sign negative \n";
103 }
else if ( r == .5 ) {
118 const double* tptr = 0;
130 if (index <= 0) index = 1;
135 std::cout <<
"index = " << index <<
" dx = " << dx <<
" h = " << h <<
"\n";
142 std::cout <<
"offset index = " << index <<
"\n";
147 }
else if (r <
Tsteps[0]) {
154 for (
int tableN = 3; tableN >= 0; tableN-- ) {
155 if ( r <
Tsteps[tableN] ) {
continue;}
157 std::cout <<
"Using table " << tableN <<
"\n";
165 if (index == 0) index = 1;
166 if (index >=
Tsizes[tableN]) index =
Tsizes[tableN] - 1;
170 std::cout <<
"index = " << index <<
" dx = " << dx <<
" h = " << h <<
"\n";
172 index = (index<<1) +
Toffsets[tableN] - 2;
192 std::cout <<
"y0: " << y0 <<
" y1: " << y1 <<
" d0: " << d0 <<
" d1: " << d1 <<
"\n";
196 double oneMinusX = 1 -
dx;
197 double oneMinusX2 = oneMinusX * oneMinusX;
199 double f0 = (2. * dx + 1.) * oneMinusX2;
200 double f1 = (3. - 2. *
dx) * x2;
201 double g0 = h * dx * oneMinusX2;
202 double g1 = - h * oneMinusX *
x2;
205 std::cout <<
"f0: " << f0 <<
" f1: " << f1 <<
" g0: " << g0 <<
" g1: " << g1 <<
"\n";
208 double answer = f0 * y0 + f1 * y1 + g0 * d0 + g1 *
d1;
211 std::cout <<
"variate is: " << sign*answer <<
"\n";
214 return sign * answer;
243 for (
int i = 1; i < 50; i++ ) {
244 double vn2 = 1.0/(guess*guess);
245 double s1 = -13*11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2*vn2;
246 s1 += 11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2;
247 s1 += -9*7*5*3 * vn2*vn2*vn2*vn2*vn2;
248 s1 += 7*5*3 * vn2*vn2*vn2*vn2;
249 s1 += -5*3 * vn2*vn2*vn2;
250 s1 += 3 * vn2*vn2 - vn2 + 1.0;
251 v = std::sqrt ( 2.0 * std::log ( s1 / (r*guess*std::sqrt(
CLHEP::twopi)) ) );
252 if (
std::abs(v-guess) < eps )
break;