including RandGaussZiggurat and RandExpZiggurat

ChangeLog

 
+2013-09-26 Lynn Garren <garren@fnal.gov> 
+
+	* Random: including RandGaussZiggurat and RandExpZiggurat
+
 2013-01-11 Lynn Garren <garren@fnal.gov> 
 
 	* Matrix: change the names of more internal variables so -Wshadow does not complain

Random/ChangeLog

+
+2013-09-26 Lynn Garren <garren@fnal.gov> 
+
+	* Random: including RandGaussZiggurat and RandExpZiggurat
+	  (code supplied by ATLAS)
+
 ==============================
 14.11.12 Release CLHEP-2.1.3.1
 ==============================

Random/Random/CMakeLists.txt

@@ -20,6 +20,7 @@ set( pkginclude_HEADERS DoubConv.hh
 			RandEngine.h
 			RandExponential.h
 			RandExponential.icc
+			RandExpZiggurat.h
 			RandFlat.h
 			RandFlat.icc
 			RandGamma.h
@@ -30,6 +31,7 @@ set( pkginclude_HEADERS DoubConv.hh
 			RandGaussQ.icc
 			RandGaussT.h
 			RandGaussT.icc
+			RandGaussZiggurat.h
 			RandGeneral.h
 			RandGeneral.icc
 			RandLandau.h

Random/Random/Makefile.am

@@ -30,6 +30,7 @@ pkginclude_HEADERS = \
 	RandEngine.h	\
 	RandExponential.h	\
 	RandExponential.icc	\
+	RandExpZiggurat.h	\
 	RandFlat.h	\
 	RandFlat.icc	\
 	RandGamma.h	\
@@ -40,6 +41,7 @@ pkginclude_HEADERS = \
 	RandGaussQ.icc	\
 	RandGaussT.h	\
 	RandGaussT.icc	\
+	RandGaussZiggurat.h	\
 	RandGeneral.h	\
 	RandGeneral.icc	\
 	RandLandau.h	\

Random/Random/RandExpZiggurat.h

+/*
+Code adapted from:
+http://www.jstatsoft.org/v05/i08/
+
+
+Original disclaimer:
+
+The ziggurat method for RNOR and REXP
+Combine the code below with the main program in which you want
+normal or exponential variates.   Then use of RNOR in any expression
+will provide a standard normal variate with mean zero, variance 1,
+while use of REXP in any expression will provide an exponential variate
+with density exp(-x),x>0.
+Before using RNOR or REXP in your main, insert a command such as
+zigset(86947731 );
+with your own choice of seed value>0, rather than 86947731.
+(If you do not invoke zigset(...) you will get all zeros for RNOR and REXP.)
+For details of the method, see Marsaglia and Tsang, "The ziggurat method
+for generating random variables", Journ. Statistical Software.
+
+*/
+
+#ifndef RandExpZiggurat_h
+#define RandExpZiggurat_h 1
+
+#include "CLHEP/Random/defs.h"
+#include "CLHEP/Random/Random.h"
+
+namespace CLHEP {
+
+/**
+ * @author <Gabriele.Cosmo@cern.ch>
+ * @ingroup random
+ */
+class RandExpZiggurat : public HepRandom {
+
+public:
+
+  inline RandExpZiggurat ( HepRandomEngine& anEngine, double mean=1.0 );
+  inline RandExpZiggurat ( HepRandomEngine* anEngine, double mean=1.0 );
+  // These constructors should be used to instantiate a RandExpZiggurat
+  // distribution object defining a local engine for it.
+  // The static generator will be skipped using the non-static methods
+  // defined below.
+  // If the engine is passed by pointer the corresponding engine object
+  // will be deleted by the RandExpZiggurat destructor.
+  // If the engine is passed by reference the corresponding engine object
+  // will not be deleted by the RandExpZiggurat destructor.
+
+  virtual ~RandExpZiggurat();
+  // Destructor
+
+  // Static methods to shoot random values using the static generator
+
+  static float shoot() {return shoot(HepRandom::getTheEngine());};
+  static float shoot( float mean ) {return shoot(HepRandom::getTheEngine(),mean);};
+
+  /* ENGINE IS INTRINSIC FLOAT
+  static double shoot() {return shoot(HepRandom::getTheEngine());};
+  static double shoot( double mean ) {return shoot(HepRandom::getTheEngine(),mean);};
+  */
+
+  static void shootArray ( const int size, float* vect, float mean=1.0 );
+  static void shootArray ( const int size, double* vect, double mean=1.0 );
+
+  //  Static methods to shoot random values using a given engine
+  //  by-passing the static generator.
+
+  static inline float shoot( HepRandomEngine* anEngine ) {return ziggurat_REXP(anEngine);};
+  static inline float shoot( HepRandomEngine* anEngine, float mean ) {return shoot(anEngine)*mean;};
+  
+  /* ENGINE IS INTRINSIC FLOAT
+  static inline double shoot( HepRandomEngine* anEngine ) {return ziggurat_REXP(anEngine);};
+
+  static inline double shoot( HepRandomEngine* anEngine, double mean ) {return shoot(anEngine)*mean;};
+  */
+  
+  static void shootArray ( HepRandomEngine* anEngine, const int size, float* vect, float mean=1.0 );
+  static void shootArray ( HepRandomEngine* anEngine, const int size, double* vect, double mean=1.0 );
+
+  //  Methods using the localEngine to shoot random values, by-passing
+  //  the static generator.
+
+  inline float fire() {return fire(defaultMean);};
+  inline float fire( float mean ) {return ziggurat_REXP(localEngine)*mean;};
+  
+  /* ENGINE IS INTRINSIC FLOAT
+  inline double fire() {return fire(defaultMean);};
+  inline double fire( double mean ) {return ziggurat_REXP(localEngine)*mean;};
+  */
+  
+  void fireArray ( const int size, float* vect );
+  void fireArray ( const int size, double* vect );
+  void fireArray ( const int size, float* vect, float mean );
+  void fireArray ( const int size, double* vect, double mean );
+  
+  virtual double operator()();
+  inline float operator()( float mean ) {return fire( mean );};
+
+  // Save and restore to/from streams
+  
+  std::ostream & put ( std::ostream & os ) const;
+  std::istream & get ( std::istream & is );
+
+  std::string name() const;
+  HepRandomEngine & engine();
+
+  static std::string distributionName() {return "RandExpZiggurat";}  
+  // Provides the name of this distribution class
+  
+  static bool ziggurat_init();
+protected:
+  //////////////////////////
+  // Ziggurat Original code:
+  //////////////////////////
+  //static unsigned long jz,jsr=123456789;
+  //
+  //#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr)
+  //#define UNI (.5 + (signed) SHR3*.2328306e-9)
+  //#define IUNI SHR3
+  //
+  //static long hz;
+  //static unsigned long iz, kn[128], ke[256];
+  //static float wn[128],fn[128], we[256],fe[256];
+  //
+  //#define RNOR (hz=SHR3, iz=hz&127, (fabs(hz)<kn[iz])? hz*wn[iz] : nfix())
+  //#define REXP (jz=SHR3, iz=jz&255, (    jz <ke[iz])? jz*we[iz] : efix())
+
+  static unsigned long kn[128], ke[256];
+  static float wn[128],fn[128], we[256],fe[256];
+
+  static bool ziggurat_is_init;
+
+  static inline unsigned long ziggurat_SHR3(HepRandomEngine* anEngine) {return (unsigned int)(*anEngine);};
+  static inline float ziggurat_UNI(HepRandomEngine* anEngine) {return anEngine->flat();};
+  static inline float ziggurat_REXP(HepRandomEngine* anEngine) {
+    unsigned long jz=ziggurat_SHR3(anEngine);
+    unsigned long iz=jz&255;
+    return (jz<ke[iz]) ? jz*we[iz] : ziggurat_efix(jz,anEngine);
+  };
+  static float ziggurat_efix(unsigned long jz,HepRandomEngine* anEngine);
+
+private:
+
+  // Private copy constructor. Defining it here disallows use.
+  RandExpZiggurat(const RandExpZiggurat& d);
+
+  HepRandomEngine* localEngine;
+  bool deleteEngine;
+  double defaultMean;
+};
+
+}  // namespace CLHEP
+
+#ifdef ENABLE_BACKWARDS_COMPATIBILITY
+//  backwards compatibility will be enabled ONLY in CLHEP 1.9
+using namespace CLHEP;
+#endif
+
+namespace CLHEP {
+
+inline RandExpZiggurat::RandExpZiggurat(HepRandomEngine & anEngine, double mean ) : localEngine(&anEngine), deleteEngine(false), defaultMean(mean) 
+{
+  if(!ziggurat_is_init) ziggurat_init();
+}
+
+inline RandExpZiggurat::RandExpZiggurat(HepRandomEngine * anEngine, double mean ) : localEngine(anEngine), deleteEngine(true), defaultMean(mean) 
+{
+  if(!ziggurat_is_init) ziggurat_init();
+}
+
+}  // namespace CLHEP
+
+#endif

Random/Random/RandGaussZiggurat.h

+/*
+Code adapted from:
+http://www.jstatsoft.org/v05/i08/
+
+
+Original disclaimer:
+
+The ziggurat method for RNOR and REXP
+Combine the code below with the main program in which you want
+normal or exponential variates.   Then use of RNOR in any expression
+will provide a standard normal variate with mean zero, variance 1,
+while use of REXP in any expression will provide an exponential variate
+with density exp(-x),x>0.
+Before using RNOR or REXP in your main, insert a command such as
+zigset(86947731 );
+with your own choice of seed value>0, rather than 86947731.
+(If you do not invoke zigset(...) you will get all zeros for RNOR and REXP.)
+For details of the method, see Marsaglia and Tsang, "The ziggurat method
+for generating random variables", Journ. Statistical Software.
+
+*/
+
+
+#ifndef RandGaussZiggurat_h
+#define RandGaussZiggurat_h 1
+
+#include "math.h"
+#include "CLHEP/Random/defs.h"
+#include "CLHEP/Random/RandGauss.h"
+
+namespace CLHEP {
+
+/**
+ * @author ATLAS
+ * @ingroup random
+ */
+class RandGaussZiggurat : public RandGauss {
+
+public:
+
+  inline RandGaussZiggurat ( HepRandomEngine& anEngine, double mean=0.0, double stdDev=1.0 );
+  inline RandGaussZiggurat ( HepRandomEngine* anEngine, double mean=0.0, double stdDev=1.0 );
+
+  // Destructor
+  virtual ~RandGaussZiggurat();
+
+  // Static methods to shoot random values using the static generator
+  
+  static inline float shoot() {return ziggurat_RNOR(HepRandom::getTheEngine());};
+  static inline float shoot( float mean, float stdDev ) {return shoot()*stdDev + mean;};
+
+  static void shootArray ( const int size, float* vect, float mean=0.0, float stdDev=1.0 );
+  static void shootArray ( const int size, double* vect, double mean=0.0, double stdDev=1.0 );
+
+  //  Static methods to shoot random values using a given engine
+  //  by-passing the static generator.
+
+  static inline float shoot( HepRandomEngine* anotherEngine ) {return ziggurat_RNOR(anotherEngine);};
+  static inline float shoot( HepRandomEngine* anotherEngine, float mean, float stdDev ) {return shoot(anotherEngine)*stdDev + mean;};
+
+  static void shootArray ( HepRandomEngine* anotherEngine, const int size, float* vect, float mean=0.0, float stdDev=1.0 );
+  static void shootArray ( HepRandomEngine* anotherEngine, const int size, double* vect, double mean=0.0, double stdDev=1.0 );
+
+  //  Instance methods using the localEngine to instead of the static 
+  //  generator, and the default mean and stdDev established at construction
+
+  inline float fire() {return ziggurat_RNOR(localEngine.get()) * defaultStdDev + defaultMean;};
+
+  inline float fire ( float mean, float stdDev ) {return ziggurat_RNOR(localEngine.get()) * stdDev + mean;};
+  
+  void fireArray  ( const int size, float* vect);
+  void fireArray  ( const int size, double* vect);
+  void fireArray  ( const int size, float* vect, float mean, float stdDev );
+  void fireArray  ( const int size, double* vect, double mean, double stdDev );
+
+  virtual double operator()();
+  virtual double operator()( double mean, double stdDev );
+
+  // Save and restore to/from streams
+  
+  std::ostream & put ( std::ostream & os ) const;
+  std::istream & get ( std::istream & is );
+
+  std::string name() const;
+  HepRandomEngine & engine();
+
+  static std::string distributionName() {return "RandGaussZiggurat";}  
+  // Provides the name of this distribution class
+  
+  static bool ziggurat_init();
+protected:  
+  //////////////////////////
+  // Ziggurat Original code:
+  //////////////////////////
+  //static unsigned long jz,jsr=123456789;
+  //
+  //#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr)
+  //#define UNI (.5 + (signed) SHR3*.2328306e-9)
+  //#define IUNI SHR3
+  //
+  //static long hz;
+  //static unsigned long iz, kn[128], ke[256];
+  //static float wn[128],fn[128], we[256],fe[256];
+  //
+  //#define RNOR (hz=SHR3, iz=hz&127, (fabs(hz)<kn[iz])? hz*wn[iz] : nfix())
+  //#define REXP (jz=SHR3, iz=jz&255, (    jz <ke[iz])? jz*we[iz] : efix())
+
+  static unsigned long kn[128], ke[256];
+  static float wn[128],fn[128], we[256],fe[256];
+
+  static bool ziggurat_is_init;
+  
+  static inline unsigned long ziggurat_SHR3(HepRandomEngine* anEngine) {return (unsigned int)(*anEngine);};
+  static inline float ziggurat_UNI(HepRandomEngine* anEngine) {return anEngine->flat();};
+  static inline float ziggurat_RNOR(HepRandomEngine* anEngine) {
+    long hz=(signed)ziggurat_SHR3(anEngine);
+    unsigned long iz=hz&127;
+    return (fabs(hz)<kn[iz]) ? hz*wn[iz] : ziggurat_nfix(hz,anEngine);
+  };
+  static float ziggurat_nfix(long hz,HepRandomEngine* anEngine);
+  
+private:
+
+  // Private copy constructor. Defining it here disallows use.
+  RandGaussZiggurat(const RandGaussZiggurat& d);
+
+  // All the engine info, and the default mean and sigma, are in the RandGauss
+  // base class.
+
+};
+
+}  // namespace CLHEP
+
+#ifdef ENABLE_BACKWARDS_COMPATIBILITY
+//  backwards compatibility will be enabled ONLY in CLHEP 1.9
+using namespace CLHEP;
+#endif
+
+namespace CLHEP {
+
+RandGaussZiggurat::RandGaussZiggurat(HepRandomEngine & anEngine, double mean,double stdDev ): RandGauss(anEngine, mean, stdDev) 
+{
+  if(!ziggurat_is_init) ziggurat_init();
+}
+
+RandGaussZiggurat::RandGaussZiggurat(HepRandomEngine * anEngine, double mean,double stdDev ): RandGauss(anEngine, mean, stdDev) 
+{
+  if(!ziggurat_is_init) ziggurat_init();
+}
+
+}  // namespace CLHEP
+
+#endif

Random/src/CMakeLists.txt

-# This include has been moved to the top CMakeLists.txt 
-# add CLHEP/Random to search path so we find gaussTables.cdat
-#include_directories ("${CMAKE_CURRENT_SOURCE_DIR}/..")
 
 clhep_build_library( Random DoubConv.cc
 			    DRand48Engine.cc
@@ -21,11 +18,13 @@ clhep_build_library( Random DoubConv.cc
 			    RandChiSquare.cc
 			    RandEngine.cc
 			    RandExponential.cc
+			    RandExpZiggurat.cc
 			    RandFlat.cc
 			    RandGamma.cc
 			    RandGauss.cc
 			    RandGaussQ.cc
 			    RandGaussT.cc
+			    RandGaussZiggurat.cc
 			    RandGeneral.cc
 			    RandLandau.cc
 			    Random.cc

Random/src/Makefile.am

@@ -26,11 +26,13 @@ libCLHEP_Random_@VERSION@_a_SOURCES = \
 	RandChiSquare.cc	\
 	RandEngine.cc	\
 	RandExponential.cc	\
+	RandExpZiggurat.cc	\
 	RandFlat.cc	\
 	RandGamma.cc	\
 	RandGauss.cc	\
 	RandGaussQ.cc	\
 	RandGaussT.cc	\
+	RandGaussZiggurat.cc	\
 	RandGeneral.cc	\
 	RandLandau.cc	\
 	Random.cc	\

Random/src/RandExpZiggurat.cc

+#include "CLHEP/Random/defs.h"
+#include "CLHEP/Random/DoubConv.hh"
+
+#include "CLHEP/Random/RandExpZiggurat.h"
+
+#include <iostream>
+#include <cmath>	// for log()
+
+namespace CLHEP {
+
+bool RandExpZiggurat::ziggurat_is_init=RandExpZiggurat::ziggurat_init();
+unsigned long RandExpZiggurat::kn[128], RandExpZiggurat::ke[256];
+float RandExpZiggurat::wn[128],RandExpZiggurat::fn[128],RandExpZiggurat::we[256],RandExpZiggurat::fe[256];
+
+std::string RandExpZiggurat::name() const {return "RandExpZiggurat";}
+
+HepRandomEngine & RandExpZiggurat::engine() {return *localEngine;}
+
+RandExpZiggurat::~RandExpZiggurat() {
+  if ( deleteEngine ) delete localEngine;
+}
+
+RandExpZiggurat::RandExpZiggurat(const RandExpZiggurat& right) : HepRandom(right),defaultMean(right.defaultMean)
+{
+  if(!ziggurat_is_init) ziggurat_init();
+}
+
+double RandExpZiggurat::operator()()
+{
+  return fire( defaultMean );
+}
+
+void RandExpZiggurat::shootArray( const int size, float* vect, float mean )
+{
+   for (int i=0; i<size; ++i) vect[i] = shoot(mean);
+}
+
+void RandExpZiggurat::shootArray( const int size, double* vect, double mean )
+{
+   for (int i=0; i<size; ++i) vect[i] = shoot(mean);
+}
+
+void RandExpZiggurat::shootArray(HepRandomEngine* anEngine, const int size, float* vect, float mean )
+{
+   for (int i=0; i<size; ++i) vect[i] = shoot(anEngine, mean);
+}
+
+void RandExpZiggurat::shootArray(HepRandomEngine* anEngine, const int size, double* vect, double mean )
+{
+   for (int i=0; i<size; ++i) vect[i] = shoot(anEngine, mean);
+}
+
+void RandExpZiggurat::fireArray( const int size, float* vect)
+{
+   for (int i=0; i<size; ++i) vect[i] = fire( defaultMean );
+}
+
+void RandExpZiggurat::fireArray( const int size, double* vect)
+{
+   for (int i=0; i<size; ++i) vect[i] = fire( defaultMean );
+}
+
+void RandExpZiggurat::fireArray( const int size, float* vect, float mean )
+{
+   for (int i=0; i<size; ++i) vect[i] = fire( mean );
+}
+
+void RandExpZiggurat::fireArray( const int size, double* vect, double mean )
+{
+   for (int i=0; i<size; ++i) vect[i] = fire( mean );
+}
+
+std::ostream & RandExpZiggurat::put ( std::ostream & os ) const {
+  int pr=os.precision(20);
+  std::vector<unsigned long> t(2);
+  os << " " << name() << "\n";
+  os << "Uvec" << "\n";
+  t = DoubConv::dto2longs(defaultMean);
+  os << defaultMean << " " << t[0] << " " << t[1] << "\n";
+  os.precision(pr);
+  return os;
+#ifdef REMOVED
+  int pr=os.precision(20);
+  os << " " << name() << "\n";
+  os << defaultMean << "\n";
+  os.precision(pr);
+  return os;
+#endif
+}
+
+std::istream & RandExpZiggurat::get ( std::istream & is ) {
+  std::string inName;
+  is >> inName;
+  if (inName != name()) {
+    is.clear(std::ios::badbit | is.rdstate());
+    std::cerr << "Mismatch when expecting to read state of a "
+    	      << name() << " distribution\n"
+	      << "Name found was " << inName
+	      << "\nistream is left in the badbit state\n";
+    return is;
+  }
+  if (possibleKeywordInput(is, "Uvec", defaultMean)) {
+    std::vector<unsigned long> t(2);
+    is >> defaultMean >> t[0] >> t[1]; defaultMean = DoubConv::longs2double(t); 
+    return is;
+  }
+  // is >> defaultMean encompassed by possibleKeywordInput
+  return is;
+}
+
+
+float RandExpZiggurat::ziggurat_efix(unsigned long jz,HepRandomEngine* anEngine)
+{ 
+  unsigned long iz=jz&255;
+  
+  float x;
+  for(;;)
+  {  
+    if(iz==0) return (7.69711-log(ziggurat_UNI(anEngine)));          /* iz==0 */
+    x=jz*we[iz]; 
+    if( fe[iz]+ziggurat_UNI(anEngine)*(fe[iz-1]-fe[iz]) < exp(-x) ) return (x);
+
+    /* initiate, try to exit for(;;) loop */
+    jz=ziggurat_SHR3(anEngine);
+    iz=(jz&255);
+    if(jz<ke[iz]) return (jz*we[iz]);
+  }
+}
+
+bool RandExpZiggurat::ziggurat_init()
+{  
+  const double m1 = 2147483648.0, m2 = 4294967296.;
+  double dn=3.442619855899,tn=dn,vn=9.91256303526217e-3, q;
+  double de=7.697117470131487, te=de, ve=3.949659822581572e-3;
+  int i;
+
+/* Set up tables for RNOR */
+  q=vn/exp(-.5*dn*dn);
+  kn[0]=(unsigned long)((dn/q)*m1);
+  kn[1]=0;
+
+  wn[0]=q/m1;
+  wn[127]=dn/m1;
+
+  fn[0]=1.;
+  fn[127]=exp(-.5*dn*dn);
+
+  for(i=126;i>=1;i--) {
+    dn=sqrt(-2.*log(vn/dn+exp(-.5*dn*dn)));
+    kn[i+1]=(unsigned long)((dn/tn)*m1);
+    tn=dn;
+    fn[i]=exp(-.5*dn*dn);
+    wn[i]=dn/m1;
+  }
+
+/* Set up tables for REXP */
+  q = ve/exp(-de);
+  ke[0]=(unsigned long)((de/q)*m2);
+  ke[1]=0;
+
+  we[0]=q/m2;
+  we[255]=de/m2;
+
+  fe[0]=1.;
+  fe[255]=exp(-de);
+
+  for(i=254;i>=1;i--) {
+    de=-log(ve/de+exp(-de));
+    ke[i+1]= (unsigned long)((de/te)*m2);
+    te=de;
+    fe[i]=exp(-de);
+    we[i]=de/m2;
+  }
+  ziggurat_is_init=true;
+  return true;
+}
+
+}  // namespace CLHEP

Random/src/RandGaussZiggurat.cc

+#include "CLHEP/Random/defs.h"
+#include "CLHEP/Random/RandGaussZiggurat.h"
+#include "CLHEP/Units/PhysicalConstants.h"
+#include <iostream>
+#include <cmath>	// for log()
+
+namespace CLHEP {
+
+bool RandGaussZiggurat::ziggurat_is_init=RandGaussZiggurat::ziggurat_init();
+//bool RandGaussZiggurat::ziggurat_is_init=false;
+unsigned long RandGaussZiggurat::kn[128], RandGaussZiggurat::ke[256];
+float RandGaussZiggurat::wn[128],RandGaussZiggurat::fn[128],RandGaussZiggurat::we[256],RandGaussZiggurat::fe[256];
+
+HepRandomEngine & RandGaussZiggurat::engine() {return RandGauss::engine();}