CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

gaussSpeed.cc
Go to the documentation of this file.
1 #include "CLHEP/Random/Randomize.h"
2 #include "CLHEP/Random/defs.h"
3 #include <iostream>
4 
5 #include "CLHEP/Random/RandGaussQ.h"
6 #include "CLHEP/Random/RandGaussT.h"
7 #include "CLHEP/Random/RandPoissonQ.h"
8 #include "CLHEP/Random/RandPoissonT.h"
9 #include "CLHEP/Random/RandBit.h"
10 #include "CLHEP/Units/PhysicalConstants.h"
11 
12 using std::cin;
13 using std::cout;
14 using std::cerr;
15 using std::endl;
16 using namespace CLHEP;
17 //#ifndef _WIN32
18 //using std::exp;
19 //#endif
20 
21 
22 // ---------
23 // RandGauss
24 //
25 // mf 12/13/04 Correction in way engines are supplied to RandBit() ctor
26 // (gcc3.3.1 exposed previously innocuous mistake)
27 // ---------
28 
29 double gammln1(double xx) {
30 
31 // Returns the value ln(Gamma(xx) for xx > 0. Full accuracy is obtained for
32 // xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first.
33 // (Adapted from Numerical Recipes in C)
34 
35  static double cof[6] = {76.18009172947146,-86.50532032941677,
36  24.01409824083091, -1.231739572450155,
37  0.1208650973866179e-2, -0.5395239384953e-5};
38  int j;
39  double x = xx - 1.0;
40  double tmp = x + 5.5;
41  tmp -= (x + 0.5) * std::log(tmp);
42  double ser = 1.000000000190015;
43 
44  for ( j = 0; j <= 5; j++ ) {
45  x += 1.0;
46  ser += cof[j]/x;
47  }
48  return -tmp + std::log(2.5066282746310005*ser);
49 }
50 
51 double gammln2(double xx) {
52 
53 // Returns the value ln(Gamma(xx) for xx > 0. Full accuracy is obtained for
54 // xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first.
55 // (Adapted from Numerical Recipes in C)
56 
57  static double cof[6] = {76.18009172947146,-86.50532032941677,
58  24.01409824083091, -1.231739572450155,
59  0.1208650973866179e-2, -0.5395239384953e-5};
60  int j;
61  double x = xx - 0.0;
62  double tmp = x + 5.5;
63  tmp -= (x + 0.5) * std::log(tmp);
64  double ser = 1.000000000190015;
65 
66  for ( j = 0; j <= 5; j++ ) {
67  x += 1.0;
68  ser += cof[j]/x;
69  }
70  return -tmp + std::log(2.5066282746310005*ser/xx);
71 }
72 #include <iomanip>
73 
74 
75 
76 int main() {
77 
78  cout << "Enter 1 for RandGauss, 2 for RandGaussQ, 3 for DualRand flat: ";
79  int choice;
80  cin >> choice;
81 
82 if (choice==1) {
83  cout << "\n--------------------------------------------\n";
84  cout << "Test of Gauss distribution speed\n\n";
85 
86  long seed;
87  cout << "Please enter an integer seed: ";
88  cin >> seed;
89 
90  int nNumbers;
91  cout << "How many numbers should we generate: ";
92  cin >> nNumbers;
93 
94  cout << "\nInstantiating distribution utilizing DualRand engine...\n";
95  DualRand eng (seed);
96  RandGauss dist (eng);
97 
98  double sum = 0;
99 
100 
101  for (int i=0; i < nNumbers; i++) {
102  sum += dist.fire();
103  }
104 
105  cout << "\n Finished: sum is " << sum << " \n";
106 }
107 
108 if (choice==2) {
109  cout << "\n--------------------------------------------\n";
110  cout << "Test of RandGaussQ distribution speed\n\n";
111 
112  long seed;
113  cout << "Please enter an integer seed: ";
114  cin >> seed;
115 
116  int nNumbers;
117  cout << "How many numbers should we generate: ";
118  cin >> nNumbers;
119 
120  cout << "\nInstantiating distribution utilizing DualRand engine...\n";
121  DualRand eng (seed);
122  RandGaussQ dist (eng);
123 
124  double sum = 0;
125 
126 
127  for (int i=0; i < nNumbers; i++) {
128  sum += dist.fire();
129  }
130 
131  cout << "\n Finished: sum is " << sum << " \n";
132 }
133 
134 if (choice==3) {
135  cout << "\n--------------------------------------------\n";
136  cout << "Test of DualRand flat speed\n\n";
137 
138  long seed;
139  cout << "Please enter an integer seed: ";
140  cin >> seed;
141 
142  int nNumbers;
143  cout << "How many numbers should we generate: ";
144  cin >> nNumbers;
145 
146  cout << "\nInstantiating distribution utilizing DualRand engine...\n";
147  DualRand eng (seed);
148 
149  double sum = 0;
150 
151 
152  for (int i=0; i < nNumbers; i++) {
153  sum += eng.flat();
154  }
155 
156  cout << "\n Finished: sum is " << sum << " \n";
157 }
158 
159 
160 #ifdef GAMMA
161  cout << "\nNow we will compute the first 20 gammas, using gammln:\n";
162 
163  double x;
164  for (x=1; x <= 20; x+=1) {
165  cout << x << std::setprecision(20) << " " << std::exp(gammln1(x))
166  << " " << std::exp(gammln2(x)) << " difference in gammln2 = " <<
167  gammln1(x)-gammln2(x) << "\n";
168  }
169 
170 
171  cout << "\nNow we will compute gamma of small numbers: \n";
172 
173  for ( x=1; x > .000000001; x *= .9 ) {
174  cout << x << std::setprecision(20) << " " <<
175  1 - std::exp(gammln1(x)) * std::exp(gammln1(2-x)) * std::sin(CLHEP::pi*(1-x)) / (CLHEP::pi*(1-x)) <<
176  " " <<
177  1 - std::exp(gammln2(x)) * std::exp(gammln1(2-x)) * std::sin(CLHEP::pi*(1-x)) / (CLHEP::pi*(1-x)) <<
178  "\n";
179  }
180 #endif // GAMMA
181 
182 #ifdef POISSON
183  cout << "\n--------------------------------------------\n";
184  cout << "Test of Poisson distribution speed\n\n";
185 
186  long seed;
187  cout << "Please enter an integer seed: ";
188  cin >> seed;
189 
190  double mu;
191  cout << "Please enter mu: ";
192  cin >> mu;
193 
194  int nNumbers;
195  cout << "How many numbers should we generate: ";
196  cin >> nNumbers;
197 
198  cout << "\nInstantiating distribution utilizing DualRand engine...\n";
199  DualRand eng (seed);
200  RandPoisson dist (eng, mu);
201  // RandFlat dist (eng);
202 
203  double sum = 0;
204 
205 
206  for (int i=0; i < nNumbers; i++) {
207  sum += dist.fire();
208 // sum += dist.quick();
209 // sum += dist.fire(mu);
210 // sum += dist.quick(mu);
211 
212  }
213 
214  cout << "\n Finished: sum is " << sum << " \n";
215 #endif // POISSON
216 
217 
218 #define MISC
219 #ifdef MISC
220 
221  DualRand e;
222 
223  // RandGauss usage modes
224 
225  cout << "testing RandGaussT::shoot(): " << RandGaussT::shoot() << "\n";
226  cout << "testing RandGaussT::shoot(&e): " << RandGaussT::shoot(&e) << "\n";
227  cout << "testing RandGaussT::shoot(100,10): " <<
228  RandGaussT::shoot(100,10) << "\n";
229  cout << "testing RandGaussT::shoot(&e,100,10): " <<
230  RandGaussT::shoot(&e,100,10) << "\n";
231  RandGaussT gt (e, 50,2);
232  cout << "testing gt.fire(): " << gt.fire() << "\n";
233  cout << "testing gt.fire(200,2): " << gt.fire(200,2) << "\n";
234 
235  cout << "testing RandGaussQ::shoot(): " << RandGaussQ::shoot() << "\n";
236  cout << "testing RandGaussQ::shoot(&e): " << RandGaussQ::shoot(&e) << "\n";
237  cout << "testing RandGaussQ::shoot(100,10): " <<
238  RandGaussQ::shoot(100,10) << "\n";
239  cout << "testing RandGaussQ::shoot(&e,100,10): " <<
240  RandGaussQ::shoot(&e,100,10) << "\n";
241  RandGaussQ qt (e, 50,2);
242  cout << "testing qt.fire(): " << qt.fire() << "\n";
243  cout << "testing qt.fire(200,2): " << qt.fire(200,2) << "\n";
244 
245  // RandPoisson usage modes
246 
247  cout << "testing RandPoissonT::shoot(): " << RandPoissonT::shoot() << "\n";
248  cout << "testing RandPoissonT::shoot(&e): "
249  << RandPoissonT::shoot(&e) << "\n";
250  cout << "testing RandPoissonT::shoot(90): " <<
251  RandPoissonT::shoot(90) << "\n";
252  cout << "testing RandPoissonT::shoot(&e,90): " <<
253  RandPoissonT::shoot(&e,90) << "\n";
254  RandPoissonT pgt (e,50);
255  cout << "testing pgt.fire(): " << pgt.fire() << "\n";
256  cout << "testing pgt.fire(20): " << pgt.fire(20) << "\n";
257 
258  cout << "testing RandPoissonQ::shoot(): " << RandPoissonQ::shoot() << "\n";
259  cout << "testing RandPoissonQ::shoot(&e): " << RandPoissonQ::shoot(&e) << "\n";
260  cout << "testing RandPoissonQ::shoot(90): " <<
261  RandPoissonQ::shoot(90) << "\n";
262  cout << "testing RandPoissonQ::shoot(&e,90): " <<
263  RandPoissonQ::shoot(&e,90) << "\n";
264  RandPoissonQ pqt (e,50);
265  cout << "testing pqt.fire(): " << pqt.fire() << "\n";
266  cout << "testing pqt.fire(20): " << pqt.fire(20) << "\n";
267 
268  // RandBit modes coming from RandFlat and bit modes
269 
270  cout << "testing RandBit::shoot(): " << RandBit::shoot() << "\n";
271  cout << "testing RandBit::shoot(&e): " << RandBit::shoot(&e) << "\n";
272  cout << "testing RandBit::shoot(10,20): " << RandBit::shoot(10,20) << "\n";
273  cout << "testing RandBit::shoot(&e,10,20): "<<
274  RandBit::shoot(&e,10,20) << "\n";
275  RandBit b ( e, 1000, 1100 );
276  cout << "testing b.fire(): " << b.fire() << "\n";
277  cout << "testing b.fire(10,20): " << b.fire(10,20) << "\n";
278  int i;
279  cout << "testing RandBit::shootBit(): ";
280  for (i=0; i<40; i++) {
281  cout << RandBit::shootBit();
282  } cout << "\n";
283  cout << "testing RandBit::shootBit(&e): ";
284  for (i=0; i<40; i++) {
285  cout << RandBit::shootBit(&e);
286  } cout << "\n";
287  cout << "testing RandBit::fireBit(): ";
288  for (i=0; i<40; i++) {
289  cout << b.fireBit();
290  } cout << "\n";
291 
292  // Timing for RandBit:
293 
294  cout << "Timing RandFlat::shootBit(): Enter N: ";
295  int N;
296  cin >> N;
297  int sum=0;
298  for (i=0; i<N; i++) {
299  sum+= RandFlat::shootBit();
300  }
301  cout << "--------- Done.............. Sum = " << sum << "\n";
302  cout << "Timing RandBit::shootBit(): Enter N: ";
303  cin >> N;
304  sum = 0;
305  for (i=0; i<N; i++) {
306  sum += RandBit::shootBit();
307  }
308  cout << "--------- Done.............. Sum = " << sum << "\n";
309 
310 #endif // MISC
311 
312  return 0;
313 }
314