Revision 601 trunk/src/linad99/combv.cpp
combv.cpp (revision 601)  

1 
/* 

1 
/**


2  2 
* $Id$ 
3  3 
* Author: Unknown 
4  4 
*/ 
...  ...  
45  45 
* \param n a number 
46  46 
* \param k a number 
47  47 
* \return log of the binomial coefficent 
48 
*/dvariable log_comb(double n, const dvariable& k) 

48 
*/ 

49 
dvariable log_comb(double n, const dvariable& k) 

49  50 
{ 
50  51 
return factln(n)factln(k)factln(nk); 
51  52 
} 
...  ...  
87  88 
} 
88  89  
89  90 
/** 
90 
* Logbinomial of two arrays


91 
* \param n an array


91 
* Logbinomial of a number with an arrays


92 
* \param n a number


92  93 
* \param r an array 
93  94 
* \return log of the binomial coefficent 
94  95 
*/ 
...  ...  
235  236 
RETURN_ARRAYS_DECREMENT(); 
236  237 
return tmp; 
237  238 
} 
239  
240 
static dvariable gammlnguts(const prevariable _z) 

241 
{ 

242 
double z = value(_z); 

243 
double zdot=1.0; 

244 
const double lpi =1.1447298858494001741434272; 

245 
const double pi =3.1415926535897932384626432; 

246 
const double lpp =0.9189385332046727417803297; 

247 
int n=7; 

248 
const double c[9]={0.99999999999980993, 

249 
676.5203681218851, 

250 
1259.1392167224028, 

251 
771.32342877765313, 

252 
176.61502916214059, 

253 
12.507343278686905, 

254 
0.13857109526572012, 

255 
9.9843695780195716e6, 

256 
1.5056327351493116e7}; 

257 
z=1.0; 

258 
double x=c[0]; 

259 
double xdot=0.0; 

260 
for (int i=1;i<=n+1;i++) 

261 
{ 

262 
double zinv=1.0/(z+i); 

263 
x+=c[i]*zinv; 

264 
//xdot=c[i]/square(z+i)*zdot; since zdot=1.0 

265 
xdot=c[i]*square(zinv); 

266 
} 

267 
double t=z+n+0.5; 

268 
double tdot=zdot; 

269 
//return lpp + (z+0.5)*log(t) t + log(x); 

270 
double ans= lpp + (z+0.5)*log(t) t + log(x); 

271 
//double ansdot=zdot*log(t) + (z+0.5)/t*tdot tdot +xdot/x; 

272 
// since tdot=1.0 

273 
// since zdot=1.0 

274 
double ansdot=log(t) + (z+0.5)/t 1.0 +xdot/x; 

275 
dvariable u; 

276 
u.v>x=ans; 

277 
gradient_structure::GRAD_STACK1>set_gradient_stack(default_evaluation, 

278 
&(u.v>x), &(_z.v>x), ansdot ); 

279 
return(u); 

280 
} 

281  
282 
dvariable gammln(const prevariable& z) 

283 
{ 

284 
const double lpi =1.1447298858494001741434272; 

285 
const double pi =3.1415926535897932384626432; 

286 
if (z<0.5) 

287 
{ 

288 
return lpi  log(sin(pi*z))  gammlnguts(1.0z); 

289 
} 

290 
else 

291 
{ 

292 
return gammlnguts(z); 

293 
} 

294 
} 
Also available in: Unified diff