CLHEP/Random/RandGaussZiggurat.h

0001 /*
0002 Code adapted from:
0003 http://www.jstatsoft.org/v05/i08/
0004
0005
0006 Original disclaimer:
0007
0008 The ziggurat method for RNOR and REXP
0009 Combine the code below with the main program in which you want
0010 normal or exponential variates.   Then use of RNOR in any expression
0011 will provide a standard normal variate with mean zero, variance 1,
0012 while use of REXP in any expression will provide an exponential variate
0013 with density exp(-x),x>0.
0014 Before using RNOR or REXP in your main, insert a command such as
0015 zigset(86947731 );
0016 with your own choice of seed value>0, rather than 86947731.
0017 (If you do not invoke zigset(...) you will get all zeros for RNOR and REXP.)
0018 For details of the method, see Marsaglia and Tsang, "The ziggurat method
0019 for generating random variables", Journ. Statistical Software.
0020
0021 */
0022
0023
0024 #ifndef RandGaussZiggurat_h
0025 #define RandGaussZiggurat_h 1
0026
0027 #include "cmath"
0028 #include "CLHEP/Random/defs.h"
0029 #include "CLHEP/Random/RandGauss.h"
0030 #include "CLHEP/Utility/thread_local.h"
0031
0032 namespace CLHEP {
0033
0034 /**
0035  * @author ATLAS
0036  * @ingroup random
0037  */
0038 class RandGaussZiggurat : public RandGauss {
0039
0040 public:
0041
0042   inline RandGaussZiggurat ( HepRandomEngine& anEngine, double mean=0.0, double stdDev=1.0 );
0043   inline RandGaussZiggurat ( HepRandomEngine* anEngine, double mean=0.0, double stdDev=1.0 );
0044
0045   // Destructor
0046   virtual ~RandGaussZiggurat();
0047
0048   // Static methods to shoot random values using the static generator
0049
0050   static inline float shoot() {return ziggurat_RNOR(HepRandom::getTheEngine());};
0051   static inline float shoot( float mean, float stdDev ) {return shoot()*stdDev + mean;};
0052
0053   static void shootArray ( const int size, float* vect, float mean=0.0, float stdDev=1.0 );
0054   static void shootArray ( const int size, double* vect, double mean=0.0, double stdDev=1.0 );
0055
0056   //  Static methods to shoot random values using a given engine
0057   //  by-passing the static generator.
0058
0059   static inline float shoot( HepRandomEngine* anotherEngine ) {return ziggurat_RNOR(anotherEngine);};
0060   static inline float shoot( HepRandomEngine* anotherEngine, float mean, float stdDev ) {return shoot(anotherEngine)*stdDev + mean;};
0061
0062   static void shootArray ( HepRandomEngine* anotherEngine, const int size, float* vect, float mean=0.0, float stdDev=1.0 );
0063   static void shootArray ( HepRandomEngine* anotherEngine, const int size, double* vect, double mean=0.0, double stdDev=1.0 );
0064
0065   //  Instance methods using the localEngine to instead of the static
0066   //  generator, and the default mean and stdDev established at construction
0067
0068   inline float fire() {return ziggurat_RNOR(localEngine.get()) * defaultStdDev + defaultMean;};
0069
0070   inline float fire ( float mean, float stdDev ) {return ziggurat_RNOR(localEngine.get()) * stdDev + mean;};
0071
0072   void fireArray  ( const int size, float* vect);
0073   void fireArray  ( const int size, double* vect);
0074   void fireArray  ( const int size, float* vect, float mean, float stdDev );
0075   void fireArray  ( const int size, double* vect, double mean, double stdDev );
0076
0077   virtual double operator()();
0078   virtual double operator()( double mean, double stdDev );
0079
0080   // Save and restore to/from streams
0081
0082   std::ostream & put ( std::ostream & os ) const;
0083   std::istream & get ( std::istream & is );
0084
0085   std::string name() const;
0086   HepRandomEngine & engine();
0087
0088   static std::string distributionName() {return "RandGaussZiggurat";}
0089   // Provides the name of this distribution class
0090
0091   static bool ziggurat_init();
0092 protected:
0093   //////////////////////////
0094   // Ziggurat Original code:
0095   //////////////////////////
0096   //static unsigned long jz,jsr=123456789;
0097   //
0098   //#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr)
0099   //#define UNI (.5 + (signed) SHR3*.2328306e-9)
0100   //#define IUNI SHR3
0101   //
0102   //static long hz;
0103   //static unsigned long iz, kn[128], ke[256];
0104   //static float wn[128],fn[128], we[256],fe[256];
0105   //
0106   //#define RNOR (hz=SHR3, iz=hz&127, (std::fabs(hz)<kn[iz])? hz*wn[iz] : nfix())
0107   //#define REXP (jz=SHR3, iz=jz&255, (    jz <ke[iz])? jz*we[iz] : efix())
0108
0109   static CLHEP_THREAD_LOCAL unsigned long kn[128], ke[256];
0110   static CLHEP_THREAD_LOCAL float wn[128],fn[128], we[256],fe[256];
0111
0112   static CLHEP_THREAD_LOCAL bool ziggurat_is_init;
0113
0114   static inline unsigned long ziggurat_SHR3(HepRandomEngine* anEngine) {return (unsigned int)(*anEngine);};
0115   static inline float ziggurat_UNI(HepRandomEngine* anEngine) {return anEngine->flat();};
0116   static inline float ziggurat_RNOR(HepRandomEngine* anEngine) {
0117     if(!ziggurat_is_init) ziggurat_init();
0118     long hz=(signed)ziggurat_SHR3(anEngine);
0119     unsigned long iz=hz&127;
0120     return ((unsigned long)std::abs(hz)<kn[iz]) ? hz*wn[iz] : ziggurat_nfix(hz,anEngine);
0121   };
0122   static float ziggurat_nfix(long hz,HepRandomEngine* anEngine);
0123
0124 private:
0125
0126   // Private copy constructor. Defining it here disallows use.
0127   RandGaussZiggurat(const RandGaussZiggurat& d);
0128
0129   // All the engine info, and the default mean and sigma, are in the RandGauss
0130   // base class.
0131
0132 };
0133
0134 }  // namespace CLHEP
0135
0136 #ifdef ENABLE_BACKWARDS_COMPATIBILITY
0137 //  backwards compatibility will be enabled ONLY in CLHEP 1.9
0138 using namespace CLHEP;
0139 #endif
0140
0141 namespace CLHEP {
0142
0143 RandGaussZiggurat::RandGaussZiggurat(HepRandomEngine & anEngine, double mean,double stdDev ): RandGauss(anEngine, mean, stdDev)
0144 {
0145 }
0146
0147 RandGaussZiggurat::RandGaussZiggurat(HepRandomEngine * anEngine, double mean,double stdDev ): RandGauss(anEngine, mean, stdDev)
0148 {
0149 }
0150
0151 }  // namespace CLHEP
0152
0153 #endif