← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to scipy/special/cephes/hyperg.c

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

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