Gammaq.cxx

Go to the documentation of this file.
00001 
00010 #include "Gammaq.h"
00011 
00012 #include <iostream>
00013 
00014 #include <cassert>
00015 #include <cmath>
00016 
00017 using std::cerr;
00018 using std::endl;
00019 using std::abs;
00020 
00021 namespace hippodraw {
00022   namespace Numeric {
00023 
00024 double gammln( double xx )
00025 {
00026   double x, y, tmp, ser;
00027   static double cof[6]={76.18009172947146,-86.50532032941677,
00028                         24.01409824083091,-1.231739572450155,
00029                         0.1208650973866179e-2,-0.5395239384953e-5};
00030   
00031   y = x = xx;
00032   tmp = x + 5.5;
00033   tmp -= ( x + 0.5 ) * log( tmp );
00034   ser = 1.000000000190015;
00035   
00036   for ( int j = 0; j <= 5; j++)
00037     ser += cof[j] / (++y);
00038   
00039   return -tmp + log( 2.5066282746310005 * ser / x);
00040 }
00041 
00042 
00043 void gser(double *gamser, double a, double x, double *gln)
00044 {
00045   int n;
00046   double sum, del, ap;
00047   
00048   *gln = gammln( a );
00049 
00050   if ( x <= 0.0 )
00051     {
00052       if (x < 0.0) cerr << "Error: x less than 0 in routine gser" << endl;
00053       *gamser=0.0;
00054       return;
00055     }
00056   else
00057     {
00058       ap  = a;
00059       del = sum =1.0 / a;
00060       for (n = 1; n <= ITMAX; n++)
00061         {
00062           ++ap;
00063           del *= x/ap;
00064           sum += del;
00065           
00066           if (abs(del) < abs(sum)*EPS)
00067             {
00068               *gamser = sum * exp( -x + a * log( x ) -( *gln ) );
00069               return;
00070             }
00071         }
00072       
00073       cerr << "ERROR: a too large, ITMAX too small in routine gser" << endl;
00074       return;
00075     }
00076 }
00077 
00078 void gcf( double *gammcf, double a, double x, double *gln )
00079 {
00080   int i;
00081   double an,b,c,d,del,h;
00082   
00083   *gln = gammln( a );
00084   b = x + 1.0 - a;
00085   c = 1.0 / FPMIN;
00086   d = 1.0 / b;
00087   h = d;
00088   
00089   for( i = 1; i <= ITMAX; i++ )
00090     {
00091       an = -i*(i-a);
00092       b += 2.0;
00093       d=an*d+b;
00094       
00095       if ( abs( d ) < FPMIN )
00096         d = FPMIN;
00097       
00098       c = b + an / c;
00099       
00100       if ( abs( c ) < FPMIN )
00101         c = FPMIN;
00102       
00103       d   = 1.0 / d;
00104       del = d * c;
00105       h  *= del;
00106       
00107       if ( abs( del - 1.0 ) < EPS)
00108         break;
00109     }
00110   
00111   if (i > ITMAX) cerr << "a too large, ITMAX too small in gcf" << endl;
00112   
00113   *gammcf = exp( -x + a * log( x ) -( *gln ) ) * h;
00114 }
00115 
00116 double gammq(double a, double x)
00117 {
00118   void gcf(double *gammcf, double a, double x, double *gln);
00119   
00120   double gamser,gammcf,gln;
00121   
00122   if (x < 0.0 || a <= 0.0)
00123     {
00124       cerr << "a = " << a << " x = " << x << endl;
00125       cerr << "Invalid arguments in routine GAMMQ" << endl;
00126       assert(0);
00127     }
00128 
00129   if (x < (a+1.0))
00130     {
00131       gser(&gamser,a,x,&gln);
00132       return 1.0-gamser;
00133     }
00134   else
00135     {
00136       gcf(&gammcf,a,x,&gln);
00137       return gammcf;
00138     }
00139 }
00140 
00141 } // end namespace Numeric
00142 } // end namespace hippodraw
00143 
00150 /*
00151 int main(void)
00152 {
00153   char txt[MAXSTR];
00154   int i,nval;
00155   float a,val,x;
00156   FILE *fp;
00157   
00158   if ((fp = fopen("fncval.dat","r")) == NULL)
00159     printf("Data file fncval.dat not found\n");
00160 
00161   std::fgets(txt,MAXSTR,fp);
00162   
00163   while ( std::strncmp(txt,"Incomplete Gamma Function",25))
00164     {
00165       std::fgets(txt,MAXSTR,fp);
00166       if ( std::feof(fp))
00167         std::printf("Data not found in fncval.dat\n");
00168     }
00169   
00170   std::fscanf(fp,"%d %*s",&nval);
00171   std::printf("\n%s\n",txt);
00172   std::printf("%4s %11s %14s %14s \n","a","x","actual","gammq(a,x)");
00173 
00174   for (i=1;i<=nval;i++)
00175     {
00176       std::fscanf(fp,"%f %f %f",&a,&x,&val);
00177       std::printf("%6.2f %12.6f %12.6f %12.6f\n", a, x, (1.0-val), gammq(a,x));
00178     }
00179   
00180   std::fclose(fp);
00181   return 0;
00182 }
00183 
00184 */
00185 
00186 #undef ITMAX
00187 #undef EPS
00188 #undef FPMIN
00189 #undef MAXSTR

Generated for HippoDraw Class Library by doxygen