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 5.2.p2)


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