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

Also available in: Unified diff