CLHEP 2.4.6.4
C++ Class Library for High Energy Physics
Loading...
Searching...
No Matches
Macros | Functions
gaussSpeed.cc File Reference
#include "CLHEP/Random/Randomize.h"
#include "CLHEP/Random/defs.h"
#include <iostream>
#include "CLHEP/Random/RandGaussQ.h"
#include "CLHEP/Random/RandGaussT.h"
#include "CLHEP/Random/RandPoissonQ.h"
#include "CLHEP/Random/RandPoissonT.h"
#include "CLHEP/Random/RandBit.h"
#include "CLHEP/Units/PhysicalConstants.h"
#include <iomanip>

Go to the source code of this file.

Macros

#define MISC
 

Functions

double gammln1 (double xx)
 
double gammln2 (double xx)
 
int main ()
 

Macro Definition Documentation

◆ MISC

#define MISC

Function Documentation

◆ gammln1()

double gammln1 ( double  xx)

Definition at line 29 of file gaussSpeed.cc.

29 {
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}

Referenced by main().

◆ gammln2()

double gammln2 ( double  xx)

Definition at line 51 of file gaussSpeed.cc.

51 {
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}

Referenced by main().

◆ main()

int main ( )

Definition at line 76 of file gaussSpeed.cc.

76 {
77
78 cout << "Enter 1 for RandGauss, 2 for RandGaussQ, 3 for DualRand flat: ";
79 int choice;
80 cin >> choice;
81
82if (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
108if (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
134if (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}
CLHEP::RandBit::shootBit
static int shootBit()
CLHEP::RandFlat::shootBit
static int shootBit()
CLHEP::RandGaussQ::shoot
static double shoot()
CLHEP::RandGaussT::shoot
static double shoot()
CLHEP::RandPoissonQ::shoot
static long shoot(double mean=1.0)
Definition: RandPoissonQ.cc:121
CLHEP::RandPoissonT::shoot
static long shoot(double mean=1.0)
Definition: RandPoissonT.cc:60
gammln2
double gammln2(double xx)
Definition: gaussSpeed.cc:51
gammln1
double gammln1(double xx)
Definition: gaussSpeed.cc:29