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

1 
/* 

1 
/**


2  2 
* $Id$ 
3  3 
* 
4  4 
* Author: David Fournier 
5 
* Copyright (c) 20082012 Regents of the University of California


5 
* Copyright (c) 2008, 2009, 2010 Regents of the University of California


6  6 
*/ 
7 
/** 

8 
* \file 

9 
* This file deals with the Incomplete Gamma Functions 

10 
* of constant types. All supporting mathematical functions 

11 
* required to compute the Inmomplete Gamma Function 

12 
* are included. They being: log gamma, 

13 
* and some polynomial evaluation functions. 

14 
*/ 

15  
16  7 
#include <fvar.hpp> 
8 
#define ITMAX 100 

9 
//#define EPS 3.0e7 

10 
#define EPS 1.0e9 

11 
#define FPMIN 1.0e30 

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

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

17  14  
18 
double igam(const double &a, const double &x); 

19 
double igamc(const double &a, const double &x); 

20 
double lgam(double x); 

21 
int mtherr(char *s, int n) 

22 
{ 

23 
return 0; 

24 
} 

15 
dvariable gamma_deviate(const prevariable& _x,const prevariable& _a) 

16 
{ 

17 
prevariable& x= (prevariable&)(_x); 

18 
prevariable& a= (prevariable&)(_a); 

25  19  
26 
/** 

27 
* Gamma Deviate 

28 
* \param _x 

29 
* \param _a 

30 
*/ 

31 
dvariable gamma_deviate(const prevariable & _x, const prevariable & _a) 

32 
{ 

33 
prevariable & x = (prevariable &) (_x); 

34 
prevariable & a = (prevariable &) (_a); 

20 
dvariable y=cumd_norm(x); 

35  21  
36 
dvariable y = cumd_norm(x);


22 
y=.9999*y+.00005;


37  23  
38 
y = .9999 * y + .00005;


24 
dvariable z=inv_cumd_gamma(y,a);


39  25  
40 
dvariable z = inv_cumd_gamma(y, a); 

26 
return z; 

27 
} 

41  28  
42 
return z; 

29  
30 
static double gammp(double a,double x) 

31 
{ 

32 
double gamser,gammcf,gln; 

33  
34 
if (x < 0.0  a <= 0.0) 

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

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

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

38 
return gamser; 

39 
} else { 

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

41 
return 1.0gammcf; 

42 
} 

43  43 
} 
44  44  
45 
/** 

46 
* A wrapper for igam 

47 
*/ 

48 
static double gammp(double a, double x) 

45 
/** Incomplete gamma function. 

46 
Continued fraction approximation. 

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

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

49 
Press, Teukolsky, Vetterling, Flannery, chapter 6 

50 
*/ 

51 
static void gcf(double& gammcf,double a,double x,double &gln) 

49  52 
{ 
50 
return igam(a, x); 

53 
int i; 

54 
double an,b,c,d,del,h; 

55  
56 
gln=gammln(a); 

57 
b=x+1.0a; 

58 
c=1.0/FPMIN; 

59 
d=1.0/b; 

60 
h=d; 
Also available in: Unified diff