ADMB Documentation  11.1.2366
 All Classes Files Functions Variables Typedefs Friends Defines
combc.cpp
Go to the documentation of this file.
00001 /*
00002   $Id: combc.cpp 1112 2013-07-12 21:41:41Z johnoel $
00003 
00004   Author: David Fournier
00005   Copyright (c) 2009-2012 ADMB Foundation
00006  */
00014 #include <fvar.hpp>
00015 
00016 double factln(double n);
00017 double gammln(double xx);
00018 
00027 double log_comb(double n, double k)
00028 {
00029   return factln(n)-factln(k)-factln(n-k);
00030 }
00031 
00038 double factln(double n)
00039 {
00040   return gammln(n+1.0);
00041 }
00042 
00052 double gammln(double xx)
00053 {
00054   double x,tmp,ser;
00055   static double cof[6]={76.18009173,-86.50532033,24.01409822,
00056     -1.231739516,0.120858003e-2,-0.536382e-5};
00057   int j;
00058   x=xx-1.0;
00059   tmp=x+5.5;
00060   tmp -= (x+0.5)*log(tmp);
00061   ser=1.0;
00062   for (j=0;j<=5;j++)
00063   {
00064     x += 1.0;
00065     ser += cof[j]/x;
00066   }
00067   return -tmp+log(2.50662827465*ser);
00068 }
00069 
00077 dvector log_comb(const dvector& n, const dvector& r)
00078 {
00079   int mmin=n.indexmin();
00080   int mmax=n.indexmax();
00081   if (mmin != r.indexmin() || mmax != r.indexmax())
00082   {
00083     cerr << "Incompatible array bounds in function "
00084       "dvector log_comb(const dvector& n, const dvector& r)" << endl;
00085     ad_exit(1);
00086   }
00087   dvector tmp(mmin,mmax);
00088   for (int i=mmin;i<=mmax;i++)
00089   {
00090     tmp(i)=log_comb(n(i),r(i));
00091   }
00092   return tmp;
00093 }
00094 
00101 dvector log_comb(double n, const dvector& r)
00102 {
00103   int mmin=r.indexmin();
00104   int mmax=r.indexmax();
00105   dvector tmp(mmin,mmax);
00106   for (int i=mmin;i<=mmax;i++)
00107   {
00108     tmp(i)=log_comb(n,r(i));
00109   }
00110   return tmp;
00111 }
00112 
00119 dvector gammln(const dvector& v)
00120 {
00121   int mmin=v.indexmin();
00122   int mmax=v.indexmax();
00123   dvector tmp(mmin,mmax);
00124   for (int i=mmin;i<=mmax;i++)
00125   {
00126     tmp(i)=gammln(v(i));
00127   }
00128   return tmp;
00129 }
00130 
00137 dvector factln(const dvector& r)
00138 {
00139     int mmin=r.indexmin();
00140     int mmax=r.indexmax();
00141     dvector tmp(mmin,mmax);
00142     for (int i=mmin;i<=mmax;i++)
00143     {
00144       tmp(i)=gammln(r(i)+1.0);
00145     }
00146     return tmp;
00147 }