Go to the documentation of this file.00001
00007 #include <fvar.hpp>
00008 #include <math.h>
00009 #define EPS 1.0e-9
00010 #define FPMIN 1.0e-30
00011
00022 double betacf(const double a, const double b, const double x, int MAXIT)
00023 {
00024 int m,m2;
00025 double aa,c,d,del,h,qab,qam,qap;
00026
00027 qab=a+b;
00028 qap=a+1.0;
00029 qam=a-1.0;
00030 c=1.0;
00031 d=1.0-qab*x/qap;
00032 if (fabs(d) < FPMIN) d=FPMIN;
00033 d=1.0/d;
00034 h=d;
00035 for (m=1;m<=MAXIT;m++) {
00036 m2=2*m;
00037 aa=m*(b-m)*x/((qam+m2)*(a+m2));
00038 d=1.0+aa*d;
00039 if (fabs(d) < FPMIN) d=FPMIN;
00040 c=1.0+aa/c;
00041 if (fabs(c) < FPMIN) c=FPMIN;
00042 d=1.0/d;
00043 h *= d*c;
00044 aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
00045 d=1.0+aa*d;
00046 if (fabs(d) < FPMIN) d=FPMIN;
00047 c=1.0+aa/c;
00048 if (fabs(c) < FPMIN) c=FPMIN;
00049 d=1.0/d;
00050 del=d*c;
00051 h *= del;
00052 if (fabs(del-1.0) < EPS) break;
00053 }
00054 if (m > MAXIT) cerr << "a or b too big, or MAXIT too small in betacf"
00055 << endl;
00056 return h;
00057 }
00058 #undef MAXIT
00059 #undef EPS
00060 #undef FPMIN
1.8.0
