CLHEP/Random/RandExpZiggurat.h

0001 // $Id:$
0002 // -*- C++ -*-
0003 //
0004 // -----------------------------------------------------------------------
0005 //                             HEP Random
0006 //                       --- RandExpZiggurat ---
0007 //                          class header file
0008 // -----------------------------------------------------------------------
0009 // This file is part of Geant4 (simulation toolkit for HEP).
0010
0011 // Class defining methods for shooting or firing Landau distributed
0012 // random values.
0013
0014 // =======================================================================
0015 //
0016 // Code adapted from:
0017 // http://www.jstatsoft.org/v05/i08/
0018 //
0019 //
0020 // Original disclaimer:
0021
0022 // The ziggurat method for RNOR and REXP
0023 // Combine the code below with the main program in which you want
0024 // normal or exponential variates.   Then use of RNOR in any expression
0025 // will provide a standard normal variate with mean zero, variance 1,
0026 // while use of REXP in any expression will provide an exponential variate
0027 // with density exp(-x),x>0.
0028 // Before using RNOR or REXP in your main, insert a command such as
0029 // zigset(86947731 );
0030 // with your own choice of seed value>0, rather than 86947731.
0031 // (If you do not invoke zigset(...) you will get all zeros for RNOR and REXP.)
0032 // For details of the method, see Marsaglia and Tsang, "The ziggurat method
0033 // for generating random variables", Journ. Statistical Software.
0034
0035 // =======================================================================
0036
0037 #ifndef RandExpZiggurat_h
0038 #define RandExpZiggurat_h 1
0039
0040 #include "CLHEP/Random/defs.h"
0041 #include "CLHEP/Random/Random.h"
0042 #include "CLHEP/Utility/memory.h"
0043 #include "CLHEP/Utility/thread_local.h"
0044
0045 namespace CLHEP {
0046
0047 /**
0048  * @author ATLAS
0049  * @ingroup random
0050  */
0051 class RandExpZiggurat : public HepRandom {
0052
0053 public:
0054
0055   inline RandExpZiggurat ( HepRandomEngine& anEngine, double mean=1.0 );
0056   inline RandExpZiggurat ( HepRandomEngine* anEngine, double mean=1.0 );
0057   // These constructors should be used to instantiate a RandExpZiggurat
0058   // distribution object defining a local engine for it.
0059   // The static generator will be skipped using the non-static methods
0060   // defined below.
0061   // If the engine is passed by pointer the corresponding engine object
0062   // will be deleted by the RandExpZiggurat destructor.
0063   // If the engine is passed by reference the corresponding engine object
0064   // will not be deleted by the RandExpZiggurat destructor.
0065
0066   virtual ~RandExpZiggurat();
0067   // Destructor
0068
0069   // Static methods to shoot random values using the static generator
0070
0071   static float shoot() {return shoot(HepRandom::getTheEngine());};
0072   static float shoot( float mean ) {return shoot(HepRandom::getTheEngine(),mean);};
0073
0074   /* ENGINE IS INTRINSIC FLOAT
0075   static double shoot() {return shoot(HepRandom::getTheEngine());};
0076   static double shoot( double mean ) {return shoot(HepRandom::getTheEngine(),mean);};
0077   */
0078
0079   static void shootArray ( const int size, float* vect, float mean=1.0 );
0080   static void shootArray ( const int size, double* vect, double mean=1.0 );
0081
0082   //  Static methods to shoot random values using a given engine
0083   //  by-passing the static generator.
0084
0085   static inline float shoot( HepRandomEngine* anEngine ) {return ziggurat_REXP(anEngine);};
0086   static inline float shoot( HepRandomEngine* anEngine, float mean ) {return shoot(anEngine)*mean;};
0087
0088   /* ENGINE IS INTRINSIC FLOAT
0089   static inline double shoot( HepRandomEngine* anEngine ) {return ziggurat_REXP(anEngine);};
0090
0091   static inline double shoot( HepRandomEngine* anEngine, double mean ) {return shoot(anEngine)*mean;};
0092   */
0093
0094   static void shootArray ( HepRandomEngine* anEngine, const int size, float* vect, float mean=1.0 );
0095   static void shootArray ( HepRandomEngine* anEngine, const int size, double* vect, double mean=1.0 );
0096
0097   //  Methods using the localEngine to shoot random values, by-passing
0098   //  the static generator.
0099
0100   inline float fire() {return fire(defaultMean);};
0101   inline float fire( float mean ) {return ziggurat_REXP(localEngine.get())*mean;};
0102
0103   /* ENGINE IS INTRINSIC FLOAT
0104   inline double fire() {return fire(defaultMean);};
0105   inline double fire( double mean ) {return ziggurat_REXP(localEngine.get())*mean;};
0106   */
0107
0108   void fireArray ( const int size, float* vect );
0109   void fireArray ( const int size, double* vect );
0110   void fireArray ( const int size, float* vect, float mean );
0111   void fireArray ( const int size, double* vect, double mean );
0112
0113   virtual double operator()();
0114   inline float operator()( float mean ) {return fire( mean );};
0115
0116   // Save and restore to/from streams
0117
0118   std::ostream & put ( std::ostream & os ) const;
0119   std::istream & get ( std::istream & is );
0120
0121   std::string name() const;
0122   HepRandomEngine & engine();
0123
0124   static std::string distributionName() {return "RandExpZiggurat";}
0125   // Provides the name of this distribution class
0126
0127   static bool ziggurat_init();
0128 protected:
0129   //////////////////////////
0130   // Ziggurat Original code:
0131   //////////////////////////
0132   //static unsigned long jz,jsr=123456789;
0133   //
0134   //#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr)
0135   //#define UNI (.5 + (signed) SHR3*.2328306e-9)
0136   //#define IUNI SHR3
0137   //
0138   //static long hz;
0139   //static unsigned long iz, kn[128], ke[256];
0140   //static float wn[128],fn[128], we[256],fe[256];
0141   //
0142   //#define RNOR (hz=SHR3, iz=hz&127, (fabs(hz)<kn[iz])? hz*wn[iz] : nfix())
0143   //#define REXP (jz=SHR3, iz=jz&255, (    jz <ke[iz])? jz*we[iz] : efix())
0144
0145   static CLHEP_THREAD_LOCAL unsigned long kn[128], ke[256];
0146   static CLHEP_THREAD_LOCAL float wn[128],fn[128], we[256],fe[256];
0147
0148   static CLHEP_THREAD_LOCAL bool ziggurat_is_init;
0149
0150   static inline unsigned long ziggurat_SHR3(HepRandomEngine* anEngine) {return (unsigned int)(*anEngine);};
0151   static inline float ziggurat_UNI(HepRandomEngine* anEngine) {return anEngine->flat();};
0152   static inline float ziggurat_REXP(HepRandomEngine* anEngine) {
0153     if(!ziggurat_is_init) ziggurat_init();
0154     unsigned long jz=ziggurat_SHR3(anEngine);
0155     unsigned long iz=jz&255;
0156     return (jz<ke[iz]) ? jz*we[iz] : ziggurat_efix(jz,anEngine);
0157   };
0158   static float ziggurat_efix(unsigned long jz,HepRandomEngine* anEngine);
0159
0160 private:
0161
0162   // Private copy constructor. Defining it here disallows use.
0163   RandExpZiggurat(const RandExpZiggurat& d);
0164
0165   std::shared_ptr<HepRandomEngine> localEngine;
0166   double defaultMean;
0167 };
0168
0169 }  // namespace CLHEP
0170
0171 #ifdef ENABLE_BACKWARDS_COMPATIBILITY
0172 //  backwards compatibility will be enabled ONLY in CLHEP 1.9
0173 using namespace CLHEP;
0174 #endif
0175
0176 namespace CLHEP {
0177
0178 inline RandExpZiggurat::RandExpZiggurat(HepRandomEngine & anEngine, double mean ) : localEngine(&anEngine, do_nothing_deleter()), defaultMean(mean)
0179 {
0180 }
0181
0182 inline RandExpZiggurat::RandExpZiggurat(HepRandomEngine * anEngine, double mean ) : localEngine(anEngine), defaultMean(mean)
0183 {
0184 }
0185
0186 }  // namespace CLHEP
0187
0188 #endif