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) 2008-2012 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.0e-7
10
#define EPS 1.0e-9
11
#define FPMIN 1.0e-30
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.0-gammcf;
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.0-a;
58
  c=1.0/FPMIN;
59
  d=1.0/b;
60
  h=d;
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff