Gammaq.cxx
Go to the documentation of this file.
1 
10 #include "Gammaq.h"
11 
12 #include <iostream>
13 
14 #include <cassert>
15 #include <cmath>
16 
17 using std::cerr;
18 using std::endl;
19 using std::abs;
20 
21 namespace hippodraw {
22  namespace Numeric {
23 
24 double gammln( double xx )
25 {
26  double x, y, tmp, ser;
27  static double cof[6]={76.18009172947146,-86.50532032941677,
28  24.01409824083091,-1.231739572450155,
29  0.1208650973866179e-2,-0.5395239384953e-5};
30 
31  y = x = xx;
32  tmp = x + 5.5;
33  tmp -= ( x + 0.5 ) * log( tmp );
34  ser = 1.000000000190015;
35 
36  for ( int j = 0; j <= 5; j++)
37  ser += cof[j] / (++y);
38 
39  return -tmp + log( 2.5066282746310005 * ser / x);
40 }
41 
42 
43 void gser(double *gamser, double a, double x, double *gln)
44 {
45  int n;
46  double sum, del, ap;
47 
48  *gln = gammln( a );
49 
50  if ( x <= 0.0 )
51  {
52  if (x < 0.0) cerr << "Error: x less than 0 in routine gser" << endl;
53  *gamser=0.0;
54  return;
55  }
56  else
57  {
58  ap = a;
59  del = sum =1.0 / a;
60  for (n = 1; n <= ITMAX; n++)
61  {
62  ++ap;
63  del *= x/ap;
64  sum += del;
65 
66  if (abs(del) < abs(sum)*EPS)
67  {
68  *gamser = sum * exp( -x + a * log( x ) -( *gln ) );
69  return;
70  }
71  }
72 
73  cerr << "ERROR: a too large, ITMAX too small in routine gser" << endl;
74  return;
75  }
76 }
77 
78 void gcf( double *gammcf, double a, double x, double *gln )
79 {
80  int i;
81  double an,b,c,d,del,h;
82 
83  *gln = gammln( a );
84  b = x + 1.0 - a;
85  c = 1.0 / FPMIN;
86  d = 1.0 / b;
87  h = d;
88 
89  for( i = 1; i <= ITMAX; i++ )
90  {
91  an = -i*(i-a);
92  b += 2.0;
93  d=an*d+b;
94 
95  if ( abs( d ) < FPMIN )
96  d = FPMIN;
97 
98  c = b + an / c;
99 
100  if ( abs( c ) < FPMIN )
101  c = FPMIN;
102 
103  d = 1.0 / d;
104  del = d * c;
105  h *= del;
106 
107  if ( abs( del - 1.0 ) < EPS)
108  break;
109  }
110 
111  if (i > ITMAX) cerr << "a too large, ITMAX too small in gcf" << endl;
112 
113  *gammcf = exp( -x + a * log( x ) -( *gln ) ) * h;
114 }
115 
116 double gammq(double a, double x)
117 {
118  void gcf(double *gammcf, double a, double x, double *gln);
119 
120  double gamser,gammcf,gln;
121 
122  if (x < 0.0 || a <= 0.0)
123  {
124  cerr << "a = " << a << " x = " << x << endl;
125  cerr << "Invalid arguments in routine GAMMQ" << endl;
126  assert(0);
127  }
128 
129  if (x < (a+1.0))
130  {
131  gser(&gamser,a,x,&gln);
132  return 1.0-gamser;
133  }
134  else
135  {
136  gcf(&gammcf,a,x,&gln);
137  return gammcf;
138  }
139 }
140 
141 } // end namespace Numeric
142 } // end namespace hippodraw
143 
150 /*
151 int main(void)
152 {
153  char txt[MAXSTR];
154  int i,nval;
155  float a,val,x;
156  FILE *fp;
157 
158  if ((fp = fopen("fncval.dat","r")) == NULL)
159  printf("Data file fncval.dat not found\n");
160 
161  std::fgets(txt,MAXSTR,fp);
162 
163  while ( std::strncmp(txt,"Incomplete Gamma Function",25))
164  {
165  std::fgets(txt,MAXSTR,fp);
166  if ( std::feof(fp))
167  std::printf("Data not found in fncval.dat\n");
168  }
169 
170  std::fscanf(fp,"%d %*s",&nval);
171  std::printf("\n%s\n",txt);
172  std::printf("%4s %11s %14s %14s \n","a","x","actual","gammq(a,x)");
173 
174  for (i=1;i<=nval;i++)
175  {
176  std::fscanf(fp,"%f %f %f",&a,&x,&val);
177  std::printf("%6.2f %12.6f %12.6f %12.6f\n", a, x, (1.0-val), gammq(a,x));
178  }
179 
180  std::fclose(fp);
181  return 0;
182 }
183 
184 */
185 
186 #undef ITMAX
187 #undef EPS
188 #undef FPMIN
189 #undef MAXSTR

Generated for HippoDraw Class Library by doxygen