CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 #include "CLHEP/Random/Randomize.h" 00002 #include "CLHEP/Random/defs.h" 00003 #include <iostream> 00004 00005 #include "CLHEP/Random/RandGaussQ.h" 00006 #include "CLHEP/Random/RandGaussT.h" 00007 #include "CLHEP/Random/RandPoissonQ.h" 00008 #include "CLHEP/Random/RandPoissonT.h" 00009 #include "CLHEP/Random/RandBit.h" 00010 #include "CLHEP/Units/PhysicalConstants.h" 00011 00012 using std::cin; 00013 using std::cout; 00014 using std::cerr; 00015 using std::endl; 00016 using namespace CLHEP; 00017 //#ifndef _WIN32 00018 //using std::exp; 00019 //#endif 00020 00021 00022 // --------- 00023 // RandGauss 00024 // 00025 // mf 12/13/04 Correction in way engines are supplied to RandBit() ctor 00026 // (gcc3.3.1 exposed previously innocuous mistake) 00027 // --------- 00028 00029 double gammln1(double xx) { 00030 00031 // Returns the value ln(Gamma(xx) for xx > 0. Full accuracy is obtained for 00032 // xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first. 00033 // (Adapted from Numerical Recipes in C) 00034 00035 static double cof[6] = {76.18009172947146,-86.50532032941677, 00036 24.01409824083091, -1.231739572450155, 00037 0.1208650973866179e-2, -0.5395239384953e-5}; 00038 int j; 00039 double x = xx - 1.0; 00040 double tmp = x + 5.5; 00041 tmp -= (x + 0.5) * log(tmp); 00042 double ser = 1.000000000190015; 00043 00044 for ( j = 0; j <= 5; j++ ) { 00045 x += 1.0; 00046 ser += cof[j]/x; 00047 } 00048 return -tmp + log(2.5066282746310005*ser); 00049 } 00050 00051 double gammln2(double xx) { 00052 00053 // Returns the value ln(Gamma(xx) for xx > 0. Full accuracy is obtained for 00054 // xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first. 00055 // (Adapted from Numerical Recipes in C) 00056 00057 static double cof[6] = {76.18009172947146,-86.50532032941677, 00058 24.01409824083091, -1.231739572450155, 00059 0.1208650973866179e-2, -0.5395239384953e-5}; 00060 int j; 00061 double x = xx - 0.0; 00062 double tmp = x + 5.5; 00063 tmp -= (x + 0.5) * log(tmp); 00064 double ser = 1.000000000190015; 00065 00066 for ( j = 0; j <= 5; j++ ) { 00067 x += 1.0; 00068 ser += cof[j]/x; 00069 } 00070 return -tmp + log(2.5066282746310005*ser/xx); 00071 } 00072 #include <iomanip> 00073 00074 00075 00076 int main() { 00077 00078 cout << "Enter 1 for RandGauss, 2 for RandGaussQ, 3 for DualRand flat: "; 00079 int choice; 00080 cin >> choice; 00081 00082 if (choice==1) { 00083 cout << "\n--------------------------------------------\n"; 00084 cout << "Test of Gauss distribution speed\n\n"; 00085 00086 long seed; 00087 cout << "Please enter an integer seed: "; 00088 cin >> seed; 00089 00090 int nNumbers; 00091 cout << "How many numbers should we generate: "; 00092 cin >> nNumbers; 00093 00094 cout << "\nInstantiating distribution utilizing DualRand engine...\n"; 00095 DualRand eng (seed); 00096 RandGauss dist (eng); 00097 00098 double sum = 0; 00099 00100 00101 for (int i=0; i < nNumbers; i++) { 00102 sum += dist.fire(); 00103 } 00104 00105 cout << "\n Finished: sum is " << sum << " \n"; 00106 } 00107 00108 if (choice==2) { 00109 cout << "\n--------------------------------------------\n"; 00110 cout << "Test of RandGaussQ distribution speed\n\n"; 00111 00112 long seed; 00113 cout << "Please enter an integer seed: "; 00114 cin >> seed; 00115 00116 int nNumbers; 00117 cout << "How many numbers should we generate: "; 00118 cin >> nNumbers; 00119 00120 cout << "\nInstantiating distribution utilizing DualRand engine...\n"; 00121 DualRand eng (seed); 00122 RandGaussQ dist (eng); 00123 00124 double sum = 0; 00125 00126 00127 for (int i=0; i < nNumbers; i++) { 00128 sum += dist.fire(); 00129 } 00130 00131 cout << "\n Finished: sum is " << sum << " \n"; 00132 } 00133 00134 if (choice==3) { 00135 cout << "\n--------------------------------------------\n"; 00136 cout << "Test of DualRand flat speed\n\n"; 00137 00138 long seed; 00139 cout << "Please enter an integer seed: "; 00140 cin >> seed; 00141 00142 int nNumbers; 00143 cout << "How many numbers should we generate: "; 00144 cin >> nNumbers; 00145 00146 cout << "\nInstantiating distribution utilizing DualRand engine...\n"; 00147 DualRand eng (seed); 00148 00149 double sum = 0; 00150 00151 00152 for (int i=0; i < nNumbers; i++) { 00153 sum += eng.flat(); 00154 } 00155 00156 cout << "\n Finished: sum is " << sum << " \n"; 00157 } 00158 00159 00160 #ifdef GAMMA 00161 cout << "\nNow we will compute the first 20 gammas, using gammln:\n"; 00162 00163 double x; 00164 for (x=1; x <= 20; x+=1) { 00165 cout << x << std::setprecision(20) << " " << exp(gammln1(x)) 00166 << " " << exp(gammln2(x)) << " difference in gammln2 = " << 00167 gammln1(x)-gammln2(x) << "\n"; 00168 } 00169 00170 00171 cout << "\nNow we will compute gamma of small numbers: \n"; 00172 00173 for ( x=1; x > .000000001; x *= .9 ) { 00174 cout << x << std::setprecision(20) << " " << 00175 1 - exp(gammln1(x)) * exp(gammln1(2-x)) * sin(CLHEP::pi*(1-x)) / (CLHEP::pi*(1-x)) << 00176 " " << 00177 1 - exp(gammln2(x)) * exp(gammln1(2-x)) * sin(CLHEP::pi*(1-x)) / (CLHEP::pi*(1-x)) << 00178 "\n"; 00179 } 00180 #endif // GAMMA 00181 00182 #ifdef POISSON 00183 cout << "\n--------------------------------------------\n"; 00184 cout << "Test of Poisson distribution speed\n\n"; 00185 00186 long seed; 00187 cout << "Please enter an integer seed: "; 00188 cin >> seed; 00189 00190 double mu; 00191 cout << "Please enter mu: "; 00192 cin >> mu; 00193 00194 int nNumbers; 00195 cout << "How many numbers should we generate: "; 00196 cin >> nNumbers; 00197 00198 cout << "\nInstantiating distribution utilizing DualRand engine...\n"; 00199 DualRand eng (seed); 00200 RandPoisson dist (eng, mu); 00201 // RandFlat dist (eng); 00202 00203 double sum = 0; 00204 00205 00206 for (int i=0; i < nNumbers; i++) { 00207 sum += dist.fire(); 00208 // sum += dist.quick(); 00209 // sum += dist.fire(mu); 00210 // sum += dist.quick(mu); 00211 00212 } 00213 00214 cout << "\n Finished: sum is " << sum << " \n"; 00215 #endif // POISSON 00216 00217 00218 #define MISC 00219 #ifdef MISC 00220 00221 DualRand e; 00222 00223 // RandGauss usage modes 00224 00225 cout << "testing RandGaussT::shoot(): " << RandGaussT::shoot() << "\n"; 00226 cout << "testing RandGaussT::shoot(&e): " << RandGaussT::shoot(&e) << "\n"; 00227 cout << "testing RandGaussT::shoot(100,10): " << 00228 RandGaussT::shoot(100,10) << "\n"; 00229 cout << "testing RandGaussT::shoot(&e,100,10): " << 00230 RandGaussT::shoot(&e,100,10) << "\n"; 00231 RandGaussT gt (e, 50,2); 00232 cout << "testing gt.fire(): " << gt.fire() << "\n"; 00233 cout << "testing gt.fire(200,2): " << gt.fire(200,2) << "\n"; 00234 00235 cout << "testing RandGaussQ::shoot(): " << RandGaussQ::shoot() << "\n"; 00236 cout << "testing RandGaussQ::shoot(&e): " << RandGaussQ::shoot(&e) << "\n"; 00237 cout << "testing RandGaussQ::shoot(100,10): " << 00238 RandGaussQ::shoot(100,10) << "\n"; 00239 cout << "testing RandGaussQ::shoot(&e,100,10): " << 00240 RandGaussQ::shoot(&e,100,10) << "\n"; 00241 RandGaussQ qt (e, 50,2); 00242 cout << "testing qt.fire(): " << qt.fire() << "\n"; 00243 cout << "testing qt.fire(200,2): " << qt.fire(200,2) << "\n"; 00244 00245 // RandPoisson usage modes 00246 00247 cout << "testing RandPoissonT::shoot(): " << RandPoissonT::shoot() << "\n"; 00248 cout << "testing RandPoissonT::shoot(&e): " 00249 << RandPoissonT::shoot(&e) << "\n"; 00250 cout << "testing RandPoissonT::shoot(90): " << 00251 RandPoissonT::shoot(90) << "\n"; 00252 cout << "testing RandPoissonT::shoot(&e,90): " << 00253 RandPoissonT::shoot(&e,90) << "\n"; 00254 RandPoissonT pgt (e,50); 00255 cout << "testing pgt.fire(): " << pgt.fire() << "\n"; 00256 cout << "testing pgt.fire(20): " << pgt.fire(20) << "\n"; 00257 00258 cout << "testing RandPoissonQ::shoot(): " << RandPoissonQ::shoot() << "\n"; 00259 cout << "testing RandPoissonQ::shoot(&e): " << RandPoissonQ::shoot(&e) << "\n"; 00260 cout << "testing RandPoissonQ::shoot(90): " << 00261 RandPoissonQ::shoot(90) << "\n"; 00262 cout << "testing RandPoissonQ::shoot(&e,90): " << 00263 RandPoissonQ::shoot(&e,90) << "\n"; 00264 RandPoissonQ pqt (e,50); 00265 cout << "testing pqt.fire(): " << pqt.fire() << "\n"; 00266 cout << "testing pqt.fire(20): " << pqt.fire(20) << "\n"; 00267 00268 // RandBit modes coming from RandFlat and bit modes 00269 00270 cout << "testing RandBit::shoot(): " << RandBit::shoot() << "\n"; 00271 cout << "testing RandBit::shoot(&e): " << RandBit::shoot(&e) << "\n"; 00272 cout << "testing RandBit::shoot(10,20): " << RandBit::shoot(10,20) << "\n"; 00273 cout << "testing RandBit::shoot(&e,10,20): "<< 00274 RandBit::shoot(&e,10,20) << "\n"; 00275 RandBit b ( e, 1000, 1100 ); 00276 cout << "testing b.fire(): " << b.fire() << "\n"; 00277 cout << "testing b.fire(10,20): " << b.fire(10,20) << "\n"; 00278 int i; 00279 cout << "testing RandBit::shootBit(): "; 00280 for (i=0; i<40; i++) { 00281 cout << RandBit::shootBit(); 00282 } cout << "\n"; 00283 cout << "testing RandBit::shootBit(&e): "; 00284 for (i=0; i<40; i++) { 00285 cout << RandBit::shootBit(&e); 00286 } cout << "\n"; 00287 cout << "testing RandBit::fireBit(): "; 00288 for (i=0; i<40; i++) { 00289 cout << b.fireBit(); 00290 } cout << "\n"; 00291 00292 // Timing for RandBit: 00293 00294 cout << "Timing RandFlat::shootBit(): Enter N: "; 00295 int N; 00296 cin >> N; 00297 int sum=0; 00298 for (i=0; i<N; i++) { 00299 sum+= RandFlat::shootBit(); 00300 } 00301 cout << "--------- Done.............. Sum = " << sum << "\n"; 00302 cout << "Timing RandBit::shootBit(): Enter N: "; 00303 cin >> N; 00304 sum = 0; 00305 for (i=0; i<N; i++) { 00306 sum += RandBit::shootBit(); 00307 } 00308 cout << "--------- Done.............. Sum = " << sum << "\n"; 00309 00310 #endif // MISC 00311 00312 return 0; 00313 } 00314