Revision 692 branches/merge-trunk-davef/src/df1b2-separable/df1b2bet.cpp

df1b2bet.cpp (revision 692)
9 9
#define EPS double(1.0e-9)
10 10
#define FPMIN double(1.0e-30)
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.0-bt*betacf(b,a,1.0-x,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=a-1.0;
56
    c=1.0;
57
    d=1.0-qab*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*(b-m)*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=a-double(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;
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff