ADMB Documentation  11.1.2535
 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   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   int i;
00040   for (i=1;i<=n+1;i++)
00041   {
00042     double zinv=1.0/(z+i);
00043     x+=c[i]*zinv;
00044     //xdot-=c[i]/square(z+i)*zdot;  since zdot=1.0
00045     xdot-=c[i]*square(zinv);
00046     x2dot+=2.0*c[i]*cube(zinv);
00047     x3dot-=6.0*c[i]*fourth(zinv);
00048   }
00049   double t=z+n+0.5;
00050   //return lpp + (z+0.5)*log(t) -t + log(x);
00051   double ans= lpp + (z+0.5)*log(t) -t + log(x);
00052   //double tdot=zdot;
00053   //double dfx=zdot*log(t) + (z+0.5)/t*tdot -tdot +xdot/x;
00054   // since tdot=1.0
00055   // since zdot=1.0
00056   double dfx=log(t) + (z+0.5)/t -1.0 +xdot/x;
00057   double d2f=2.0/t -(z+0.5)/square(t)+x2dot/x
00058     -square(xdot)/square(x);
00059   double d3f=-3.0/square(t) + 2.0*(z+0.5)/cube(t)+ x3dot/x
00060     -3.0*x2dot*xdot/square(x)+2.0*cube(xdot)/cube(x);
00061 
00062   double * xd=zz.get_u_dot();
00063   double * zd=u.get_u_dot();
00064   *u.get_u()=ans;
00065   for (i=0;i<df1b2variable::nvar;i++)
00066   {
00067     *zd++ =dfx * *xd++;
00068   }
00069 
00070   if (!df1b2_gradlist::no_derivatives)
00071     f1b2gradlist->write_pass1(&zz,&u,dfx,d2f,d3f);
00072   return(u);
00073 }
00074 
00079 df1b2variable gammln(const df1b2variable& z)
00080 {
00081   const double lpi =1.1447298858494001741434272;
00082   const double pi =3.1415926535897932384626432;
00083   if (value(z)<0.5)
00084   {
00085     return lpi - log(sin(pi*z)) - gammlnguts(1.0-z);
00086   }
00087   else
00088   {
00089     return gammlnguts(z);
00090   }
00091 }
00092 
00097 df1b2vector gammln(const df1b2vector&  z){
00098   int from=z.indexmin();
00099   int to=z.indexmax();
00100   df1b2vector ret(from,to);
00101   for(int i=from; i<=to; ++i){
00102     ret(i)=gammln(z(i));
00103  }
00104  return(ret);
00105 }
00106 
00111 df1b2variable log_comb(const df1b2variable& n, double k)
00112 {
00113   return factln(n)-factln(k)-factln(n-k);
00114 }
00115 
00120 df1b2variable log_comb(const df1b2variable& n, const df1b2variable& k)
00121 {
00122   return factln(n)-factln(k)-factln(n-k);
00123 }
00124 
00129 df1b2variable log_comb(double n, const df1b2variable& k)
00130 {
00131   return factln(n)-factln(k)-factln(n-k);
00132 }
00133 
00138 df1b2variable factln(const df1b2variable& n)
00139 {
00140   return gammln(n+1.0);
00141 }