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

1 
/* 

1 
/**


2  2 
* $Id$ 
3  3 
* 
4  4 
* Author: David Fournier 
...  ...  
14  14 
*/ 
15  15  
16  16 
#include <fvar.hpp> 
17 
#define ITMAX 100 

18 
#define EPS 1.e9 

19 
//#define EPS 3.0e7 

20 
#define FPMIN 1.0e30 

21 
void gcf(double& gammcf,double a,double x,double &gln); 

22 
void gser(double& gamser,double a,double x,double& gln); 

17  23  
18 
#ifdef NANS 

19 
#undef NANS 

20 
#endif 

21 
#ifdef INFINITIES 

22 
#undef INFINITIES 

23 
#endif 

24  
25 
namespace Cephes 

24 
double gammp(double a,double x) 

26  25 
{ 
27 
extern const double A[]; 

28 
extern const double B[]; 

29 
extern const double C[]; 

30 
extern int sgngam; 

31 
extern const double MAXLOG; 

32 
extern const double MAXNUM; 

33 
#ifndef PI 

34 
extern const double PI; 

35 
#endif 

36 
extern const double LOGPI; 

37 
extern const double big; 

38 
extern const double biginv; 

39 
extern const double MACHEP; 

40 
extern const double MYINF; 

41 
extern const double SQTPI; 

42 
extern const double MAXLGM; 

43 
extern const double LS2PI; 

44 
extern const double MAXSTIR; 

26 
double gamser,gammcf,gln; 

45  27  
46 
double polevl(double x, const double *_coef, int N); 

47 
double p1evl(double x, const double *_coef, int N); 

48 
/** 

49 
* \ingroup gammafunc 

50 
* Polynomial evaluation 

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

52 
* \param _coef The coefficents of the polynomial 

53 
* \param N \f$N\f$ The degree of the polynomial 

54 
* \return The polynomial evaluated at \f$x\f$ 

55 
* 

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

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

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

59 
*/ 

60 
double polevl(double x, const double *_coef, int N) 

61 
{ 

62 
double *coef = (double *) (_coef); 

63 
double ans; 

64 
int i; 

65 
double *p; 

28 
if (x < 0.0  a <= 0.0) 

29 
cerr << "Invalid arguments in routine gammp" << endl; 

30 
if (x < (a+1.0)) { 

31 
gser(gamser,a,x,gln); 

32 
return gamser; 

33 
} else { 

34 
gcf(gammcf,a,x,gln); 

35 
return 1.0gammcf; 

36 
} 

37 
} 

66  38  
67 
p = coef;


68 
ans = *p++;


69 
i = N;


39 
double cumd_gamma(double x,double a)


40 
{


41 
double gamser,gammcf,gln;


70  42  
71 
do 

72 
ans = ans * x + *p++; 

73 
while (i); 

74  
75 
return (ans); 

76 
} 

77  
78 
/** 

79 
* \ingroup gammafunc 

80 
* Polynomial evaluation when leading coefficent is 1 

81 
* (i.e. leading term is \f$x^N\f$) 

82 
* \param x \f$x\f$ the point to be evaluated 
Also available in: Unified diff