Revision 601 trunk/src/linad99/dfgammp.cpp
dfgammp.cpp (revision 601)  

1 
/* 

1 
/**


2  2 
* $Id$ 
3  3 
* 
4  4 
* Author: David Fournier 
...  ...  
19  19 
#else 
20  20 
# include <fvar.hpp> 
21  21 
#endif 
22 
#include <admodel.h>


23 
#include <df11fun.h>


24 
#include <df12fun.h>


25  
22 
#define ITMAX 100


23 
#define EPS 1.0e9


24 
//#define EPS 3.0e7


25 
#define FPMIN 1.0e30 

26  26 
double get_values(double x,double y,int print_switch); 
27 
dvariable igam(const dvariable & _a, const dvariable & _x); 

28 
dvariable igamc(const dvariable & _a, const dvariable & _x); 

29 
dvariable lgam(const prevariable& v1); 

30 
dvariable private_lgam(const dvariable& v); 

31 
df3_one_variable lgam(const df3_one_variable& _x); 

32  27  
33 
extern int mtherr(char* s,int n); 

34  
35 
namespace Cephes 

28 
/** Incomplete gamma function. 

29 
Continued fraction approximation. 

30 
\n\n The implementation of this algorithm was inspired by 

31 
"Numerical Recipes in C", 2nd edition, 

32 
Press, Teukolsky, Vetterling, Flannery, chapter 6 

33 
*/ 

34 
void gcf(const dvariable& _gammcf,const dvariable& a, 

35 
const dvariable& x,const dvariable& _gln) 

36  36 
{ 
37 
extern const double A[]; 

38 
extern const double B[]; 

39 
extern const double C[]; 

40 
extern const double Q[]; 

41 
extern const double P[]; 

42 
extern const double STIR[]; 

43 
extern int sgngam; 

44 
extern const double MAXLOG; 

45 
extern const double MAXNUM; 

46 
extern const double LOGPI; 

47 
extern const double big; 

48 
extern const double biginv; 

49 
extern const double MACHEP; 

50 
extern const double MYINF; 

51 
extern const double SQTPI; 

52 
extern const double MAXLGM; 

53 
extern const double LS2PI; 

54 
extern const double MAXSTIR; 

37 
ADUNCONST(dvariable,gln) 

38 
ADUNCONST(dvariable,gammcf) 

39 
int i; 

40 
dvariable an,b,c,d,del,h; 

55  41  
56 
dvariable polevl(const dvariable & x, const void *_coef, int N); 

57 
dvariable p1evl(const dvariable & x, const void *_coef, int N); 

58  
59 
/** 

60 
* \ingroup gammafunc 

61 
* Stirling's formula (approximation to large factorials) 

62 
* \param _x \f$x\f$ 

63 
* \return Sterling's approximation to \f$x!\f$ 

64 
* 

65 
* \n\n Cephes Math Library Release 2.1: December, 1988 

66 
* Copyright 1984, 1987, 1988 by Stephen L. Moshier 

67 
* Direct inquiries to 30 Frost Street, Cambridge, MA 02140 

68 
*/ 

69 
dvariable stirf(const dvariable & _x) 

70 
{ 

71 
ADUNCONST(dvariable, x) dvariable y, w, v; 

72  
73 
w = 1.0 / x; 

74 
w = 1.0 + w * polevl(w, STIR, 4); 

75 
y = exp(x); 

76 
if (value(x) > MAXSTIR) 

77 
{ /* Avoid overflow in pow() */ 

78 
v = pow(x, 0.5 * x  0.25); 

79 
y = v * (v / y); 

80 
} else 

81 
{ 

82 
y = pow(x, x  0.5) / y; 

83 
} 

84 
y = SQTPI * y * w; 

85 
return (y); 

86 
} 

87  
88 
/** 

89 
* \ingroup gammafunc 

90 
* Polynomial evaluation 

91 
* \param x \f$x\f$ the point to be evaluated 

92 
* \param _coef The coefficents of the polynomial 
Also available in: Unified diff