ADMB Documentation  11.1x.2706
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2f20.cpp
Go to the documentation of this file.
00001 
00011 #include <df1b2fun.h>
00012 
00017 df1b2variable gammlnguts(const df1b2variable _zz)
00018 {
00019   ADUNCONST(df1b2variable,zz)
00020   df1b2variable u;
00021   double  z = value(_zz);
00022   //double zdot=1.0;
00023   const double lpp =0.9189385332046727417803297;
00024   unsigned int n=7;
00025   const double c[9]={0.99999999999980993,
00026     676.5203681218851,
00027     -1259.1392167224028,
00028      771.32342877765313,
00029     -176.61502916214059,
00030     12.507343278686905,
00031      -0.13857109526572012,
00032     9.9843695780195716e-6,
00033     1.5056327351493116e-7};
00034   z-=1.0;
00035   double x=c[0];
00036   double xdot=0.0;
00037   double x2dot=0.0;
00038   double x3dot=0.0;
00039   for (unsigned int i=1;i<=n+1;i++)
00040   {
00041     double zinv=1.0/(z+i);
00042     x+=c[i]*zinv;
00043     //xdot-=c[i]/square(z+i)*zdot;  since zdot=1.0
00044     xdot-=c[i]*square(zinv);
00045     x2dot+=2.0*c[i]*cube(zinv);
00046     x3dot-=6.0*c[i]*fourth(zinv);
00047   }
00048   double t=z+n+0.5;
00049   //return lpp + (z+0.5)*log(t) -t + log(x);
00050   double ans= lpp + (z+0.5)*log(t) -t + log(x);
00051   //double tdot=zdot;
00052   //double dfx=zdot*log(t) + (z+0.5)/t*tdot -tdot +xdot/x;
00053   // since tdot=1.0
00054   // since zdot=1.0
00055   double dfx=log(t) + (z+0.5)/t -1.0 +xdot/x;
00056   double d2f=2.0/t -(z+0.5)/square(t)+x2dot/x
00057     -square(xdot)/square(x);
00058   double d3f=-3.0/square(t) + 2.0*(z+0.5)/cube(t)+ x3dot/x
00059     -3.0*x2dot*xdot/square(x)+2.0*cube(xdot)/cube(x);
00060 
00061   double * xd=zz.get_u_dot();
00062   double * zd=u.get_u_dot();
00063   *u.get_u()=ans;
00064   for (unsigned int i=0;i<df1b2variable::nvar;i++)
00065   {
00066     *zd++ =dfx * *xd++;
00067   }
00068 
00069   if (!df1b2_gradlist::no_derivatives)
00070     f1b2gradlist->write_pass1(&zz,&u,dfx,d2f,d3f);
00071   return(u);
00072 }
00073 
00078 df1b2variable gammln(const df1b2variable& z)
00079 {
00080   const double lpi =1.1447298858494001741434272;
00081   const double pi =3.1415926535897932384626432;
00082   if (value(z)<0.5)
00083   {
00084     return lpi - log(sin(pi*z)) - gammlnguts(1.0-z);
00085   }
00086   else
00087   {
00088     return gammlnguts(z);
00089   }
00090 }
00091 
00096 df1b2vector gammln(const df1b2vector&  z){
00097   int from=z.indexmin();
00098   int to=z.indexmax();
00099   df1b2vector ret(from,to);
00100   for(int i=from; i<=to; ++i){
00101     ret(i)=gammln(z(i));
00102  }
00103  return(ret);
00104 }
00105 
00110 df1b2variable log_comb(const df1b2variable& n, double k)
00111 {
00112   return factln(n)-factln(k)-factln(n-k);
00113 }
00114 
00119 df1b2variable log_comb(const df1b2variable& n, const df1b2variable& k)
00120 {
00121   return factln(n)-factln(k)-factln(n-k);
00122 }
00123 
00128 df1b2variable log_comb(double n, const df1b2variable& k)
00129 {
00130   return factln(n)-factln(k)-factln(n-k);
00131 }
00132 
00137 df1b2variable factln(const df1b2variable& n)
00138 {
00139   return gammln(n+1.0);
00140 }