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 
CLHEP::RandPoissonT
Definition: Matrix/CLHEP/Random/RandPoissonT.h:41
CLHEP::RandPoissonQ::fire
long fire()
Definition: RandPoissonQ.cc:134
CLHEP::RandGaussT::shoot
static double shoot()
CLHEP::RandGaussQ::fire
double fire()
b
@ b
Definition: testCategories.cc:125
CLHEP::RandPoisson::fire
long fire()
Definition: RandPoisson.cc:214
mu
the goal is to keep the overall false rejection probability down at the to level For each validated we discuss which of course is by necessity relative timing We take the time for a single random via one of the fastest good and at any rate the ratios will vary by around depending on the processor and memory configuration used A timing for a distribution of units would mean no time used beyond the uniform random Summary Distribution Validated Validation Rejected Past N RandGauss RandGaussT RandGaussQ bins stepwise bins RandPoisson RandPoissonT mu< 100 N=50, 000, 000 ------- mu > RandPoissonQ mu< 100 N=50, 000, 000 -------(same as RandPoissonT) mu > RandGauss shoot() etc 2.5 units Validation tests applied and a very accurate series is substituted for the table method shoot() etc 1.7 units Validation tests applied and below about **a series approximation is but we have applied this test with N up to and never get anywhere near the rejection point Analytical considerations indicate that the validation test would not reject until O(10 **24) samples were inspected. ---------------------------------------------------------- 2. RandGeneral Method since we wish to have good resolution power even if just one of the mu values is it would be unwise to dial this down any further Validation for several values o mu)
Definition: validation.doc:263
CLHEP::RandGaussQ
Definition: Matrix/CLHEP/Random/RandGaussQ.h:32
CLHEP::DualRand::flat
double flat()
Definition: DualRand.cc:107
CLHEP::RandGaussQ::shoot
static double shoot()
CLHEP::RandPoissonQ::shoot
static long shoot(double m=1.0)
Definition: RandPoissonQ.cc:118
CLHEP::RandGaussT::fire
double fire()
N
the goal is to keep the overall false rejection probability down at the to level For each validated we discuss which of course is by necessity relative timing We take the time for a single random via one of the fastest good and at any rate the ratios will vary by around depending on the processor and memory configuration used A timing for a distribution of units would mean no time used beyond the uniform random Summary Distribution Validated Validation Rejected Past N RandGauss N
Definition: validation.doc:48
CLHEP::RandBit
Definition: Matrix/CLHEP/Random/RandBit.h:40
CLHEP::RandBit::shootBit
static int shootBit()
CLHEP
Definition: ClhepVersion.h:13
CLHEP::RandPoissonT::shoot
static long shoot(double m=1.0)
Definition: RandPoissonT.cc:56
CLHEP::RandGaussT
Definition: Matrix/CLHEP/Random/RandGaussT.h:41
gammln2
double gammln2(double xx)
Definition: gaussSpeed.cc:51
j
long j
Definition: JamesRandomSeeding.txt:28
CLHEP::DualRand
Definition: Matrix/CLHEP/Random/DualRand.h:51
gammln1
double gammln1(double xx)
Definition: gaussSpeed.cc:29
i
long i
Definition: JamesRandomSeeding.txt:27
CLHEP::RandGauss
Definition: Matrix/CLHEP/Random/RandGauss.h:42
CLHEP::RandPoissonT::fire
long fire()
Definition: RandPoissonT.cc:73
CLHEP::RandPoissonQ
Definition: Matrix/CLHEP/Random/RandPoissonQ.h:33
CLHEP::RandPoisson
Definition: Matrix/CLHEP/Random/RandPoisson.h:42
CLHEP::RandFlat::shoot
static double shoot()
Definition: RandFlat.cc:60
x
any side effects of that construction would occur twice The semantics of throw x
Definition: whyZMthrowRethrows.txt:37
main
int main()
Definition: gaussSpeed.cc:76
CLHEP::RandGauss::fire
double fire()
CLHEP::RandFlat::shootBit
static int shootBit()