Revision 795 branches/replacement/src/linad99/combc.cpp

combc.cpp (revision 795)
4 4
 * Author: David Fournier
5 5
 * Copyright (c) 2009 ADMB Foundation
6 6
 */
7
/**
8
 * \file
9
 * This file has routines for the
10
 * log-factorial, log-binomial, and log-gamma functions
11
 */
12

  
7 13
#include <fvar.hpp>
8 14

  
9

  
10 15
double factln(double n);
11 16
double gammln(double xx);
12 17
double lgam(const double xx);
13 18

  
19
/**
20
 * Log of the binomial coefficent
21
 * i.e log of 'n choose k'
22
 * \param n a number
23
 * \param k a number
24
 * \return log of the binomial coefficent
25
 */
14 26
double log_comb(double n, double k)
15 27
{
16 28
  return factln(n)-factln(k)-factln(n-k);
17 29
}
18 30

  
31
/**
32
 * Log-factorial \f$\ln(n!)\f$
33
 * \param n a number
34
 * \return log of the factorial
35
 */
19 36
double factln(double n)
20 37
{
21 38
  return gammln(n+1.0);
22 39
}
23 40

  
24
/* Log gamma function.
25
    \param xx \f$x\f$
26
    \return \f$\ln\bigr(\Gamma(x)\bigl)\f$
27

  
28
    \n\n The implementation of this algorithm was inspired by
29
    "Numerical Recipes in C", 2nd edition,
30
    Press, Teukolsky, Vetterling, Flannery, chapter 6
31

  
32
    \deprecated Scheduled for replacement by 2010.
33
*/
34
/*
41
/**
42
 * A wrapper for igam
43
 */
35 44
double gammln(double xx)
36 45
{
37
  double x,tmp,ser;
38
  static double cof[6]={76.18009173,-86.50532033,24.01409822,
39
    -1.231739516,0.120858003e-2,-0.536382e-5};
40
  int j;
41
  x=xx-1.0;
42
  tmp=x+5.5;
43
  tmp -= (x+0.5)*log(tmp);
44
  ser=1.0;
45
  for (j=0;j<=5;j++) 
46
  {
47
    x += 1.0;
48
    ser += cof[j]/x;
49
  }
50
  return -tmp+log(2.50662827465*ser);
46
  return lgam(xx);
51 47
}
52
*/
53 48

  
54
double gammln(double xx)
49
/**
50
 * Log-binomial of two arrays
51
 * \param n an array
52
 * \param r an array
53
 * \return log of the binomial coefficent
54
 */
55
dvector log_comb(_CONST dvector& n,_CONST dvector& r)
55 56
{
56
  return lgam(xx);
57
}
58

  
59

  
60

  
61
  dvector log_comb(_CONST dvector& n,_CONST dvector& r)
57
  int mmin=n.indexmin();
58
  int mmax=n.indexmax();
59
  if (mmin != r.indexmin() || mmax != r.indexmax()) 
62 60
  {
63
    int mmin=n.indexmin();
64
    int mmax=n.indexmax();
65
    if (mmin != r.indexmin() || mmax != r.indexmax()) 
66
    {
67
      cerr << "Incompatible array bounds in function "
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff