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(n-k);
51 52
}
......
87 88
} 
88 89

  
89 90
/**
90
 * Log-binomial of two arrays
91
 * \param n an array
91
 * Log-binomial 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.9843695780195716e-6, 
256
    1.5056327351493116e-7};
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.0-z);
289
  }
290
  else
291
  {
292
    return gammlnguts(z);
293
  }
294
}

Also available in: Unified diff