sysrap/tests/S4MTRandGaussQ.cc

0001
0002 #include <cmath>
0003 #include "S4MTRandGaussQ.h"
0004 #include "Randomize.hh"
0005
0006 double S4MTRandGaussQ::shoot()
0007 {
0008     double u = G4UniformRand() ;
0009     return transformQuick(u);
0010 }
0011
0012 //
0013 // Table of errInts, for use with transform(r) and quickTransform(r)
0014 //
0015
0016 // Since all these are this is static to this compilation unit only, the
0017 // info is establised a priori and not at each invocation.
0018
0019 // The main data is of course the gaussQTables table; the rest is all
0020 // bookkeeping to know what the tables mean.
0021
0022 #define Table0size   250
0023 #define Table1size  1000
0024 #define TableSize   (Table0size+Table1size)
0025
0026 #define Table0step  (2.0E-6)
0027 #define Table1step  (5.0E-4)
0028
0029 #define Table0scale   (1.0/Table1step)
0030
0031 #define Table0offset 0
0032 #define Table1offset (Table0size)
0033
0034 static const float gaussTables [TableSize] = {
0035 #include "gaussQTables.cdat"
0036 };
0037
0038
0039
0040 double S4MTRandGaussQ::transformQuick (double r)
0041 {
0042   double sign = +1.0; // We always compute a negative number of
0043                         // sigmas.  For r > 0 we will multiply by
0044                         // sign = -1 to return a positive number.
0045   if ( r > .5 ) {
0046     r = 1-r;
0047     sign = -1.0;
0048   }
0049
0050   int index;
0051   double  dx;
0052
0053   if ( r >= Table1step ) {
0054     index = int((Table1size<<1) * r); // 1 to Table1size
0055     if (index == Table1size) return 0.0;
0056     dx = (Table1size<<1) * r - index;  // fraction of way to next bin
0057     index += Table1offset-1;
0058   } else if ( r > Table0step ) {
0059     double rr = r * Table0scale;
0060     index = int(Table0size * rr);    // 1 to Table0size
0061     dx = Table0size * rr - index;      // fraction of way to next bin
0062     index += Table0offset-1;
0063   } else {                             // r <= Table0step - not in tables
0064     return sign*transformSmall(r);
0065   }
0066
0067   double y0 = gaussTables [index++];
0068   double y1 = gaussTables [index];
0069
0070   return (float) (sign * ( y1 * dx + y0 * (1.0-dx) ));
0071
0072 } // transformQuick()
0073
0074
0075
0076 double S4MTRandGaussQ::transformSmall (double r)
0077 {
0078   // Solve for -v in the asymtotic formula
0079   //
0080   // errInt (-v) =  exp(-v*v/2)         1     1*3    1*3*5
0081   //                ------------ * (1 - ---- + ---- - ----- + ... )
0082   //                v*sqrt(2*pi)        v**2   v**4   v**6
0083
0084   // The value of r (=errInt(-v)) supplied is going to less than 2.0E-13,
0085   // which is such that v < -7.25.  Since the value of r is meaningful only
0086   // to an absolute error of 1E-16 (double precision accuracy for a number
0087   // which on the high side could be of the form 1-epsilon), computing
0088   // v to more than 3-4 digits of accuracy is suspect; however, to ensure
0089   // smoothness with the table generator (which uses quite a few terms) we
0090   // also use terms up to 1*3*5* ... *13/v**14, and insist on accuracy of
0091   // solution at the level of 1.0e-7.
0092
0093   // This routine is called less than one time in a million firings, so
0094   // speed is of no concern.  As a matter of technique, we terminate the
0095   // iterations in case they would be infinite, but this should not happen.
0096
0097   double eps = 1.0e-7;
0098   double guess = 7.5;
0099   double v;
0100
0101   for ( int i = 1; i < 50; i++ ) {
0102     double vn2 = 1.0/(guess*guess);
0103     double s = -13*11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2*vn2;
0104              s +=    11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2;
0105              s +=      -9*7*5*3 * vn2*vn2*vn2*vn2*vn2;
0106              s +=         7*5*3 * vn2*vn2*vn2*vn2;
0107              s +=          -5*3 * vn2*vn2*vn2;
0108              s +=             3 * vn2*vn2    - vn2  +    1.0;
0109
0110     v = std::sqrt ( 2.0 * std::log ( s / (r*guess*std::sqrt(2.*M_PI)) ) );
0111     if ( std::fabs(v-guess) < eps ) break;
0112     guess = v;
0113   }
0114   return -v;
0115
0116 } // transformSmall()
0117
0118