← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/special/cephes/igami.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:
Lib/special/cephes/igami.c
1
 
/*                                                      igami()
2
 
 *
3
 
 *      Inverse of complemented imcomplete Gamma integral
4
 
 *
5
 
 *
6
 
 *
7
 
 * SYNOPSIS:
8
 
 *
9
 
 * double a, x, p, igami();
10
 
 *
11
 
 * x = igami( a, p );
12
 
 *
13
 
 * DESCRIPTION:
14
 
 *
15
 
 * Given p, the function finds x such that
16
 
 *
17
 
 *  igamc( a, x ) = p.
18
 
 *
19
 
 * Starting with the approximate value
20
 
 *
21
 
 *         3
22
 
 *  x = a t
23
 
 *
24
 
 *  where
25
 
 *
26
 
 *  t = 1 - d - ndtri(p) sqrt(d)
27
 
 * 
28
 
 * and
29
 
 *
30
 
 *  d = 1/9a,
31
 
 *
32
 
 * the routine performs up to 10 Newton iterations to find the
33
 
 * root of igamc(a,x) - p = 0.
34
 
 *
35
 
 * ACCURACY:
36
 
 *
37
 
 * Tested at random a, p in the intervals indicated.
38
 
 *
39
 
 *                a        p                      Relative error:
40
 
 * arithmetic   domain   domain     # trials      peak         rms
41
 
 *    IEEE     0.5,100   0,0.5       100000       1.0e-14     1.7e-15
42
 
 *    IEEE     0.01,0.5  0,0.5       100000       9.0e-14     3.4e-15
43
 
 *    IEEE    0.5,10000  0,0.5        20000       2.3e-13     3.8e-14
44
 
 */
45
 
 
46
 
/*
47
 
Cephes Math Library Release 2.3:  March, 1995
48
 
Copyright 1984, 1987, 1995 by Stephen L. Moshier
49
 
*/
50
 
 
51
 
#include "mconf.h"
52
 
#include <stdio.h>
53
 
 
54
 
extern double MACHEP, MAXNUM, MAXLOG, MINLOG, NAN;
55
 
#ifndef ANSIPROT
56
 
double igamc(), ndtri(), exp(), fabs(), log(), sqrt(), lgam();
57
 
#endif
58
 
 
59
 
double igami( a, y0 )
60
 
double a, y0;
61
 
{
62
 
double x0, x1, x, yl, yh, y, d, lgm, dithresh;
63
 
int i, dir;
64
 
 
65
 
/* bound the solution */
66
 
x0 = MAXNUM;
67
 
yl = 0;
68
 
x1 = 0;
69
 
yh = 1.0;
70
 
dithresh = 5.0 * MACHEP;
71
 
 
72
 
if ((y0<0.0) || (y0>1.0) || (a<=0)) {
73
 
   mtherr("igami", DOMAIN);
74
 
   return(NAN);
75
 
}
76
 
 
77
 
if (y0==0.0) {
78
 
  return(MAXNUM);
79
 
}
80
 
 
81
 
if (y0==1.0){
82
 
   return 0.0;
83
 
}
84
 
 
85
 
/* approximation to inverse function */
86
 
d = 1.0/(9.0*a);
87
 
y = ( 1.0 - d - ndtri(y0) * sqrt(d) );
88
 
x = a * y * y * y;
89
 
 
90
 
lgm = lgam(a);
91
 
 
92
 
for( i=0; i<10; i++ )
93
 
        {
94
 
        if( x > x0 || x < x1 )
95
 
                goto ihalve;
96
 
        y = igamc(a,x);
97
 
        if( y < yl || y > yh )
98
 
                goto ihalve;
99
 
        if( y < y0 )
100
 
                {
101
 
                x0 = x;
102
 
                yl = y;
103
 
                }
104
 
        else
105
 
                {
106
 
                x1 = x;
107
 
                yh = y;
108
 
                }
109
 
/* compute the derivative of the function at this point */
110
 
        d = (a - 1.0) * log(x) - x - lgm;
111
 
        if( d < -MAXLOG )
112
 
                goto ihalve;
113
 
        d = -exp(d);
114
 
/* compute the step to the next approximation of x */
115
 
        d = (y - y0)/d;
116
 
        if( fabs(d/x) < MACHEP )
117
 
                goto done;
118
 
        x = x - d;
119
 
        }
120
 
 
121
 
/* Resort to interval halving if Newton iteration did not converge. */
122
 
ihalve:
123
 
 
124
 
d = 0.0625;
125
 
if( x0 == MAXNUM )
126
 
        {
127
 
        if( x <= 0.0 )
128
 
                x = 1.0;
129
 
        while( x0 == MAXNUM )
130
 
                {
131
 
                x = (1.0 + d) * x;
132
 
                y = igamc( a, x );
133
 
                if( y < y0 )
134
 
                        {
135
 
                        x0 = x;
136
 
                        yl = y;
137
 
                        break;
138
 
                        }
139
 
                d = d + d;
140
 
                }
141
 
        }
142
 
d = 0.5;
143
 
dir = 0;
144
 
 
145
 
for( i=0; i<400; i++ )
146
 
        {
147
 
        x = x1  +  d * (x0 - x1);
148
 
        y = igamc( a, x );
149
 
        lgm = (x0 - x1)/(x1 + x0);
150
 
        if( fabs(lgm) < dithresh )
151
 
                break;
152
 
        lgm = (y - y0)/y0;
153
 
        if( fabs(lgm) < dithresh )
154
 
                break;
155
 
        if( x <= 0.0 )
156
 
                break;
157
 
        if( y >= y0 )
158
 
                {
159
 
                x1 = x;
160
 
                yh = y;
161
 
                if( dir < 0 )
162
 
                        {
163
 
                        dir = 0;
164
 
                        d = 0.5;
165
 
                        }
166
 
                else if( dir > 1 )
167
 
                        d = 0.5 * d + 0.5; 
168
 
                else
169
 
                        d = (y0 - yl)/(yh - yl);
170
 
                dir += 1;
171
 
                }
172
 
        else
173
 
                {
174
 
                x0 = x;
175
 
                yl = y;
176
 
                if( dir > 0 )
177
 
                        {
178
 
                        dir = 0;
179
 
                        d = 0.5;
180
 
                        }
181
 
                else if( dir < -1 )
182
 
                        d = 0.5 * d;
183
 
                else
184
 
                        d = (y0 - yl)/(yh - yl);
185
 
                dir -= 1;
186
 
                }
187
 
        }
188
 
if( x == 0.0 )
189
 
        mtherr( "igami", UNDERFLOW );
190
 
 
191
 
done:
192
 
return( x );
193
 
}