source code for betai computation

[mvb]

Hello, anyone can suggest me where I can find some good C++ code for the
computaion of betai function?

Thank you in advance,

Luigi


P.S. Any chance to have it in C++0x?

 

[Victor Bazarov]

Hello, anyone can suggest me where I can find some good C++ code for
the computaion of betai function?

Have you tried googling? Have you tried asking in a math NG?

P.S. Any chance to have it in C++0x?

There are several optional library extensions planned for it, which
will contain some common mathematical functions, like gamma or what
have you. Check the Standard Committee's page for the latest draft
document and see for yourself.

V

 

[mvb]

Thank you very much, Victor


Luigi

 

[Lionel B]

Hello, anyone can suggest me where I can find some good C++ code for the
computaion of betai function?

I think the GSL (GNU Scientific Library) has code for that. It's C rather than
C++, strictly speaking, but C++-friendly, so to speak:

http://www.gnu.org/software/gsl/

http://www.gnu.org/software/gsl/manual/html_node/Incomplete-Beta-Function.html

(is that the one you're after?)

The source code is available (it's GNU hence GPL licence).

 

[Kai-Uwe Bux]

Hello, anyone can suggest me where I can find some good C++ code for the
computaion of betai function?

The following is a straight forward implementation of Algorithm 179 (Comm.
ACM Vol 6, Nr 6; 1963). It follows the published algorithm as closely as
possible in C++. That might imply that the code is not "good", but at least
it can be trusted to some degree. The only change of significance is that
it uses gamma() instead of factorial().


typedef long double real;

real regularized_beta ( real x, real p, real q,
real const precision = 0.000000001L ) {
// Algorithm 179
assert( 0.0L <= x && x <= 1.0L );
assert( 0.0L < p );
assert( 0.0L < q );

if ( x == 0.0L ) {
return ( x );
}
if ( x == 1.0L ) {
return ( x );
}

if ( x < 0.5L ) {
return ( 1.0L - regularized_beta( 1.0L - x, q, p ) );
}

real finsum = 0.0L;
real term = 1.0L;
real temp = 1.0L - x;
real qrecur = q;
real index = q;
for ( index -= 1.0L; index > 0; index -= 1.0L ) {
qrecur = index;
term = term * ( qrecur+1.0L) / ( temp * ( p + qrecur ) );
finsum += term;
}
real infsum = 1.0L;
term = 1.0L;
index = 0.0L;
for ( index += 1.0; ( term / infsum ) > epsilon; index += 1.0 ) {
term = term * x * ( index - qrecur ) * ( p + index - 1.0L )
/ ( index * ( p + index ) );
infsum += term;
}
temp = gamma( qrecur );
real temp1 = temp;
term = gamma( qrecur + p );
real term1 = term;
for ( index = qrecur; index + 0.5L < q; index += 1.0L ) {
temp1 = temp1 * index;
term1 = term1 * ( index + p );
}
temp = pow(x,p) *
( infsum * term / (p*temp)
+
finsum * term1 * pow(1-x,q) / ( q *temp1 )) / gamma(p);
return ( temp );
}


Best

Kai-Uwe Bux