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.e-9
19
//#define EPS 3.0e-7
20
#define FPMIN 1.0e-30
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.0-gammcf;
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
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff