← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/special/cephes/igami.c

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

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
}