Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
cdf_30.c
Go to the documentation of this file.00001 #include "cdflib.h"
00002 void cumfnc(double *f,double *dfn,double *dfd,double *pnonc,
00003 double *cum,double *ccum)
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065 {
00066 #define qsmall(x) (int)(sum < 1.0e-20 || (x) < eps*sum)
00067 #define half 0.5e0
00068 #define done 1.0e0
00069 static double eps = 1.0e-4;
00070 static double dsum,dummy,prod,xx,yy,adn,aup,b,betdn,betup,centwt,dnterm,sum,
00071 upterm,xmult,xnonc;
00072 static int i,icent,ierr;
00073 static double T1,T2,T3,T4,T5,T6;
00074
00075
00076
00077
00078 if(!(*f <= 0.0e0)) goto S10;
00079 *cum = 0.0e0;
00080 *ccum = 1.0e0;
00081 return;
00082 S10:
00083 if(!(*pnonc < 1.0e-10)) goto S20;
00084
00085
00086
00087
00088 cumf(f,dfn,dfd,cum,ccum);
00089 return;
00090 S20:
00091 xnonc = *pnonc/2.0e0;
00092
00093
00094
00095 icent = xnonc;
00096 if(icent == 0) icent = 1;
00097
00098
00099
00100 T1 = (double)(icent+1);
00101 centwt = exp(-xnonc+(double)icent*log(xnonc)-alngam(&T1));
00102
00103
00104
00105
00106
00107 prod = *dfn**f;
00108 dsum = *dfd+prod;
00109 yy = *dfd/dsum;
00110 if(yy > half) {
00111 xx = prod/dsum;
00112 yy = done-xx;
00113 }
00114 else xx = done-yy;
00115 T2 = *dfn*half+(double)icent;
00116 T3 = *dfd*half;
00117 bratio(&T2,&T3,&xx,&yy,&betdn,&dummy,&ierr);
00118 adn = *dfn/2.0e0+(double)icent;
00119 aup = adn;
00120 b = *dfd/2.0e0;
00121 betup = betdn;
00122 sum = centwt*betdn;
00123
00124
00125
00126 xmult = centwt;
00127 i = icent;
00128 T4 = adn+b;
00129 T5 = adn+1.0e0;
00130 dnterm = exp(alngam(&T4)-alngam(&T5)-alngam(&b)+adn*log(xx)+b*log(yy));
00131 S30:
00132 if(qsmall(xmult*betdn) || i <= 0) goto S40;
00133 xmult *= ((double)i/xnonc);
00134 i -= 1;
00135 adn -= 1.0;
00136 dnterm = (adn+1.0)/((adn+b)*xx)*dnterm;
00137 betdn += dnterm;
00138 sum += (xmult*betdn);
00139 goto S30;
00140 S40:
00141 i = icent+1;
00142
00143
00144
00145 xmult = centwt;
00146 if(aup-1.0+b == 0) upterm = exp(-alngam(&aup)-alngam(&b)+(aup-1.0)*log(xx)+
00147 b*log(yy));
00148 else {
00149 T6 = aup-1.0+b;
00150 upterm = exp(alngam(&T6)-alngam(&aup)-alngam(&b)+(aup-1.0)*log(xx)+b*
00151 log(yy));
00152 }
00153 goto S60;
00154 S50:
00155 if(qsmall(xmult*betup)) goto S70;
00156 S60:
00157 xmult *= (xnonc/(double)i);
00158 i += 1;
00159 aup += 1.0;
00160 upterm = (aup+b-2.0e0)*xx/(aup-1.0)*upterm;
00161 betup -= upterm;
00162 sum += (xmult*betup);
00163 goto S50;
00164 S70:
00165 *cum = sum;
00166 *ccum = 0.5e0+(0.5e0-*cum);
00167 return;
00168 #undef qsmall
00169 #undef half
00170 #undef done
00171 }