gaussSpeed.cc

Go to the documentation of this file.

#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"
 
using std::cin;
using std::cout;
using std::cerr;
using std::endl;
using namespace CLHEP;
//#ifndef _WIN32
//using std::exp;
//#endif
 
␌
// ---------
// RandGauss
//
// mf 12/13/04  Correction in way engines are supplied to RandBit() ctor
//          (gcc3.3.1 exposed previously innocuous mistake)
// ---------
 
double gammln1(double xx) {
 
// Returns the value ln(Gamma(xx) for xx > 0.  Full accuracy is obtained for 
// xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first.
// (Adapted from Numerical Recipes in C)
 
  static double cof[6] = {76.18009172947146,-86.50532032941677,
                             24.01409824083091, -1.231739572450155,
                             0.1208650973866179e-2, -0.5395239384953e-5};
  int j;
  double x = xx - 1.0;
  double tmp = x + 5.5;
  tmp -= (x + 0.5) * std::log(tmp);
  double ser = 1.000000000190015;
 
  for ( j = 0; j <= 5; j++ ) {
    x += 1.0;
    ser += cof[j]/x;
  }
  return -tmp + std::log(2.5066282746310005*ser);
}
 
double gammln2(double xx) {
 
// Returns the value ln(Gamma(xx) for xx > 0.  Full accuracy is obtained for 
// xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first.
// (Adapted from Numerical Recipes in C)
 
  static double cof[6] = {76.18009172947146,-86.50532032941677,
                             24.01409824083091, -1.231739572450155,
                             0.1208650973866179e-2, -0.5395239384953e-5};
  int j;
  double x = xx - 0.0;
  double tmp = x + 5.5;
  tmp -= (x + 0.5) * std::log(tmp);
  double ser = 1.000000000190015;
 
  for ( j = 0; j <= 5; j++ ) {
    x += 1.0;
    ser += cof[j]/x;
  }
  return -tmp + std::log(2.5066282746310005*ser/xx);
}
#include <iomanip>
 
 
␌
int main() {
 
  cout << "Enter 1 for RandGauss, 2 for RandGaussQ, 3 for DualRand flat: ";
  int choice;
  cin >> choice;
 
if (choice==1) {
  cout << "\n--------------------------------------------\n";
  cout << "Test of Gauss distribution speed\n\n";
 
  long seed;
  cout << "Please enter an integer seed:        ";
  cin >> seed;
 
  int nNumbers;
  cout << "How many numbers should we generate: ";
  cin >> nNumbers;
 
  cout << "\nInstantiating distribution utilizing DualRand engine...\n";
  DualRand eng (seed);
  RandGauss dist (eng);
 
  double sum = 0;
 
 
  for (int i=0; i < nNumbers; i++) {
    sum += dist.fire();
  }
 
  cout << "\n Finished:  sum is " << sum << " \n";
}
 
if (choice==2) {
  cout << "\n--------------------------------------------\n";
  cout << "Test of RandGaussQ distribution speed\n\n";
 
  long seed;
  cout << "Please enter an integer seed:        ";
  cin >> seed;
 
  int nNumbers;
  cout << "How many numbers should we generate: ";
  cin >> nNumbers;
 
  cout << "\nInstantiating distribution utilizing DualRand engine...\n";
  DualRand eng (seed);
  RandGaussQ dist (eng);
 
  double sum = 0;
 
 
  for (int i=0; i < nNumbers; i++) {
    sum += dist.fire();
  }
 
  cout << "\n Finished:  sum is " << sum << " \n";
}
 
if (choice==3) {
  cout << "\n--------------------------------------------\n";
  cout << "Test of DualRand flat speed\n\n";
 
  long seed;
  cout << "Please enter an integer seed:        ";
  cin >> seed;
 
  int nNumbers;
  cout << "How many numbers should we generate: ";
  cin >> nNumbers;
 
  cout << "\nInstantiating distribution utilizing DualRand engine...\n";
  DualRand eng (seed);
 
  double sum = 0;
 
 
  for (int i=0; i < nNumbers; i++) {
    sum += eng.flat();
  }
 
  cout << "\n Finished:  sum is " << sum << " \n";
}
␌
 
#ifdef GAMMA
  cout << "\nNow we will compute the first 20 gammas, using gammln:\n";
 
  double x;
  for (x=1; x <= 20; x+=1) {
    cout << x << std::setprecision(20) << "    " << std::exp(gammln1(x)) 
    << "    " << std::exp(gammln2(x)) << " difference in gammln2 = " << 
    gammln1(x)-gammln2(x) << "\n";
  }
 
 
  cout << "\nNow we will compute gamma of small numbers: \n";
 
  for ( x=1; x > .000000001; x *= .9 ) {
    cout << x << std::setprecision(20) << "    " <<
    1 - std::exp(gammln1(x)) * std::exp(gammln1(2-x)) * std::sin(CLHEP::pi*(1-x)) / (CLHEP::pi*(1-x)) <<
    "    " <<
    1 - std::exp(gammln2(x)) * std::exp(gammln1(2-x)) * std::sin(CLHEP::pi*(1-x)) / (CLHEP::pi*(1-x)) <<
    "\n";
  }
#endif // GAMMA
 
#ifdef POISSON
  cout << "\n--------------------------------------------\n";
  cout << "Test of Poisson distribution speed\n\n";
 
  long seed;
  cout << "Please enter an integer seed:        ";
  cin >> seed;
 
  double mu;
  cout << "Please enter mu:        ";
  cin >> mu;
 
  int nNumbers;
  cout << "How many numbers should we generate: ";
  cin >> nNumbers;
 
  cout << "\nInstantiating distribution utilizing DualRand engine...\n";
  DualRand eng (seed);
  RandPoisson dist (eng, mu);
  // RandFlat dist (eng);
 
  double sum = 0;
 
 
  for (int i=0; i < nNumbers; i++) {
    sum += dist.fire(); 
//    sum += dist.quick();
//    sum += dist.fire(mu);
//    sum += dist.quick(mu);
 
  }
 
  cout << "\n Finished:  sum is " << sum << " \n";
#endif // POISSON
 
␌
#define MISC
#ifdef MISC
 
  DualRand e;
 
  // RandGauss usage modes
 
  cout << "testing RandGaussT::shoot():  "  << RandGaussT::shoot()  << "\n";
  cout << "testing RandGaussT::shoot(&e): "  << RandGaussT::shoot(&e) << "\n";
  cout << "testing RandGaussT::shoot(100,10): " << 
            RandGaussT::shoot(100,10) << "\n";
  cout << "testing RandGaussT::shoot(&e,100,10): " << 
            RandGaussT::shoot(&e,100,10) << "\n";
  RandGaussT gt (e, 50,2);
  cout << "testing gt.fire(): " << gt.fire() << "\n";
  cout << "testing gt.fire(200,2): " << gt.fire(200,2) << "\n";
 
  cout << "testing RandGaussQ::shoot():  "  << RandGaussQ::shoot()  << "\n";
  cout << "testing RandGaussQ::shoot(&e): "  << RandGaussQ::shoot(&e) << "\n";
  cout << "testing RandGaussQ::shoot(100,10): " << 
            RandGaussQ::shoot(100,10) << "\n";
  cout << "testing RandGaussQ::shoot(&e,100,10): " << 
            RandGaussQ::shoot(&e,100,10) << "\n";
  RandGaussQ qt (e, 50,2);
  cout << "testing qt.fire(): " << qt.fire() << "\n";
  cout << "testing qt.fire(200,2): " << qt.fire(200,2) << "\n";
 
  // RandPoisson usage modes
 
  cout << "testing RandPoissonT::shoot():  "  << RandPoissonT::shoot()  << "\n";
  cout << "testing RandPoissonT::shoot(&e): "  
                    << RandPoissonT::shoot(&e) << "\n";
  cout << "testing RandPoissonT::shoot(90): " << 
            RandPoissonT::shoot(90) << "\n";
  cout << "testing RandPoissonT::shoot(&e,90): " << 
            RandPoissonT::shoot(&e,90) << "\n";
  RandPoissonT pgt (e,50);
  cout << "testing pgt.fire(): " << pgt.fire() << "\n";
  cout << "testing pgt.fire(20): " << pgt.fire(20) << "\n";
 
  cout << "testing RandPoissonQ::shoot():  "  << RandPoissonQ::shoot()  << "\n";
  cout << "testing RandPoissonQ::shoot(&e): " << RandPoissonQ::shoot(&e) << "\n";
  cout << "testing RandPoissonQ::shoot(90): " << 
            RandPoissonQ::shoot(90) << "\n";
  cout << "testing RandPoissonQ::shoot(&e,90): " << 
            RandPoissonQ::shoot(&e,90) << "\n";
  RandPoissonQ pqt (e,50);
  cout << "testing pqt.fire(): " << pqt.fire() << "\n";
  cout << "testing pqt.fire(20): " << pqt.fire(20) << "\n";
 
  // RandBit modes coming from RandFlat and bit modes
 
  cout << "testing RandBit::shoot():  " << RandBit::shoot()  << "\n";
  cout << "testing RandBit::shoot(&e): " << RandBit::shoot(&e) << "\n";
  cout << "testing RandBit::shoot(10,20): "  << RandBit::shoot(10,20)   << "\n";
  cout << "testing RandBit::shoot(&e,10,20): "<< 
                    RandBit::shoot(&e,10,20) << "\n";
  RandBit b ( e, 1000, 1100 );
  cout << "testing b.fire(): " << b.fire() << "\n";
  cout << "testing b.fire(10,20): " << b.fire(10,20) << "\n";
  int i;
  cout << "testing RandBit::shootBit():  "; 
  for (i=0; i<40; i++) {
    cout << RandBit::shootBit();
  } cout << "\n";
  cout << "testing RandBit::shootBit(&e):  "; 
  for (i=0; i<40; i++) {
    cout << RandBit::shootBit(&e);
  } cout << "\n";
  cout << "testing RandBit::fireBit():  ";
  for (i=0; i<40; i++) {
    cout << b.fireBit();
  } cout << "\n";
 
  // Timing for RandBit:
 
  cout << "Timing RandFlat::shootBit():  Enter N: ";
  int N;
  cin >> N;
  int sum=0;
  for (i=0; i<N; i++) {
    sum+= RandFlat::shootBit();
  }
  cout << "--------- Done.............. Sum = " << sum <<  "\n";
  cout << "Timing RandBit::shootBit():  Enter N: ";
  cin >>  N;
  sum = 0;
  for (i=0; i<N; i++) {
    sum += RandBit::shootBit();
  }
  cout << "--------- Done.............. Sum = " << sum <<  "\n";
 
#endif // MISC
 
  return 0;
}