Revision 601 trunk/src/df1b2separable/df3gammp.cpp
df3gammp.cpp (revision 601)  

1 
/* 

1 
/**


2  2 
* $Id$ 
3  3 
* 
4  4 
* Author: David Fournier 
5  5 
* Copyright (c) 20092012 ADMB Foundation 
6  6 
*/ 
7 
/** 

8 
* \file 

9 
* Description not yet available. 

10 
*/ 

11  7 
#include <df1b2fun.h> 
12 
#include <df13fun.h> 

13 
#include "mconf.h" 

8 
#define ITMAX 100 

9 
#define EPS 1.0e9 

10 
//#define EPS 3.0e7 

11 
#define FPMIN 1.0e30 

12 
double get_values(double x,double y,int print_switch); 

14  13  
15 
#ifdef NANS 

16 
#undef NANS 

17 
#endif 

18 
#ifdef INFINITIES 

19 
#undef INFINITIES 

20 
#endif 

21  14  
22 
double get_values(double x, double y, int print_switch); 

23 
df3_two_variable igam(const df3_two_variable & _a, 

24 
const df3_two_variable & _x); 

25 
df3_two_variable igamc(const df3_two_variable & _a, 

26 
const df3_two_variable & _x); 

27  
28 
namespace Cephes 

15 
df1b2variable log_negbinomial_density(double x,const df1b2variable& _xmu, 

16 
const df1b2variable& _xtau) 

29  17 
{ 
18 
ADUNCONST(df1b2variable,xmu) 

19 
ADUNCONST(df1b2variable,xtau) 

20 
init_df3_two_variable mu(xmu); 

21 
init_df3_two_variable tau(xtau); 

22 
*mu.get_u_x()=1.0; 

23 
*tau.get_u_y()=1.0; 

24 
if (value(tau)1.0<0.0) 

25 
{ 

26 
cerr << "tau <=1 in log_negbinomial_density " << endl; 

27 
ad_exit(1); 

28 
} 

29 
df3_two_variable r=mu/(1.e120+(tau1.0)); 

30 
df3_two_variable tmp; 

31 
tmp=gammln(x+r)gammln(r) gammln(x+1) 

32 
+r*log(r)+x*log(mu)(r+x)*log(r+mu); 

33 
df1b2variable tmp1; 

34 
tmp1=tmp; 

35 
return tmp1; 

36 
} 

30  37  
31 
int sgngam = 0; 

32 
extern const double MAXLOG = 200; 

33 
extern const double MAXNUM = 1.7976931348623158E+308; 

34 
extern const double LOGPI = 1.14472988584940017414; 

35 
extern const double big = 4.503599627370496e15; 

36 
extern const double biginv = 2.22044604925031308085e16; 

37 
extern const double MACHEP = 2.22045e16; 

38 
extern const double MYINF = 1.7976931348623158E+308; 

38 
/** Log gamma function. 

39 
\param xx \f$x\f$ 

40 
\return \f$\ln\bigr(\Gamma(x)\bigl)\f$ 

39  41  
40 
df3_two_variable polevl(const df3_two_variable & x, const double *_coef, int N); 

41 
df3_two_variable p1evl(const df3_two_variable & x, const double *_coef, int N); 

42  
43 
/* Coefficents used for the rational function in gamma(x) */ 

44 
#ifdef UNK 

45 
extern const double P[] = { 

46 
1.60119522476751861407E4, 

47 
1.19135147006586384913E3, 

48 
1.04213797561761569935E2, 

49 
4.76367800457137231464E2, 

50 
2.07448227648435975150E1, 

51 
4.94214826801497100753E1, 

52 
9.99999999999999996796E1 

53 
}; 

54  
55 
extern const double Q[] = { 

56 
2.31581873324120129819E5, 

57 
5.39605580493303397842E4, 

58 
4.45641913851797240494E3, 

59 
1.18139785222060435552E2, 

60 
3.58236398605498653373E2, 

61 
2.34591795718243348568E1, 

62 
7.14304917030273074085E2, 

63 
1.00000000000000000320E0 

64 
}; 

65  
66 
#define MAXGAM 171.624376956302725 
Also available in: Unified diff