source: sasview/src/sas/models/c_extension/cephes/incbet.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: 6.8 KB
Line 
1/*                                                      incbet.c
2 *
3 *      Incomplete beta integral
4 *
5 *
6 * SYNOPSIS:
7 *
8 * double a, b, x, y, incbet();
9 *
10 * y = incbet( a, b, x );
11 *
12 *
13 * DESCRIPTION:
14 *
15 * Returns incomplete beta integral of the arguments, evaluated
16 * from zero to x.  The function is defined as
17 *
18 *                  x
19 *     -            -
20 *    | (a+b)      | |  a-1     b-1
21 *  -----------    |   t   (1-t)   dt.
22 *   -     -     | |
23 *  | (a) | (b)   -
24 *                 0
25 *
26 * The domain of definition is 0 <= x <= 1.  In this
27 * implementation a and b are restricted to positive values.
28 * The integral from x to 1 may be obtained by the symmetry
29 * relation
30 *
31 *    1 - incbet( a, b, x )  =  incbet( b, a, 1-x ).
32 *
33 * The integral is evaluated by a continued fraction expansion
34 * or, when b*x is small, by a power series.
35 *
36 * ACCURACY:
37 *
38 * Tested at uniformly distributed random points (a,b,x) with a and b
39 * in "domain" and x between 0 and 1.
40 *                                        Relative error
41 * arithmetic   domain     # trials      peak         rms
42 *    IEEE      0,5         10000       6.9e-15     4.5e-16
43 *    IEEE      0,85       250000       2.2e-13     1.7e-14
44 *    IEEE      0,1000      30000       5.3e-12     6.3e-13
45 *    IEEE      0,10000    250000       9.3e-11     7.1e-12
46 *    IEEE      0,100000    10000       8.7e-10     4.8e-11
47 * Outputs smaller than the IEEE gradual underflow threshold
48 * were excluded from these statistics.
49 *
50 * ERROR MESSAGES:
51 *   message         condition      value returned
52 * incbet domain      x<0, x>1          0.0
53 * incbet underflow                     0.0
54 */
55
56
57/*
58Cephes Math Library, Release 2.8:  June, 2000
59Copyright 1984, 1995, 2000 by Stephen L. Moshier
60*/
61
62#include "mconf.h"
63
64#ifdef DEC
65#define MAXGAM 34.84425627277176174
66#else
67#define MAXGAM 171.624376956302725
68#endif
69
70extern double MACHEP, MINLOG, MAXLOG;
71#ifdef ANSIPROT
72extern double gamma ( double );
73extern double lgam ( double );
74extern double exp ( double );
75extern double log ( double );
76extern double pow ( double, double );
77extern double fabs ( double );
78static double incbcf(double, double, double);
79static double incbd(double, double, double);
80static double pseries(double, double, double);
81#else
82double gamma(), lgam(), exp(), log(), pow(), fabs();
83static double incbcf(), incbd(), pseries();
84#endif
85
86static double big = 4.503599627370496e15;
87static double biginv =  2.22044604925031308085e-16;
88
89
90double incbet( aa, bb, xx )
91double aa, bb, xx;
92{
93double a, b, t, x, xc, w, y;
94int flag;
95
96if( aa <= 0.0 || bb <= 0.0 )
97        goto domerr;
98
99if( (xx <= 0.0) || ( xx >= 1.0) )
100        {
101        if( xx == 0.0 )
102                return(0.0);
103        if( xx == 1.0 )
104                return( 1.0 );
105domerr:
106        mtherr( "incbet", DOMAIN );
107        return( 0.0 );
108        }
109
110flag = 0;
111if( (bb * xx) <= 1.0 && xx <= 0.95)
112        {
113        t = pseries(aa, bb, xx);
114                goto done;
115        }
116
117w = 1.0 - xx;
118
119/* Reverse a and b if x is greater than the mean. */
120if( xx > (aa/(aa+bb)) )
121        {
122        flag = 1;
123        a = bb;
124        b = aa;
125        xc = xx;
126        x = w;
127        }
128else
129        {
130        a = aa;
131        b = bb;
132        xc = w;
133        x = xx;
134        }
135
136if( flag == 1 && (b * x) <= 1.0 && x <= 0.95)
137        {
138        t = pseries(a, b, x);
139        goto done;
140        }
141
142/* Choose expansion for better convergence. */
143y = x * (a+b-2.0) - (a-1.0);
144if( y < 0.0 )
145        w = incbcf( a, b, x );
146else
147        w = incbd( a, b, x ) / xc;
148
149/* Multiply w by the factor
150     a      b   _             _     _
151    x  (1-x)   | (a+b) / ( a | (a) | (b) ) .   */
152
153y = a * log(x);
154t = b * log(xc);
155if( (a+b) < MAXGAM && fabs(y) < MAXLOG && fabs(t) < MAXLOG )
156        {
157        t = pow(xc,b);
158        t *= pow(x,a);
159        t /= a;
160        t *= w;
161        t *= gamma(a+b) / (gamma(a) * gamma(b));
162        goto done;
163        }
164/* Resort to logarithms.  */
165y += t + lgam(a+b) - lgam(a) - lgam(b);
166y += log(w/a);
167if( y < MINLOG )
168        t = 0.0;
169else
170        t = exp(y);
171
172done:
173
174if( flag == 1 )
175        {
176        if( t <= MACHEP )
177                t = 1.0 - MACHEP;
178        else
179                t = 1.0 - t;
180        }
181return( t );
182}
183
184/* Continued fraction expansion #1
185 * for incomplete beta integral
186 */
187
188static double incbcf( a, b, x )
189double a, b, x;
190{
191double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
192double k1, k2, k3, k4, k5, k6, k7, k8;
193double r, t, ans, thresh;
194int n;
195
196k1 = a;
197k2 = a + b;
198k3 = a;
199k4 = a + 1.0;
200k5 = 1.0;
201k6 = b - 1.0;
202k7 = k4;
203k8 = a + 2.0;
204
205pkm2 = 0.0;
206qkm2 = 1.0;
207pkm1 = 1.0;
208qkm1 = 1.0;
209ans = 1.0;
210r = 1.0;
211n = 0;
212thresh = 3.0 * MACHEP;
213do
214        {
215       
216        xk = -( x * k1 * k2 )/( k3 * k4 );
217        pk = pkm1 +  pkm2 * xk;
218        qk = qkm1 +  qkm2 * xk;
219        pkm2 = pkm1;
220        pkm1 = pk;
221        qkm2 = qkm1;
222        qkm1 = qk;
223
224        xk = ( x * k5 * k6 )/( k7 * k8 );
225        pk = pkm1 +  pkm2 * xk;
226        qk = qkm1 +  qkm2 * xk;
227        pkm2 = pkm1;
228        pkm1 = pk;
229        qkm2 = qkm1;
230        qkm1 = qk;
231
232        if( qk != 0 )
233                r = pk/qk;
234        if( r != 0 )
235                {
236                t = fabs( (ans - r)/r );
237                ans = r;
238                }
239        else
240                t = 1.0;
241
242        if( t < thresh )
243                goto cdone;
244
245        k1 += 1.0;
246        k2 += 1.0;
247        k3 += 2.0;
248        k4 += 2.0;
249        k5 += 1.0;
250        k6 -= 1.0;
251        k7 += 2.0;
252        k8 += 2.0;
253
254        if( (fabs(qk) + fabs(pk)) > big )
255                {
256                pkm2 *= biginv;
257                pkm1 *= biginv;
258                qkm2 *= biginv;
259                qkm1 *= biginv;
260                }
261        if( (fabs(qk) < biginv) || (fabs(pk) < biginv) )
262                {
263                pkm2 *= big;
264                pkm1 *= big;
265                qkm2 *= big;
266                qkm1 *= big;
267                }
268        }
269while( ++n < 300 );
270
271cdone:
272return(ans);
273}
274
275
276/* Continued fraction expansion #2
277 * for incomplete beta integral
278 */
279
280static double incbd( a, b, x )
281double a, b, x;
282{
283double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
284double k1, k2, k3, k4, k5, k6, k7, k8;
285double r, t, ans, z, thresh;
286int n;
287
288k1 = a;
289k2 = b - 1.0;
290k3 = a;
291k4 = a + 1.0;
292k5 = 1.0;
293k6 = a + b;
294k7 = a + 1.0;;
295k8 = a + 2.0;
296
297pkm2 = 0.0;
298qkm2 = 1.0;
299pkm1 = 1.0;
300qkm1 = 1.0;
301z = x / (1.0-x);
302ans = 1.0;
303r = 1.0;
304n = 0;
305thresh = 3.0 * MACHEP;
306do
307        {
308       
309        xk = -( z * k1 * k2 )/( k3 * k4 );
310        pk = pkm1 +  pkm2 * xk;
311        qk = qkm1 +  qkm2 * xk;
312        pkm2 = pkm1;
313        pkm1 = pk;
314        qkm2 = qkm1;
315        qkm1 = qk;
316
317        xk = ( z * k5 * k6 )/( k7 * k8 );
318        pk = pkm1 +  pkm2 * xk;
319        qk = qkm1 +  qkm2 * xk;
320        pkm2 = pkm1;
321        pkm1 = pk;
322        qkm2 = qkm1;
323        qkm1 = qk;
324
325        if( qk != 0 )
326                r = pk/qk;
327        if( r != 0 )
328                {
329                t = fabs( (ans - r)/r );
330                ans = r;
331                }
332        else
333                t = 1.0;
334
335        if( t < thresh )
336                goto cdone;
337
338        k1 += 1.0;
339        k2 -= 1.0;
340        k3 += 2.0;
341        k4 += 2.0;
342        k5 += 1.0;
343        k6 += 1.0;
344        k7 += 2.0;
345        k8 += 2.0;
346
347        if( (fabs(qk) + fabs(pk)) > big )
348                {
349                pkm2 *= biginv;
350                pkm1 *= biginv;
351                qkm2 *= biginv;
352                qkm1 *= biginv;
353                }
354        if( (fabs(qk) < biginv) || (fabs(pk) < biginv) )
355                {
356                pkm2 *= big;
357                pkm1 *= big;
358                qkm2 *= big;
359                qkm1 *= big;
360                }
361        }
362while( ++n < 300 );
363cdone:
364return(ans);
365}
366
367/* Power series for incomplete beta integral.
368   Use when b*x is small and x not too close to 1.  */
369
370static double pseries( a, b, x )
371double a, b, x;
372{
373double s, t, u, v, n, t1, z, ai;
374
375ai = 1.0 / a;
376u = (1.0 - b) * x;
377v = u / (a + 1.0);
378t1 = v;
379t = u;
380n = 2.0;
381s = 0.0;
382z = MACHEP * ai;
383while( fabs(v) > z )
384        {
385        u = (n - b) * x / n;
386        t *= u;
387        v = t / (a + n);
388        s += v; 
389        n += 1.0;
390        }
391s += t1;
392s += ai;
393
394u = a * log(x);
395if( (a+b) < MAXGAM && fabs(u) < MAXLOG )
396        {
397        t = gamma(a+b)/(gamma(a)*gamma(b));
398        s = s * t * pow(x,a);
399        }
400else
401        {
402        t = lgam(a+b) - lgam(a) - lgam(b) + u + log(s);
403        if( t < MINLOG )
404                s = 0.0;
405        else
406        s = exp(t);
407        }
408return(s);
409}
Note: See TracBrowser for help on using the repository browser.