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 
* logfactorial, logbinomial, and loggamma 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(nk); 
17  29 
} 
18  30  
31 
/** 

32 
* Logfactorial \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.120858003e2,0.536382e5}; 

40 
int j; 

41 
x=xx1.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 
* Logbinomial 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 " 
Also available in: Unified diff