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