Doxygen Source Code Documentation
cdf_09.c File Reference
#include "cdflib.h"Go to the source code of this file.
Functions | |
| double | bpser (double *a, double *b, double *x, double *eps) |
Function Documentation
|
||||||||||||||||||||
|
Definition at line 2 of file cdf_09.c. References a, algdiv(), betaln(), bpser(), c, fifdmax1(), fifdmin1(), gam1(), gamln1(), and i. Referenced by bpser(), and bratio().
00009 {
00010 static double bpser,a0,apb,b0,c,n,sum,t,tol,u,w,z;
00011 static int i,m;
00012 /*
00013 ..
00014 .. Executable Statements ..
00015 */
00016 bpser = 0.0e0;
00017 if(*x == 0.0e0) return bpser;
00018 /*
00019 -----------------------------------------------------------------------
00020 COMPUTE THE FACTOR X**A/(A*BETA(A,B))
00021 -----------------------------------------------------------------------
00022 */
00023 a0 = fifdmin1(*a,*b);
00024 if(a0 < 1.0e0) goto S10;
00025 z = *a*log(*x)-betaln(a,b);
00026 bpser = exp(z)/ *a;
00027 goto S100;
00028 S10:
00029 b0 = fifdmax1(*a,*b);
00030 if(b0 >= 8.0e0) goto S90;
00031 if(b0 > 1.0e0) goto S40;
00032 /*
00033 PROCEDURE FOR A0 .LT. 1 AND B0 .LE. 1
00034 */
00035 bpser = pow(*x,*a);
00036 if(bpser == 0.0e0) return bpser;
00037 apb = *a+*b;
00038 if(apb > 1.0e0) goto S20;
00039 z = 1.0e0+gam1(&apb);
00040 goto S30;
00041 S20:
00042 u = *a+*b-1.e0;
00043 z = (1.0e0+gam1(&u))/apb;
00044 S30:
00045 c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z;
00046 bpser *= (c*(*b/apb));
00047 goto S100;
00048 S40:
00049 /*
00050 PROCEDURE FOR A0 .LT. 1 AND 1 .LT. B0 .LT. 8
00051 */
00052 u = gamln1(&a0);
00053 m = b0-1.0e0;
00054 if(m < 1) goto S60;
00055 c = 1.0e0;
00056 for(i=1; i<=m; i++) {
00057 b0 -= 1.0e0;
00058 c *= (b0/(a0+b0));
00059 }
00060 u = log(c)+u;
00061 S60:
00062 z = *a*log(*x)-u;
00063 b0 -= 1.0e0;
00064 apb = a0+b0;
00065 if(apb > 1.0e0) goto S70;
00066 t = 1.0e0+gam1(&apb);
00067 goto S80;
00068 S70:
00069 u = a0+b0-1.e0;
00070 t = (1.0e0+gam1(&u))/apb;
00071 S80:
00072 bpser = exp(z)*(a0/ *a)*(1.0e0+gam1(&b0))/t;
00073 goto S100;
00074 S90:
00075 /*
00076 PROCEDURE FOR A0 .LT. 1 AND B0 .GE. 8
00077 */
00078 u = gamln1(&a0)+algdiv(&a0,&b0);
00079 z = *a*log(*x)-u;
00080 bpser = a0/ *a*exp(z);
00081 S100:
00082 if(bpser == 0.0e0 || *a <= 0.1e0**eps) return bpser;
00083 /*
00084 -----------------------------------------------------------------------
00085 COMPUTE THE SERIES
00086 -----------------------------------------------------------------------
00087 */
00088 sum = n = 0.0e0;
00089 c = 1.0e0;
00090 tol = *eps/ *a;
00091 S110:
00092 n += 1.0e0;
00093 c *= ((0.5e0+(0.5e0-*b/n))**x);
00094 w = c/(*a+n);
00095 sum += w;
00096 if(fabs(w) > tol) goto S110;
00097 bpser *= (1.0e0+*a*sum);
00098 return bpser;
00099 } /* END */
|
