source: sasview/src/sas/models/c_extension/cephes/hyperg.c @ c81f9b6

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since c81f9b6 was 79492222, checked in by krzywon, 10 years ago

Changed the file and folder names to remove all SANS references.

  • Property mode set to 100644
File size: 7.1 KB
Line 
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 *    DEC       0,30         2000       1.2e-15     1.3e-16
44 qtst1:
45 21800   max =  1.4200E-14   rms =  1.0841E-15  ave = -5.3640E-17
46 ltstd:
47 25500   max = 1.2759e-14   rms = 3.7155e-16  ave = 1.5384e-18
48 *    IEEE      0,30        30000       1.8e-14     1.1e-15
49 *
50 * Larger errors can be observed when b is near a negative
51 * integer or zero.  Certain combinations of arguments yield
52 * serious cancellation error in the power series summation
53 * and also are not in the region of near convergence of the
54 * asymptotic series.  An error message is printed if the
55 * self-estimated relative error is greater than 1.0e-12.
56 *
57 */
58
59/*                                                      hyperg.c */
60
61
62/*
63Cephes Math Library Release 2.8:  June, 2000
64Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
65*/
66
67#include "mconf.h"
68
69#ifdef ANSIPROT
70extern double exp ( double );
71extern double log ( double );
72extern double gamma ( double );
73extern double lgam ( double );
74extern double fabs ( double );
75double hyp2f0 ( double, double, double, int, double * );
76static double hy1f1p(double, double, double, double *);
77static double hy1f1a(double, double, double, double *);
78double hyperg (double, double, double);
79#else
80double exp(), log(), gamma(), lgam(), fabs(), hyp2f0();
81static double hy1f1p();
82static double hy1f1a();
83double hyperg();
84#endif
85extern double MAXNUM, MACHEP;
86
87double hyperg( a, b, x)
88double a, b, x;
89{
90double asum, psum, acanc, pcanc, temp;
91
92/* See if a Kummer transformation will help */
93temp = b - a;
94if( fabs(temp) < 0.001 * fabs(a) )
95        return( exp(x) * hyperg( temp, b, -x )  );
96
97
98psum = hy1f1p( a, b, x, &pcanc );
99if( pcanc < 1.0e-15 )
100        goto done;
101
102
103/* try asymptotic series */
104
105asum = hy1f1a( a, b, x, &acanc );
106
107
108/* Pick the result with less estimated error */
109
110if( acanc < pcanc )
111        {
112        pcanc = acanc;
113        psum = asum;
114        }
115
116done:
117if( pcanc > 1.0e-12 )
118        mtherr( "hyperg", PLOSS );
119
120return( psum );
121}
122
123
124
125
126/* Power series summation for confluent hypergeometric function         */
127
128
129static double hy1f1p( a, b, x, err )
130double a, b, x;
131double *err;
132{
133double n, a0, sum, t, u, temp;
134double an, bn, maxt, pcanc;
135
136
137/* set up for power series summation */
138an = a;
139bn = b;
140a0 = 1.0;
141sum = 1.0;
142n = 1.0;
143t = 1.0;
144maxt = 0.0;
145
146
147while( t > MACHEP )
148        {
149        if( bn == 0 )                   /* check bn first since if both */
150                {
151                mtherr( "hyperg", SING );
152                return( MAXNUM );       /* an and bn are zero it is     */
153                }
154        if( an == 0 )                   /* a singularity                */
155                return( sum );
156        if( n > 200 )
157                goto pdone;
158        u = x * ( an / (bn * n) );
159
160        /* check for blowup */
161        temp = fabs(u);
162        if( (temp > 1.0 ) && (maxt > (MAXNUM/temp)) )
163                {
164                pcanc = 1.0;    /* estimate 100% error */
165                goto blowup;
166                }
167
168        a0 *= u;
169        sum += a0;
170        t = fabs(a0);
171        if( t > maxt )
172                maxt = t;
173/*
174        if( (maxt/fabs(sum)) > 1.0e17 )
175                {
176                pcanc = 1.0;
177                goto blowup;
178                }
179*/
180        an += 1.0;
181        bn += 1.0;
182        n += 1.0;
183        }
184
185pdone:
186
187/* estimate error due to roundoff and cancellation */
188if( sum != 0.0 )
189        maxt /= fabs(sum);
190maxt *= MACHEP;         /* this way avoids multiply overflow */
191pcanc = fabs( MACHEP * n  +  maxt );
192
193blowup:
194
195*err = pcanc;
196
197return( sum );
198}
199
200
201/*                                                      hy1f1a()        */
202/* asymptotic formula for hypergeometric function:
203 *
204 *        (    -a                         
205 *  --    ( |z|                           
206 * |  (b) ( -------- 2f0( a, 1+a-b, -1/x )
207 *        (  --                           
208 *        ( |  (b-a)                     
209 *
210 *
211 *                                x    a-b                     )
212 *                               e  |x|                        )
213 *                             + -------- 2f0( b-a, 1-a, 1/x ) )
214 *                                --                           )
215 *                               |  (a)                        )
216 */
217
218static double hy1f1a( a, b, x, err )
219double a, b, x;
220double *err;
221{
222double h1, h2, t, u, temp, acanc, asum, err1, err2;
223
224if( x == 0 )
225        {
226        acanc = 1.0;
227        asum = MAXNUM;
228        goto adone;
229        }
230temp = log( fabs(x) );
231t = x + temp * (a-b);
232u = -temp * a;
233
234if( b > 0 )
235        {
236        temp = lgam(b);
237        t += temp;
238        u += temp;
239        }
240
241h1 = hyp2f0( a, a-b+1, -1.0/x, 1, &err1 );
242
243temp = exp(u) / gamma(b-a);
244h1 *= temp;
245err1 *= temp;
246
247h2 = hyp2f0( b-a, 1.0-a, 1.0/x, 2, &err2 );
248
249if( a < 0 )
250        temp = exp(t) / gamma(a);
251else
252        temp = exp( t - lgam(a) );
253
254h2 *= temp;
255err2 *= temp;
256
257if( x < 0.0 )
258        asum = h1;
259else
260        asum = h2;
261
262acanc = fabs(err1) + fabs(err2);
263
264
265if( b < 0 )
266        {
267        temp = gamma(b);
268        asum *= temp;
269        acanc *= fabs(temp);
270        }
271
272
273if( asum != 0.0 )
274        acanc /= fabs(asum);
275
276acanc *= 30.0;  /* fudge factor, since error of asymptotic formula
277                 * often seems this much larger than advertised */
278
279adone:
280
281
282*err = acanc;
283return( asum );
284}
285
286/*                                                      hyp2f0()        */
287
288double hyp2f0( a, b, x, type, err )
289double a, b, x;
290int type;       /* determines what converging factor to use */
291double *err;
292{
293double a0, alast, t, tlast, maxt;
294double n, an, bn, u, sum, temp;
295
296an = a;
297bn = b;
298a0 = 1.0e0;
299alast = 1.0e0;
300sum = 0.0;
301n = 1.0e0;
302t = 1.0e0;
303tlast = 1.0e9;
304maxt = 0.0;
305
306do
307        {
308        if( an == 0 )
309                goto pdone;
310        if( bn == 0 )
311                goto pdone;
312
313        u = an * (bn * x / n);
314
315        /* check for blowup */
316        temp = fabs(u);
317        if( (temp > 1.0 ) && (maxt > (MAXNUM/temp)) )
318                goto error;
319
320        a0 *= u;
321        t = fabs(a0);
322
323        /* terminating condition for asymptotic series */
324        if( t > tlast )
325                goto ndone;
326
327        tlast = t;
328        sum += alast;   /* the sum is one term behind */
329        alast = a0;
330
331        if( n > 200 )
332                goto ndone;
333
334        an += 1.0e0;
335        bn += 1.0e0;
336        n += 1.0e0;
337        if( t > maxt )
338                maxt = t;
339        }
340while( t > MACHEP );
341
342
343pdone:  /* series converged! */
344
345/* estimate error due to roundoff and cancellation */
346*err = fabs(  MACHEP * (n + maxt)  );
347
348alast = a0;
349goto done;
350
351ndone:  /* series did not converge */
352
353/* The following "Converging factors" are supposed to improve accuracy,
354 * but do not actually seem to accomplish very much. */
355
356n -= 1.0;
357x = 1.0/x;
358
359switch( type )  /* "type" given as subroutine argument */
360{
361case 1:
362        alast *= ( 0.5 + (0.125 + 0.25*b - 0.5*a + 0.25*x - 0.25*n)/x );
363        break;
364
365case 2:
366        alast *= 2.0/3.0 - b + 2.0*a + x - n;
367        break;
368
369default:
370        ;
371}
372
373/* estimate error due to roundoff, cancellation, and nonconvergence */
374*err = MACHEP * (n + maxt)  +  fabs ( a0 );
375
376
377done:
378sum += alast;
379return( sum );
380
381/* series blew up: */
382error:
383*err = MAXNUM;
384mtherr( "hyperg", TLOSS );
385return( sum );
386}
Note: See TracBrowser for help on using the repository browser.