hyperg.c
Go to the documentation of this file.
1 /* hyperg.c
2  *
3  * Confluent hypergeometric function
4  *
5  *
6  *
7  * SYNOPSIS:
8  *
9  * double a, b, x, y, hyperg();
10  *
11  * y = hyperg( a, b, x );
12  *
13  *
14  *
15  * DESCRIPTION:
16  *
17  * Computes the confluent hypergeometric function
18  *
19  * 1 2
20  * a x a(a+1) x
21  * F ( a,b;x ) = 1 + ---- + --------- + ...
22  * 1 1 b 1! b(b+1) 2!
23  *
24  * Many higher transcendental functions are special cases of
25  * this power series.
26  *
27  * As is evident from the formula, b must not be a negative
28  * integer or zero unless a is an integer with 0 >= a > b.
29  *
30  * The routine attempts both a direct summation of the series
31  * and an asymptotic expansion. In each case error due to
32  * roundoff, cancellation, and nonconvergence is estimated.
33  * The result with smaller estimated error is returned.
34  *
35  *
36  *
37  * ACCURACY:
38  *
39  * Tested at random points (a, b, x), all three variables
40  * ranging from 0 to 30.
41  * Relative error:
42  * arithmetic domain # trials peak rms
43  * IEEE 0,30 30000 1.8e-14 1.1e-15
44  *
45  * Larger errors can be observed when b is near a negative
46  * integer or zero. Certain combinations of arguments yield
47  * serious cancellation error in the power series summation
48  * and also are not in the region of near convergence of the
49  * asymptotic series. An error message is printed if the
50  * self-estimated relative error is greater than 1.0e-12.
51  *
52  */
53 
54 /*
55  * Cephes Math Library Release 2.8: June, 2000
56  * Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
57  */
58 
59 #include "mconf.h"
60 #include <float.h>
61 
62 extern double MACHEP;
63 
64 
65 /* the `type` parameter determines what converging factor to use */
66 static double hyp2f0(double a, double b, double x, int type, double *err)
67 {
68  double a0, alast, t, tlast, maxt;
69  double n, an, bn, u, sum, temp;
70 
71  an = a;
72  bn = b;
73  a0 = 1.0e0;
74  alast = 1.0e0;
75  sum = 0.0;
76  n = 1.0e0;
77  t = 1.0e0;
78  tlast = 1.0e9;
79  maxt = 0.0;
80 
81  do {
82  if (an == 0)
83  goto pdone;
84  if (bn == 0)
85  goto pdone;
86 
87  u = an * (bn * x / n);
88 
89  /* check for blowup */
90  temp = fabs(u);
91  if ((temp > 1.0) && (maxt > (DBL_MAX / temp)))
92  goto error;
93 
94  a0 *= u;
95  t = fabs(a0);
96 
97  /* terminating condition for asymptotic series:
98  * the series is divergent (if a or b is not a negative integer),
99  * but its leading part can be used as an asymptotic expansion
100  */
101  if (t > tlast)
102  goto ndone;
103 
104  tlast = t;
105  sum += alast; /* the sum is one term behind */
106  alast = a0;
107 
108  if (n > 200)
109  goto ndone;
110 
111  an += 1.0e0;
112  bn += 1.0e0;
113  n += 1.0e0;
114  if (t > maxt)
115  maxt = t;
116  }
117  while (t > MACHEP);
118 
119 
120  pdone: /* series converged! */
121 
122  /* estimate error due to roundoff and cancellation */
123  *err = fabs(MACHEP * (n + maxt));
124 
125  alast = a0;
126  goto done;
127 
128  ndone: /* series did not converge */
129 
130  /* The following "Converging factors" are supposed to improve accuracy,
131  * but do not actually seem to accomplish very much. */
132 
133  n -= 1.0;
134  x = 1.0 / x;
135 
136  switch (type) { /* "type" given as subroutine argument */
137  case 1:
138  alast *=
139  (0.5 + (0.125 + 0.25 * b - 0.5 * a + 0.25 * x - 0.25 * n) / x);
140  break;
141 
142  case 2:
143  alast *= 2.0 / 3.0 - b + 2.0 * a + x - n;
144  break;
145 
146  default:
147  ;
148  }
149 
150  /* estimate error due to roundoff, cancellation, and nonconvergence */
151  *err = MACHEP * (n + maxt) + fabs(a0);
152 
153  done:
154  sum += alast;
155  return (sum);
156 
157  /* series blew up: */
158  error:
159  *err = INFINITY;
160  sf_error("hyperg", SF_ERROR_NO_RESULT, NULL);
161  return (sum);
162 }
163 
164 
165 /* asymptotic formula for hypergeometric function:
166  *
167  * ( -a
168  * -- ( |z|
169  * | (b) ( -------- 2f0( a, 1+a-b, -1/x )
170  * ( --
171  * ( | (b-a)
172  *
173  *
174  * x a-b )
175  * e |x| )
176  * + -------- 2f0( b-a, 1-a, 1/x ) )
177  * -- )
178  * | (a) )
179  */
180 
181 static double hy1f1a(double a, double b, double x, double *err)
182 {
183  double h1, h2, t, u, temp, acanc, asum, err1, err2;
184 
185  if (x == 0) {
186  acanc = 1.0;
187  asum = INFINITY;
188  goto adone;
189  }
190  temp = log(fabs(x));
191  t = x + temp * (a - b);
192  u = -temp * a;
193 
194  if (b > 0) {
195  temp = lgam(b);
196  t += temp;
197  u += temp;
198  }
199 
200  h1 = hyp2f0(a, a - b + 1, -1.0 / x, 1, &err1);
201 
202  temp = exp(u) / gamma(b - a);
203  h1 *= temp;
204  err1 *= temp;
205 
206  h2 = hyp2f0(b - a, 1.0 - a, 1.0 / x, 2, &err2);
207 
208  if (a < 0)
209  temp = exp(t) / gamma(a);
210  else
211  temp = exp(t - lgam(a));
212 
213  h2 *= temp;
214  err2 *= temp;
215 
216  if (x < 0.0)
217  asum = h1;
218  else
219  asum = h2;
220 
221  acanc = fabs(err1) + fabs(err2);
222 
223  if (b < 0) {
224  temp = gamma(b);
225  asum *= temp;
226  acanc *= fabs(temp);
227  }
228 
229 
230  if (asum != 0.0)
231  acanc /= fabs(asum);
232 
233  if (acanc != acanc)
234  /* nan */
235  acanc = 1.0;
236 
237  if (asum == INFINITY || asum == -INFINITY)
238  /* infinity */
239  acanc = 0;
240 
241  acanc *= 30.0; /* fudge factor, since error of asymptotic formula
242  * often seems this much larger than advertised */
243 
244  adone:
245  *err = acanc;
246  return (asum);
247 }
248 
249 
250 /* Power series summation for confluent hypergeometric function */
251 static double hy1f1p(double a, double b, double x, double *err)
252 {
253  double n, a0, sum, t, u, temp, maxn;
254  double an, bn, maxt;
255  double y, c, sumc;
256 
257 
258  /* set up for power series summation */
259  an = a;
260  bn = b;
261  a0 = 1.0;
262  sum = 1.0;
263  c = 0.0;
264  n = 1.0;
265  t = 1.0;
266  maxt = 0.0;
267  *err = 1.0;
268 
269  maxn = 200.0 + 2 * fabs(a) + 2 * fabs(b);
270 
271  while (t > MACHEP) {
272  if (bn == 0) { /* check bn first since if both */
273  sf_error("hyperg", SF_ERROR_SINGULAR, NULL);
274  return (INFINITY); /* an and bn are zero it is */
275  }
276  if (an == 0) /* a singularity */
277  return (sum);
278  if (n > maxn) {
279  /* too many terms; take the last one as error estimate */
280  c = fabs(c) + fabs(t) * 50.0;
281  goto pdone;
282  }
283  u = x * (an / (bn * n));
284 
285  /* check for blowup */
286  temp = fabs(u);
287  if ((temp > 1.0) && (maxt > (DBL_MAX / temp))) {
288  *err = 1.0; /* blowup: estimate 100% error */
289  return sum;
290  }
291 
292  a0 *= u;
293 
294  y = a0 - c;
295  sumc = sum + y;
296  c = (sumc - sum) - y;
297  sum = sumc;
298 
299  t = fabs(a0);
300 
301  an += 1.0;
302  bn += 1.0;
303  n += 1.0;
304  }
305 
306  pdone:
307 
308  /* estimate error due to roundoff and cancellation */
309  if (sum != 0.0) {
310  *err = fabs(c / sum);
311  }
312  else {
313  *err = fabs(c);
314  }
315 
316  if (*err != *err) {
317  /* nan */
318  *err = 1.0;
319  }
320 
321  return (sum);
322 }
323 
324 
325 
326 double hyperg(double a, double b, double x)
327 {
328  double asum, psum, acanc, pcanc, temp;
329 
330  /* See if a Kummer transformation will help */
331  temp = b - a;
332  if (fabs(temp) < 0.001 * fabs(a))
333  return (exp(x) * hyperg(temp, b, -x));
334 
335 
336  /* Try power & asymptotic series, starting from the one that is likely OK */
337  if (fabs(x) < 10 + fabs(a) + fabs(b)) {
338  psum = hy1f1p(a, b, x, &pcanc);
339  if (pcanc < 1.0e-15)
340  goto done;
341  asum = hy1f1a(a, b, x, &acanc);
342  }
343  else {
344  psum = hy1f1a(a, b, x, &pcanc);
345  if (pcanc < 1.0e-15)
346  goto done;
347  asum = hy1f1p(a, b, x, &acanc);
348  }
349 
350  /* Pick the result with less estimated error */
351 
352  if (acanc < pcanc) {
353  pcanc = acanc;
354  psum = asum;
355  }
356 
357  done:
358  if (pcanc > 1.0e-12)
359  sf_error("hyperg", SF_ERROR_LOSS, NULL);
360 
361  return (psum);
362 }
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
asum
RealScalar EIGEN_BLAS_FUNC() asum(int *n, RealScalar *px, int *incx)
Definition: level1_real_impl.h:14
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
b
Scalar * b
Definition: benchVecAdd.cpp:17
SF_ERROR_NO_RESULT
@ SF_ERROR_NO_RESULT
Definition: sf_error.h:15
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
type
Definition: pytypes.h:1525
log
const EIGEN_DEVICE_FUNC LogReturnType log() const
Definition: ArrayCwiseUnaryOps.h:128
exp
const EIGEN_DEVICE_FUNC ExpReturnType exp() const
Definition: ArrayCwiseUnaryOps.h:97
lgam
double lgam(double x)
Definition: gamma.c:275
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:119
n
int n
Definition: BiCGSTAB_simple.cpp:1
hy1f1a
static double hy1f1a(double a, double b, double x, double *err)
Definition: hyperg.c:181
hy1f1p
static double hy1f1p(double a, double b, double x, double *err)
Definition: hyperg.c:251
gamma
#define gamma
Definition: mconf.h:85
MACHEP
double MACHEP
Definition: const.c:54
y
Scalar * y
Definition: level1_cplx_impl.h:124
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
mconf.h
error
static double error
Definition: testRot3.cpp:37
hyp2f0
static double hyp2f0(double a, double b, double x, int type, double *err)
Definition: hyperg.c:66
sf_error
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
Definition: sf_error.c:41
SF_ERROR_SINGULAR
@ SF_ERROR_SINGULAR
Definition: sf_error.h:10
NULL
#define NULL
Definition: ccolamd.c:609
hyperg
double hyperg(double a, double b, double x)
Definition: hyperg.c:326
align_3::t
Point2 t(10, 10)
SF_ERROR_LOSS
@ SF_ERROR_LOSS
Definition: sf_error.h:14


gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:02:26