Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/externals/clhep/src/RandExpZiggurat.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /externals/clhep/src/RandExpZiggurat.cc (Version 11.3.0) and /externals/clhep/src/RandExpZiggurat.cc (Version 10.4.p1)


  1 #include "CLHEP/Random/DoubConv.h"                  1 #include "CLHEP/Random/DoubConv.h"
  2                                                     2 
  3 #include "CLHEP/Random/RandExpZiggurat.h"           3 #include "CLHEP/Random/RandExpZiggurat.h"
  4                                                     4 
  5 #include <cmath>  // for std::log()            << 
  6 #include <iostream>                                 5 #include <iostream>
  7 #include <vector>                              <<   6 #include <cmath>  // for std::log()
  8                                                     7 
  9 namespace CLHEP {                                   8 namespace CLHEP {
 10                                                     9 
 11 CLHEP_THREAD_LOCAL unsigned long RandExpZiggur     10 CLHEP_THREAD_LOCAL unsigned long RandExpZiggurat::kn[128], RandExpZiggurat::ke[256];
 12 CLHEP_THREAD_LOCAL float RandExpZiggurat::wn[1     11 CLHEP_THREAD_LOCAL float RandExpZiggurat::wn[128],RandExpZiggurat::fn[128],RandExpZiggurat::we[256],RandExpZiggurat::fe[256];
 13 CLHEP_THREAD_LOCAL bool RandExpZiggurat::ziggu     12 CLHEP_THREAD_LOCAL bool RandExpZiggurat::ziggurat_is_init = false;
 14                                                    13 
 15 std::string RandExpZiggurat::name() const {ret     14 std::string RandExpZiggurat::name() const {return "RandExpZiggurat";}
 16                                                    15 
 17 HepRandomEngine & RandExpZiggurat::engine() {r     16 HepRandomEngine & RandExpZiggurat::engine() {return *localEngine;}
 18                                                    17 
 19 RandExpZiggurat::~RandExpZiggurat() {              18 RandExpZiggurat::~RandExpZiggurat() {
 20 }                                                  19 }
 21                                                    20 
 22 RandExpZiggurat::RandExpZiggurat(const RandExp     21 RandExpZiggurat::RandExpZiggurat(const RandExpZiggurat& right) : HepRandom(right),defaultMean(right.defaultMean)
 23 {                                                  22 {
 24 }                                                  23 }
 25                                                    24 
 26 double RandExpZiggurat::operator()()               25 double RandExpZiggurat::operator()()
 27 {                                                  26 {
 28   return fire( defaultMean );                      27   return fire( defaultMean );
 29 }                                                  28 }
 30                                                    29 
 31 void RandExpZiggurat::shootArray( const int si     30 void RandExpZiggurat::shootArray( const int size, float* vect, float mean )
 32 {                                                  31 {
 33    for (int i=0; i<size; ++i) vect[i] = shoot(     32    for (int i=0; i<size; ++i) vect[i] = shoot(mean);
 34 }                                                  33 }
 35                                                    34 
 36 void RandExpZiggurat::shootArray( const int si     35 void RandExpZiggurat::shootArray( const int size, double* vect, double mean )
 37 {                                                  36 {
 38    for (int i=0; i<size; ++i) vect[i] = shoot(     37    for (int i=0; i<size; ++i) vect[i] = shoot(mean);
 39 }                                                  38 }
 40                                                    39 
 41 void RandExpZiggurat::shootArray(HepRandomEngi     40 void RandExpZiggurat::shootArray(HepRandomEngine* anEngine, const int size, float* vect, float mean )
 42 {                                                  41 {
 43    for (int i=0; i<size; ++i) vect[i] = shoot(     42    for (int i=0; i<size; ++i) vect[i] = shoot(anEngine, mean);
 44 }                                                  43 }
 45                                                    44 
 46 void RandExpZiggurat::shootArray(HepRandomEngi     45 void RandExpZiggurat::shootArray(HepRandomEngine* anEngine, const int size, double* vect, double mean )
 47 {                                                  46 {
 48    for (int i=0; i<size; ++i) vect[i] = shoot(     47    for (int i=0; i<size; ++i) vect[i] = shoot(anEngine, mean);
 49 }                                                  48 }
 50                                                    49 
 51 void RandExpZiggurat::fireArray( const int siz     50 void RandExpZiggurat::fireArray( const int size, float* vect)
 52 {                                                  51 {
 53    for (int i=0; i<size; ++i) vect[i] = fire(      52    for (int i=0; i<size; ++i) vect[i] = fire( defaultMean );
 54 }                                                  53 }
 55                                                    54 
 56 void RandExpZiggurat::fireArray( const int siz     55 void RandExpZiggurat::fireArray( const int size, double* vect)
 57 {                                                  56 {
 58    for (int i=0; i<size; ++i) vect[i] = fire(      57    for (int i=0; i<size; ++i) vect[i] = fire( defaultMean );
 59 }                                                  58 }
 60                                                    59 
 61 void RandExpZiggurat::fireArray( const int siz     60 void RandExpZiggurat::fireArray( const int size, float* vect, float mean )
 62 {                                                  61 {
 63    for (int i=0; i<size; ++i) vect[i] = fire(      62    for (int i=0; i<size; ++i) vect[i] = fire( mean );
 64 }                                                  63 }
 65                                                    64 
 66 void RandExpZiggurat::fireArray( const int siz     65 void RandExpZiggurat::fireArray( const int size, double* vect, double mean )
 67 {                                                  66 {
 68    for (int i=0; i<size; ++i) vect[i] = fire(      67    for (int i=0; i<size; ++i) vect[i] = fire( mean );
 69 }                                                  68 }
 70                                                    69 
 71 std::ostream & RandExpZiggurat::put ( std::ost     70 std::ostream & RandExpZiggurat::put ( std::ostream & os ) const {
 72   int pr=os.precision(20);                         71   int pr=os.precision(20);
 73   std::vector<unsigned long> t(2);                 72   std::vector<unsigned long> t(2);
 74   os << " " << name() << "\n";                     73   os << " " << name() << "\n";
 75   os << "Uvec" << "\n";                            74   os << "Uvec" << "\n";
 76   t = DoubConv::dto2longs(defaultMean);            75   t = DoubConv::dto2longs(defaultMean);
 77   os << defaultMean << " " << t[0] << " " << t     76   os << defaultMean << " " << t[0] << " " << t[1] << "\n";
 78   os.precision(pr);                                77   os.precision(pr);
 79   return os;                                       78   return os;
 80 #ifdef REMOVED                                     79 #ifdef REMOVED
 81   int pr=os.precision(20);                         80   int pr=os.precision(20);
 82   os << " " << name() << "\n";                     81   os << " " << name() << "\n";
 83   os << defaultMean << "\n";                       82   os << defaultMean << "\n";
 84   os.precision(pr);                                83   os.precision(pr);
 85   return os;                                       84   return os;
 86 #endif                                             85 #endif
 87 }                                                  86 }
 88                                                    87 
 89 std::istream & RandExpZiggurat::get ( std::ist     88 std::istream & RandExpZiggurat::get ( std::istream & is ) {
 90   std::string inName;                              89   std::string inName;
 91   is >> inName;                                    90   is >> inName;
 92   if (inName != name()) {                          91   if (inName != name()) {
 93     is.clear(std::ios::badbit | is.rdstate());     92     is.clear(std::ios::badbit | is.rdstate());
 94     std::cerr << "Mismatch when expecting to r     93     std::cerr << "Mismatch when expecting to read state of a "
 95             << name() << " distribution\n"         94             << name() << " distribution\n"
 96         << "Name found was " << inName             95         << "Name found was " << inName
 97         << "\nistream is left in the badbit st     96         << "\nistream is left in the badbit state\n";
 98     return is;                                     97     return is;
 99   }                                                98   }
100   if (possibleKeywordInput(is, "Uvec", default     99   if (possibleKeywordInput(is, "Uvec", defaultMean)) {
101     std::vector<unsigned long> t(2);              100     std::vector<unsigned long> t(2);
102     is >> defaultMean >> t[0] >> t[1]; default    101     is >> defaultMean >> t[0] >> t[1]; defaultMean = DoubConv::longs2double(t); 
103     return is;                                    102     return is;
104   }                                               103   }
105   // is >> defaultMean encompassed by possible    104   // is >> defaultMean encompassed by possibleKeywordInput
106   return is;                                      105   return is;
107 }                                                 106 }
108                                                   107 
109                                                   108 
110 float RandExpZiggurat::ziggurat_efix(unsigned     109 float RandExpZiggurat::ziggurat_efix(unsigned long jz,HepRandomEngine* anEngine)
111 {                                                 110 { 
112   if(!ziggurat_is_init) ziggurat_init();          111   if(!ziggurat_is_init) ziggurat_init();
113                                                   112 
114   unsigned long iz=jz&255;                        113   unsigned long iz=jz&255;
115                                                   114   
116   float x;                                        115   float x;
117   for(;;)                                         116   for(;;)
118   {                                               117   {  
119     if(iz==0) return (7.69711-std::log(ziggura    118     if(iz==0) return (7.69711-std::log(ziggurat_UNI(anEngine)));          /* iz==0 */
120     x=jz*we[iz];                                  119     x=jz*we[iz]; 
121     if( fe[iz]+ziggurat_UNI(anEngine)*(fe[iz-1    120     if( fe[iz]+ziggurat_UNI(anEngine)*(fe[iz-1]-fe[iz]) < std::exp(-x) ) return (x);
122                                                   121 
123     /* initiate, try to exit for(;;) loop */      122     /* initiate, try to exit for(;;) loop */
124     jz=ziggurat_SHR3(anEngine);                   123     jz=ziggurat_SHR3(anEngine);
125     iz=(jz&255);                                  124     iz=(jz&255);
126     if(jz<ke[iz]) return (jz*we[iz]);             125     if(jz<ke[iz]) return (jz*we[iz]);
127   }                                               126   }
128 }                                                 127 }
129                                                   128 
130 bool RandExpZiggurat::ziggurat_init()             129 bool RandExpZiggurat::ziggurat_init()
131 {                                                 130 {  
132   const double rzm1 = 2147483648.0, rzm2 = 429    131   const double rzm1 = 2147483648.0, rzm2 = 4294967296.;
133   double dn=3.442619855899,tn=dn,vn=9.91256303    132   double dn=3.442619855899,tn=dn,vn=9.91256303526217e-3, q;
134   double de=7.697117470131487, te=de, ve=3.949    133   double de=7.697117470131487, te=de, ve=3.949659822581572e-3;
135   int i;                                          134   int i;
136                                                   135 
137 /* Set up tables for RNOR */                      136 /* Set up tables for RNOR */
138   q=vn/std::exp(-.5*dn*dn);                       137   q=vn/std::exp(-.5*dn*dn);
139   kn[0]=(unsigned long)((dn/q)*rzm1);             138   kn[0]=(unsigned long)((dn/q)*rzm1);
140   kn[1]=0;                                        139   kn[1]=0;
141                                                   140 
142   wn[0]=q/rzm1;                                   141   wn[0]=q/rzm1;
143   wn[127]=dn/rzm1;                                142   wn[127]=dn/rzm1;
144                                                   143 
145   fn[0]=1.;                                       144   fn[0]=1.;
146   fn[127]=std::exp(-.5*dn*dn);                    145   fn[127]=std::exp(-.5*dn*dn);
147                                                   146 
148   for(i=126;i>=1;i--) {                           147   for(i=126;i>=1;i--) {
149     dn=std::sqrt(-2.*std::log(vn/dn+std::exp(-    148     dn=std::sqrt(-2.*std::log(vn/dn+std::exp(-.5*dn*dn)));
150     kn[i+1]=(unsigned long)((dn/tn)*rzm1);        149     kn[i+1]=(unsigned long)((dn/tn)*rzm1);
151     tn=dn;                                        150     tn=dn;
152     fn[i]=std::exp(-.5*dn*dn);                    151     fn[i]=std::exp(-.5*dn*dn);
153     wn[i]=dn/rzm1;                                152     wn[i]=dn/rzm1;
154   }                                               153   }
155                                                   154 
156 /* Set up tables for REXP */                      155 /* Set up tables for REXP */
157   q = ve/std::exp(-de);                           156   q = ve/std::exp(-de);
158   ke[0]=(unsigned long)((de/q)*rzm2);             157   ke[0]=(unsigned long)((de/q)*rzm2);
159   ke[1]=0;                                        158   ke[1]=0;
160                                                   159 
161   we[0]=q/rzm2;                                   160   we[0]=q/rzm2;
162   we[255]=de/rzm2;                                161   we[255]=de/rzm2;
163                                                   162 
164   fe[0]=1.;                                       163   fe[0]=1.;
165   fe[255]=std::exp(-de);                          164   fe[255]=std::exp(-de);
166                                                   165 
167   for(i=254;i>=1;i--) {                           166   for(i=254;i>=1;i--) {
168     de=-std::log(ve/de+std::exp(-de));            167     de=-std::log(ve/de+std::exp(-de));
169     ke[i+1]= (unsigned long)((de/te)*rzm2);       168     ke[i+1]= (unsigned long)((de/te)*rzm2);
170     te=de;                                        169     te=de;
171     fe[i]=std::exp(-de);                          170     fe[i]=std::exp(-de);
172     we[i]=de/rzm2;                                171     we[i]=de/rzm2;
173   }                                               172   }
174   ziggurat_is_init=true;                          173   ziggurat_is_init=true;
175   return true;                                    174   return true;
176 }                                                 175 }
177                                                   176 
178 }  // namespace CLHEP                             177 }  // namespace CLHEP
179                                                   178