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 */ |