ADMB Documentation  11.1.2490
 All Classes Files Functions Variables Typedefs Friends Defines
dfgammp.cpp
Go to the documentation of this file.
00001 
00017 #include <fvar.hpp>
00018 #define ITMAX 100
00019 #define EPS 1.0e-9
00020 //#define EPS 3.0e-7
00021 #define FPMIN 1.0e-30
00022 double get_values(double x,double y,int print_switch);
00023 
00030 void gcf(const dvariable& _gammcf,const dvariable& a,
00031   const dvariable& x,const dvariable& _gln)
00032 {
00033   ADUNCONST(dvariable,gln)
00034   ADUNCONST(dvariable,gammcf)
00035   int i;
00036   dvariable an,b,c,d,del,h;
00037 
00038   gln=gammln(a);
00039   b=x+1.0-a;
00040   c=1.0/FPMIN;
00041   d=1.0/b;
00042   h=d;
00043   for (i=1;i<=ITMAX;i++) {
00044     an = -i*(i-a);
00045     b += 2.0;
00046     d=an*d+b;
00047     if (fabs(value(d)) < FPMIN) d=FPMIN;
00048     c=b+an/c;
00049     if (fabs(value(c)) < FPMIN) c=FPMIN;
00050     d=1.0/d;
00051     del=d*c;
00052     h *= del;
00053     if (fabs(value(del)-1.0) < EPS) break;
00054   }
00055   if (i > ITMAX)
00056     cerr << "a too large, ITMAX too small in gcf" << endl;
00057   gammcf=exp(-x+a*log(x)-(gln))*h;
00058 }
00059 
00066 void gser(const dvariable& _gamser,const dvariable& a,
00067   const dvariable& x,const dvariable& _gln)
00068 {
00069   ADUNCONST(dvariable,gln)
00070   ADUNCONST(dvariable,gamser)
00071   int n;
00072   dvariable sum,del,ap;
00073 
00074   gln=gammln(a);
00075 
00076   if (value(x) <= 0.0) {
00077     if (value(x) < 0.0)
00078       cerr << "x less than 0 in routine gser" << endl;
00079     gamser=0.0;
00080     return;
00081   }
00082   else
00083   {
00084     ap=a;
00085     del=sum=1.0/a;
00086     for (n=1;n<=ITMAX;n++) {
00087       ap+=1.0;
00088       del *= x/ap;
00089       sum += del;
00090       if (fabs(value(del)) < fabs(value(sum))*EPS) {
00091         gamser=sum*exp(-x+a*log(x)-(gln));
00092         return;
00093       }
00094     }
00095     cerr << "a too large, ITMAX too small in routine gser" << endl;
00096     return;
00097   }
00098 }
00099 
00100 dvariable cumd_gamma(const dvariable& x, const dvariable& a)
00101 {
00102   dvariable gamser,gammcf,gln;
00103 
00104   if (value(x) < 0.0 || value(a) <= 0.0)
00105     cerr << "Invalid arguments in routine gammp" << endl;
00106   if (value(x) < (value(a)+1.0)) {
00107     gser(gamser,a,x,gln);
00108     return gamser;
00109   } else {
00110     gcf(gammcf,a,x,gln);
00111     return 1.0-gammcf;
00112   }
00113 }