QMcBeaver: MathFunctions.cpp Source File

00001 //            QMcBeaver
00002 //
00003 //         Constructed by 
00004 //
00005 //     Michael Todd Feldmann 
00006 //              and 
00007 //   David Randall "Chip" Kent IV
00008 //
00009 // Copyright 2000-2.  All rights reserved.
00010 //
00011 // drkent@users.sourceforge.net mtfeldmann@users.sourceforge.net
00012 
00013 #include "MathFunctions.h"
00014 
00015 double MathFunctions::erf(double x)
00016 {
00017   // Series on [0,1]
00018   int ERFC_COEF_LENGTH = 14;
00019   double ERFC_COEF[] = 
00020   {
00021     -.490461212346918080399845440334e-1,
00022     -.142261205103713642378247418996e0,
00023     .100355821875997955757546767129e-1,
00024     -.576876469976748476508270255092e-3,
00025     .274199312521960610344221607915e-4,
00026     -.110431755073445076041353812959e-5,
00027     .384887554203450369499613114982e-7,
00028     -.118085825338754669696317518016e-8,
00029     .323342158260509096464029309534e-10,
00030     -.799101594700454875816073747086e-12,
00031     .179907251139614556119672454866e-13,
00032     -.371863548781869263823168282095e-15,
00033     .710359900371425297116899083947e-17,
00034     -.126124551191552258324954248533e-18
00035   };
00036 
00037   double  ans;
00038   double  y = fabs(x);
00039   
00040   if (y <= 1.49012e-08) 
00041     {
00042       // 1.49012e-08 = Math.sqrt(2*EPSILON_SMALL)
00043       ans = 2 * x / 1.77245385090551602729816748334;
00044     } 
00045   else if (y <= 1) 
00046     {
00047       ans = x * (1 + csevl(2 * x * x - 1, ERFC_COEF, ERFC_COEF_LENGTH));
00048     } 
00049   else if (y < 6.013687357) 
00050     {
00051       // 6.013687357 = 
00052       // Math.sqrt(-Math.log(1.77245385090551602729816748334 * EPSILON_SMALL))
00053       ans = sign(1 - erfc(y), x);
00054     } 
00055   else 
00056     {
00057       ans = sign(1, x);
00058     }
00059 
00060   return ans;
00061 }
00062 
00063 double MathFunctions::erfc(double x)
00064 {
00065   // Series on [0,1]
00066   int ERFC_COEF_LENGTH = 14;
00067   static double ERFC_COEF[] = 
00068   {
00069     -.490461212346918080399845440334e-1,
00070     -.142261205103713642378247418996e0,
00071     .100355821875997955757546767129e-1,
00072     -.576876469976748476508270255092e-3,
00073     .274199312521960610344221607915e-4,
00074     -.110431755073445076041353812959e-5,
00075     .384887554203450369499613114982e-7,
00076     -.118085825338754669696317518016e-8,
00077     .323342158260509096464029309534e-10,
00078     -.799101594700454875816073747086e-12,
00079     .179907251139614556119672454866e-13,
00080     -.371863548781869263823168282095e-15,
00081     .710359900371425297116899083947e-17,
00082     -.126124551191552258324954248533e-18
00083   };
00084 
00085   // Series on [0.25,1.00]
00086   int ERFC2_COEF_LENGTH = 27;
00087   static double ERFC2_COEF[] = 
00088   {
00089     -.69601346602309501127391508262e-1,
00090     -.411013393626208934898221208467e-1,
00091     .391449586668962688156114370524e-2,
00092     -.490639565054897916128093545077e-3,
00093     .715747900137703638076089414183e-4,
00094     -.115307163413123283380823284791e-4,
00095     .199467059020199763505231486771e-5,
00096     -.364266647159922287393611843071e-6,
00097     .694437261000501258993127721463e-7,
00098     -.137122090210436601953460514121e-7,
00099     .278838966100713713196386034809e-8,
00100     -.581416472433116155186479105032e-9,
00101     .123892049175275318118016881795e-9,
00102     -.269063914530674343239042493789e-10,
00103     .594261435084791098244470968384e-11,
00104     -.133238673575811957928775442057e-11,
00105     .30280468061771320171736972433e-12,
00106     -.696664881494103258879586758895e-13,
00107     .162085454105392296981289322763e-13,
00108     -.380993446525049199987691305773e-14,
00109     .904048781597883114936897101298e-15,
00110     -.2164006195089607347809812047e-15,
00111     .522210223399585498460798024417e-16,
00112     -.126972960236455533637241552778e-16,
00113     .310914550427619758383622741295e-17,
00114     -.766376292032038552400956671481e-18,
00115     .190081925136274520253692973329e-18
00116   };
00117 
00118   // Series on [0,0.25]
00119   int ERFCC_COEF_LENGTH = 29;
00120   static double ERFCC_COEF[] = 
00121   {
00122     .715179310202924774503697709496e-1,
00123     -.265324343376067157558893386681e-1,
00124     .171115397792085588332699194606e-2,
00125     -.163751663458517884163746404749e-3,
00126     .198712935005520364995974806758e-4,
00127     -.284371241276655508750175183152e-5,
00128     .460616130896313036969379968464e-6,
00129     -.822775302587920842057766536366e-7,
00130     .159214187277090112989358340826e-7,
00131     -.329507136225284321486631665072e-8,
00132     .72234397604005554658126115389e-9,
00133     -.166485581339872959344695966886e-9,
00134     .401039258823766482077671768814e-10,
00135     -.100481621442573113272170176283e-10,
00136     .260827591330033380859341009439e-11,
00137     -.699111056040402486557697812476e-12,
00138     .192949233326170708624205749803e-12,
00139     -.547013118875433106490125085271e-13,
00140     .158966330976269744839084032762e-13,
00141     -.47268939801975548392036958429e-14,
00142     .14358733767849847867287399784e-14,
00143     -.444951056181735839417250062829e-15,
00144     .140481088476823343737305537466e-15,
00145     -.451381838776421089625963281623e-16,
00146     .147452154104513307787018713262e-16,
00147     -.489262140694577615436841552532e-17,
00148     .164761214141064673895301522827e-17,
00149     -.562681717632940809299928521323e-18,
00150     .194744338223207851429197867821e-18
00151   };
00152 
00153   double ans;
00154   double y = fabs(x);
00155   
00156   if (x <= -6.013687357) 
00157     {
00158       // -6.013687357 = 
00159       // -Math.sqrt(-Math.log(1.77245385090551602729816748334 * EPSILON_SMALL))
00160       ans = 2;
00161     } 
00162   else if (y < 1.49012e-08) 
00163     {
00164       // 1.49012e-08 = Math.sqrt(2*EPSILON_SMALL)
00165       ans = 1 - 2*x/1.77245385090551602729816748334;
00166     } 
00167   else 
00168     {
00169       double ysq = y*y;
00170     
00171       if (y < 1) 
00172         {
00173           ans = 1 - x*(1+csevl(2*ysq-1,ERFC_COEF,ERFC_COEF_LENGTH));
00174         } 
00175       else if (y <= 4.0) 
00176         {
00177           ans = exp(-ysq)/y*(0.5+csevl((8.0/ysq-5.0)/3.0,ERFC2_COEF,
00178                                        ERFC2_COEF_LENGTH));
00179           if (x < 0)  ans = 2.0 - ans;
00180           if (x < 0)  ans = 2.0 - ans;
00181           if (x < 0)  ans = 2.0 - ans;
00182         } 
00183       else 
00184         {
00185           ans = exp(-ysq)/y*(0.5+csevl(8.0/ysq-1,ERFCC_COEF,
00186                                        ERFCC_COEF_LENGTH));
00187           if (x < 0)  ans = 2.0 - ans;
00188         }
00189     }
00190 
00191   return ans;
00192 }
00193 
00194 void MathFunctions::fitU(double Z, double Frk, double r, double delta,
00195                          double & zeta, double & a, double & I)
00196 {
00197   //Now we solve for a and z (looking at eqns 41 & 42)
00198   double t1 = Frk + Z;
00199   double discriminant = t1*(r*t1 - 4.0)/r;
00200   zeta = 0.0;
00201   if(discriminant > 0.0)
00202     {
00203       double root1 = Z - Frk;
00204       double root2 = Z - Frk;
00205       double sq = sqrt(discriminant);     
00206       root1 -= sq;
00207       root2 += sq;
00208       root1 /= 2.0;
00209       root2 /= 2.0;
00210       
00211       if(root1*root2 <= 0.0)
00212         {
00213           //return the postive root
00214           zeta = root1 > root2 ? root1 : root2;
00215         } else if(root1 <= 0.0)
00216           {
00217             //both were negative
00218             zeta = 0.0;      
00219           } else {
00220             //both positive, return the smallest
00221             zeta = root1 < root2 ? root1 : root2;
00222           }
00223     }
00224   
00225   if(zeta <= 1e-50)
00226     {
00227       if(Frk < 0.0)
00228         zeta = -Frk;
00229       else
00230         zeta = 1.0;
00231     }
00232   
00233   a = -(Frk + zeta)/(-1.0 + r*(Frk + zeta));  
00234 
00235   /*
00236     This is from eqns on pg 158 (In the appendix).
00237     Note: there is a typo in the formulas. Umrigar says the nodes
00238     are at -a, but as is easily verified, the nodes are actually at -1/a.
00239     
00240     Because we're integrating the absolute value of the function, we need
00241     to separate this into two integrals if the integrand goes through a node.
00242   */
00243   double node = -1.0/a;
00244   if(fabs(a) > 1e-50 && r/delta < node && node < r*delta)
00245     {
00246       I = fabs(accel_G(1.5,zeta*r*delta,-zeta/a,a,zeta)) +
00247           fabs(accel_G(1.5,-zeta/a,zeta*r/delta,a,zeta));
00248     } else {
00249       I = fabs(accel_G(1.5,zeta*r*delta,zeta*r/delta,a,zeta));
00250     }
00251 
00252   /*
00253   printf("a=%20.10f zeta=%20.10f r=%20.10f root=%20.10f L=%20.10f U=%20.10f I=%20.10f\n",
00254          a,zeta,r,node,r/delta,r*delta,I);
00255   */
00256 }
00257 
00258 double MathFunctions::accel_G(double n, double x, double y,
00259                               double a, double zeta)
00260 {
00261   double den = pow(zeta,-n-1);
00262   double t1 = (gamma_inc(n,x)   - gamma_inc(n,y)  )*den*zeta;
00263   double t2 = (gamma_inc(n+1,x) - gamma_inc(n+1,y))*den;
00264   return t1 + a*t2;
00265 }
00266 
00267 double MathFunctions::csevl(double x, double * coef, int lengthCoef)
00268 {
00269   double b0, b1, b2, twox;
00270   int    i;
00271 
00272   b1 = 0.0;
00273   b0 = 0.0;
00274   b2 = 0.0;
00275   twox = 2.0*x;
00276 
00277   for (i = lengthCoef-1;  i >= 0;  i--) 
00278     {
00279       b2 = b1;
00280       b1 = b0;
00281       b0 = twox*b1 - b2 + coef[i];
00282     }
00283 
00284   return 0.5*(b0-b2);
00285 }
00286 
00287 double MathFunctions::sign(double x, double y)
00288 {
00289   double abs_x = ((x < 0) ? -x : x);
00290   return (y < 0.0) ? -abs_x : abs_x;
00291 }
00292 
00293 #define ITMAX 100
00294 #define EPS 3.0e-7
00295 #define FPMIN 1.0e-30
00296 
00297 double MathFunctions::F_gamma(double v, double t)
00298 {
00299   /* Computes the function 
00300       Integrate[u^(2 v) Exp(-t u^2), {u, 0, 1}]
00301       = (1/2) t^(-1/2 - v) * Gamma(1/2 + v, t)
00302   */
00303   if (t < 0.000001) 
00304     return 1.0 / (2.0 * v + 1.0); /* limiting case */
00305   else
00306     return 0.5 * pow(t, -0.5 - v) * gamma_inc(0.5 + v, t);
00307 }
00308 
00309 double MathFunctions::gamma_inc(double a, double x)
00310 {
00311   /* From Numerical Recipes in Fortran, gammp.f 
00312      Returns the incomplete gamma function 
00313         Integrate[Exp[-t] t^(a - 1), {t, 0, x}] / Gamma[a].
00314      Picks the series or continued fraction algorithm depending on input. 
00315   */
00316 
00317   if (x < 0 || a < 0) {printf("Bad argument in gamma function.\n"); exit(1);}
00318   if (x < a + 1)
00319     return gamma_series(a, x);
00320   else
00321     return gamma_cf(a, x);
00322 }
00323 
00324 double MathFunctions::gamma_log(double x2)
00325 {
00326   /* From Numerical Recipes in C, gammln.c
00327      Returns the value Log(gamma(x)), x > 0
00328   */
00329   
00330   double x, y, tmp, ser;
00331   static double cof[] = {76.18009172947146, -86.50532032941677, 
00332                          24.01409824083091, -1.231739572450155,
00333                          0.1208650973866179E-2, -0.5395239384953E-5, 
00334                          2.5066282746310005};
00335   int j;
00336 
00337   y = x = x2;
00338   tmp = x + 5.5;
00339   tmp -= (x + 0.5) * log(tmp);
00340   ser = 1.000000000190015;
00341 
00342   for (j = 0; j < 6; j++)
00343     {
00344       y += 1.0;
00345       ser +=  cof[j] / y;
00346     }
00347   return -tmp + log(2.5066282746310005 * ser/x);
00348 }
00349 
00350 double MathFunctions::gamma_series(double a, double x)
00351 {
00352   /* From Numerical Recipes in C, gser.c
00353      Returns the incomplete gamma function
00354         Integrate[Exp[-t] t^(a - 1), {t, 0, x}].
00355      Uses a series expansion.
00356   */
00357 
00358   if (x <= 0.0) {printf("x less than 0 in routine gamma_series.\n"); exit(1);}
00359  
00360   int n; 
00361   double sum, del, ap;
00362 
00363   ap = a;
00364   del = sum = 1.0 / a;
00365   for (n = 1; n <= ITMAX; n++)
00366   {
00367     ++ap;
00368     del *= x / ap;
00369     sum += del;
00370     if (fabs(del) < fabs(sum) * EPS)
00371       return sum * exp(-x + a * log(x));
00372   }
00373   printf("a too large, ITMAX too small in routine gamma_series.\n"); exit(1);
00374   return 0.0;
00375 }
00376 
00377 double MathFunctions::gamma_cf(double a, double x)
00378 {
00379   /* From Numerical Recipes in C, gcf.c
00380      Returns the incomplete gamma function
00381         Integrate[Exp[-t] t^(a - 1), {t, 0, x}] / Gamma[a].
00382      Uses a continued fraction expansion.
00383   */
00384 
00385   double an, b, c, d, del, h;
00386 
00387   b = x + 1.0 - a;
00388   c = 1.0 / FPMIN;
00389   d = 1.0 / b;
00390   h = d;
00391   int i;
00392   for (i = 1; i <= ITMAX; i++)
00393   {
00394     an = -i * (i - a);
00395     b += 2.0;
00396     d = an * d + b;
00397     if (fabs(d) < FPMIN) d = FPMIN;
00398     c = b + an / c;
00399     if (fabs(c) < FPMIN) c = FPMIN;
00400     d = 1.0 / d;
00401     del = d * c;
00402     h *= del;
00403     if (fabs(del - 1.0) < EPS) break;
00404   }
00405   if (i > ITMAX) 
00406     {
00407       printf("a too large, ITMAX too small in gamma_cf\n"); 
00408       exit(1);
00409     }
00410   return exp(gamma_log(a)) - exp(-x + a * log(x)) * h;
00411 }
00412 
00413 double MathFunctions::legendre(int l, double x)
00414 {
00415   if ( l == 0 )
00416     return 1.0;
00417   if ( l == 1 )
00418     return x;
00419   
00420   return ( (2.0*l-1.0)*x*legendre(l-1, x) - (l-1.0)*legendre(l-2, x) ) / l;
00421 }
00422 
00423 double MathFunctions::rij(Array2D<double> &position, int i, int j)
00424 {
00425   double r;
00426   r = sqrt((position(i,0)-position(j,0))*
00427            (position(i,0)-position(j,0))+
00428            (position(i,1)-position(j,1))*
00429            (position(i,1)-position(j,1))+
00430            (position(i,2)-position(j,2))*
00431            (position(i,2)-position(j,2)));
00432   return r;
00433 }
00434 
00435 double MathFunctions::rij(Array2D<double> &positioni,
00436                           Array2D<double> &positionj,
00437                           int i, int j)
00438 {
00439   double r;
00440   r = sqrt((positioni(i,0)-positionj(j,0))*
00441            (positioni(i,0)-positionj(j,0))+
00442            (positioni(i,1)-positionj(j,1))*
00443            (positioni(i,1)-positionj(j,1))+
00444            (positioni(i,2)-positionj(j,2))*
00445            (positioni(i,2)-positionj(j,2)));
00446   return r;
00447 }
00448 
00449 #define EXP_A (1048576/M_LN2) 
00450 #define EXP_C 60801     
00451 double MathFunctions::fast_exp(double y)
00452 {
00453   union 
00454   { 
00455     double d; 
00456 #ifdef LITTLE_ENDIAN 
00457     struct { int j, i; } n; 
00458 #else 
00459     struct { int i, j; } n; 
00460 #endif 
00461   } 
00462   _eco; 
00463   _eco.n.i = (int)(EXP_A*(y)) + (1072693248 - EXP_C); 
00464   _eco.n.j = 0; 
00465   return  _eco.d; 
00466 }