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