Revision 692 branches/mergetrunkdavef/src/df1b2separable/df1b2bet.cpp
df1b2bet.cpp (revision 692)  

9  9 
#define EPS double(1.0e9) 
10  10 
#define FPMIN double(1.0e30) 
11  11 
df1b2variable betacf(const df1b2variable& a,const df1b2variable& b, 
12 
double x, int MAXIT); 

12 
double x, int MAXIT=100);


13  13  
14  14  
15  15 
//df1b2variable betai(const df1b2variable& a,const df1b2variable& b, 
16  16 
// double x, int maxit=100); 
17 
df3_two_variable betacf(df3_two_variable& a, df3_two_variable& b, 

18 
double x,int MAXIT=100); 

17  19  
18  20 
/** Incomplete beta function for df1b2variable objects. 
19  21 
\param a \f$a\f$ 
...  ...  
41  43 
return 1.0bt*betacf(b,a,1.0x,maxit)/b; 
42  44 
} 
43  45  
46  
47 
df3_two_variable betacf(df3_two_variable& a, df3_two_variable& b, 

48 
double x,int MAXIT) 

49 
{ 

50 
int m,m2; 

51 
df3_two_variable aa,c,d,del,h,qab,qam,qap; 

52 


53 
qab=a+b; 

54 
qap=a+1.0; 

55 
qam=a1.0; 

56 
c=1.0; 

57 
d=1.0qab*x/qap; 

58 
if (fabs(value(d)) < FPMIN) d=FPMIN; 

59 
d=1.0/d; 

60 
h=d; 

61 
for (m=1;m<=MAXIT;m++) 

62 
{ 

63 
m2=2*m; 

64 
aa=m*(bm)*x/((qam+m2)*(a+m2)); 

65 
d=1.0+aa*d; 

66 
if (fabs(value(d)) < FPMIN) d=FPMIN; 

67 
c=1.0+aa/c; 

68 
if (fabs(value(c)) < FPMIN) c=FPMIN; 

69 
d=1.0/d; 

70 
h *= d*c; 

71 
aa = (a+m)*(qab+m)*x/((a+m2)*(qap+m2)); 

72 
d=1.0+aa*d; 

73 
if (fabs(value(d)) < FPMIN) d=FPMIN; 

74 
c=1.0+aa/c; 

75 
if (fabs(value(c)) < FPMIN) c=FPMIN; 

76 
d=1.0/d; 

77 


78 
del=d*c; 

79 
h *= del; 

80 
if (fabs(value(del)1.0) < EPS) break; 

81 
} 

82 
if (m > MAXIT) 

83 
{ 

84 
cerr << "a or b too big, or MAXIT too small in betacf" 

85 
<< endl; 

86 
ad_exit(1); 

87 
} 

88 
return h; 

89 
} 

90 


91  
44  92 
/** Incomplete beta function for df1b2variable objects. 
45  93 
Evaluates the continued fraction for imcomplete beta function. 
46  94 
\param _a \f$a\f$ 
...  ...  
53  101 
"Numerical Recipes in C", 2nd edition, 
54  102 
Press, Teukolsky, Vetterling, Flannery, chapter 2 
55  103 
*/ 
56 
df1b2variable betacf(const df1b2variable& a,const df1b2variable& b,


57 
double x, int MAXIT)


104 
df1b2variable betacf(const df1b2variable& _xa,const df1b2variable& _xb,


105 
double x,int MAXIT)


58  106 
{ 
59 
int m,m2; 

60 
df1b2variable aa,c,d,del,h,qab,qam,qap; 

107 
ADUNCONST(df1b2variable,xa) 

108 
ADUNCONST(df1b2variable,xb) 

109 
init_df3_two_variable a(xa); 

110 
init_df3_two_variable b(xb); 

111 
df3_two_variable z=betacf(a,b,x,MAXIT); 

112 
df1b2variable tmp; 

113 
tmp=z; 

114 
return tmp; 

115 
} 

61  116  
62 
qab=a+b; 

63 
qap=a+double(1.0); 

64 
qam=adouble(1.0); 

65 
c=double(1.0); 

66 
d=double(1.0)qab*x/qap; 

67 
if (fabs(value(d)) < FPMIN) d=FPMIN; 

68 
d=double(1.0)/d; 

69 
h=d; 
Also available in: Unified diff