source: sasview/src/sans/models/c_extension/cephes/incbi.c @ 79d5b6c

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 79d5b6c was 230f479, checked in by Mathieu Doucet <doucetm@…>, 11 years ago

Rename C source dir for models (minor updates)

  • Property mode set to 100644
File size: 5.0 KB
Line 
1/*                                                      incbi()
2 *
3 *      Inverse of imcomplete beta integral
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * double a, b, x, y, incbi();
10 *
11 * x = incbi( a, b, y );
12 *
13 *
14 *
15 * DESCRIPTION:
16 *
17 * Given y, the function finds x such that
18 *
19 *  incbet( a, b, x ) = y .
20 *
21 * The routine performs interval halving or Newton iterations to find the
22 * root of incbet(a,b,x) - y = 0.
23 *
24 *
25 * ACCURACY:
26 *
27 *                      Relative error:
28 *                x     a,b
29 * arithmetic   domain  domain  # trials    peak       rms
30 *    IEEE      0,1    .5,10000   50000    5.8e-12   1.3e-13
31 *    IEEE      0,1   .25,100    100000    1.8e-13   3.9e-15
32 *    IEEE      0,1     0,5       50000    1.1e-12   5.5e-15
33 *    VAX       0,1    .5,100     25000    3.5e-14   1.1e-15
34 * With a and b constrained to half-integer or integer values:
35 *    IEEE      0,1    .5,10000   50000    5.8e-12   1.1e-13
36 *    IEEE      0,1    .5,100    100000    1.7e-14   7.9e-16
37 * With a = .5, b constrained to half-integer or integer values:
38 *    IEEE      0,1    .5,10000   10000    8.3e-11   1.0e-11
39 */
40
41
42/*
43Cephes Math Library Release 2.8:  June, 2000
44Copyright 1984, 1996, 2000 by Stephen L. Moshier
45*/
46
47#include "mconf.h"
48
49extern double MACHEP, MAXNUM, MAXLOG, MINLOG;
50#ifdef ANSIPROT
51extern double ndtri ( double );
52extern double exp ( double );
53extern double fabs ( double );
54extern double log ( double );
55extern double sqrt ( double );
56extern double lgam ( double );
57extern double incbet ( double, double, double );
58#else
59double ndtri(), exp(), fabs(), log(), sqrt(), lgam(), incbet();
60#endif
61
62double incbi( aa, bb, yy0 )
63double aa, bb, yy0;
64{
65double a, b, y0, d, y, x, x0, x1, lgm, yp, di, dithresh, yl, yh, xt;
66int i, rflg, dir, nflg;
67
68
69i = 0;
70if( yy0 <= 0 )
71        return(0.0);
72if( yy0 >= 1.0 )
73        return(1.0);
74x0 = 0.0;
75yl = 0.0;
76x1 = 1.0;
77yh = 1.0;
78nflg = 0;
79
80if( aa <= 1.0 || bb <= 1.0 )
81        {
82        dithresh = 1.0e-6;
83        rflg = 0;
84        a = aa;
85        b = bb;
86        y0 = yy0;
87        x = a/(a+b);
88        y = incbet( a, b, x );
89        goto ihalve;
90        }
91else
92        {
93        dithresh = 1.0e-4;
94        }
95/* approximation to inverse function */
96
97yp = -ndtri(yy0);
98
99if( yy0 > 0.5 )
100        {
101        rflg = 1;
102        a = bb;
103        b = aa;
104        y0 = 1.0 - yy0;
105        yp = -yp;
106        }
107else
108        {
109        rflg = 0;
110        a = aa;
111        b = bb;
112        y0 = yy0;
113        }
114
115lgm = (yp * yp - 3.0)/6.0;
116x = 2.0/( 1.0/(2.0*a-1.0)  +  1.0/(2.0*b-1.0) );
117d = yp * sqrt( x + lgm ) / x
118        - ( 1.0/(2.0*b-1.0) - 1.0/(2.0*a-1.0) )
119        * (lgm + 5.0/6.0 - 2.0/(3.0*x));
120d = 2.0 * d;
121if( d < MINLOG )
122        {
123        x = 1.0;
124        goto under;
125        }
126x = a/( a + b * exp(d) );
127y = incbet( a, b, x );
128yp = (y - y0)/y0;
129if( fabs(yp) < 0.2 )
130        goto newt;
131
132/* Resort to interval halving if not close enough. */
133ihalve:
134
135dir = 0;
136di = 0.5;
137for( i=0; i<100; i++ )
138        {
139        if( i != 0 )
140                {
141                x = x0  +  di * (x1 - x0);
142                if( x == 1.0 )
143                        x = 1.0 - MACHEP;
144                if( x == 0.0 )
145                        {
146                        di = 0.5;
147                        x = x0  +  di * (x1 - x0);
148                        if( x == 0.0 )
149                                goto under;
150                        }
151                y = incbet( a, b, x );
152                yp = (x1 - x0)/(x1 + x0);
153                if( fabs(yp) < dithresh )
154                        goto newt;
155                yp = (y-y0)/y0;
156                if( fabs(yp) < dithresh )
157                        goto newt;
158                }
159        if( y < y0 )
160                {
161                x0 = x;
162                yl = y;
163                if( dir < 0 )
164                        {
165                        dir = 0;
166                        di = 0.5;
167                        }
168                else if( dir > 3 )
169                        di = 1.0 - (1.0 - di) * (1.0 - di);
170                else if( dir > 1 )
171                        di = 0.5 * di + 0.5; 
172                else
173                        di = (y0 - y)/(yh - yl);
174                dir += 1;
175                if( x0 > 0.75 )
176                        {
177                        if( rflg == 1 )
178                                {
179                                rflg = 0;
180                                a = aa;
181                                b = bb;
182                                y0 = yy0;
183                                }
184                        else
185                                {
186                                rflg = 1;
187                                a = bb;
188                                b = aa;
189                                y0 = 1.0 - yy0;
190                                }
191                        x = 1.0 - x;
192                        y = incbet( a, b, x );
193                        x0 = 0.0;
194                        yl = 0.0;
195                        x1 = 1.0;
196                        yh = 1.0;
197                        goto ihalve;
198                        }
199                }
200        else
201                {
202                x1 = x;
203                if( rflg == 1 && x1 < MACHEP )
204                        {
205                        x = 0.0;
206                        goto done;
207                        }
208                yh = y;
209                if( dir > 0 )
210                        {
211                        dir = 0;
212                        di = 0.5;
213                        }
214                else if( dir < -3 )
215                        di = di * di;
216                else if( dir < -1 )
217                        di = 0.5 * di;
218                else
219                        di = (y - y0)/(yh - yl);
220                dir -= 1;
221                }
222        }
223mtherr( "incbi", PLOSS );
224if( x0 >= 1.0 )
225        {
226        x = 1.0 - MACHEP;
227        goto done;
228        }
229if( x <= 0.0 )
230        {
231under:
232        mtherr( "incbi", UNDERFLOW );
233        x = 0.0;
234        goto done;
235        }
236
237newt:
238
239if( nflg )
240        goto done;
241nflg = 1;
242lgm = lgam(a+b) - lgam(a) - lgam(b);
243
244for( i=0; i<8; i++ )
245        {
246        /* Compute the function at this point. */
247        if( i != 0 )
248                y = incbet(a,b,x);
249        if( y < yl )
250                {
251                x = x0;
252                y = yl;
253                }
254        else if( y > yh )
255                {
256                x = x1;
257                y = yh;
258                }
259        else if( y < y0 )
260                {
261                x0 = x;
262                yl = y;
263                }
264        else
265                {
266                x1 = x;
267                yh = y;
268                }
269        if( x == 1.0 || x == 0.0 )
270                break;
271        /* Compute the derivative of the function at this point. */
272        d = (a - 1.0) * log(x) + (b - 1.0) * log(1.0-x) + lgm;
273        if( d < MINLOG )
274                goto done;
275        if( d > MAXLOG )
276                break;
277        d = exp(d);
278        /* Compute the step to the next approximation of x. */
279        d = (y - y0)/d;
280        xt = x - d;
281        if( xt <= x0 )
282                {
283                y = (x - x0) / (x1 - x0);
284                xt = x0 + 0.5 * y * (x - x0);
285                if( xt <= 0.0 )
286                        break;
287                }
288        if( xt >= x1 )
289                {
290                y = (x1 - x) / (x1 - x0);
291                xt = x1 - 0.5 * y * (x1 - x);
292                if( xt >= 1.0 )
293                        break;
294                }
295        x = xt;
296        if( fabs(d/x) < 128.0 * MACHEP )
297                goto done;
298        }
299/* Did not converge.  */
300dithresh = 256.0 * MACHEP;
301goto ihalve;
302
303done:
304
305if( rflg )
306        {
307        if( x <= MACHEP )
308                x = 1.0 - MACHEP;
309        else
310                x = 1.0 - x;
311        }
312return( x );
313}
Note: See TracBrowser for help on using the repository browser.