Doxygen Source Code Documentation

Main Page Alphabetical List Data Structures File List Data Fields Globals Search
cdf_09.c

00001 #include "cdflib.h"
00002 double bpser(double *a,double *b,double *x,double *eps)
00003 /*
00004 -----------------------------------------------------------------------
00005      POWER SERIES EXPANSION FOR EVALUATING IX(A,B) WHEN B .LE. 1
00006      OR B*X .LE. 0.7.  EPS IS THE TOLERANCE USED.
00007 -----------------------------------------------------------------------
00008 */
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 */
AFNI/NIfTI Server

Sections

Personal tools

Navigation

Doxygen Source Code Documentation

cdf_09.c
