#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define MAXIT 100
#define EPS 3.0e-7
#define FPMIN 1.0e-80

double fprob( double F, int df1, int df2 );
double betai(double a, double b, double x);
double betacf (double a, double b, double x);
double gammln (double xx);
void nrerror(char error_text[]);


double fprob( double F, int df1, int df2 ) {
  /*  fprintf(stderr, "fprob F %lf df1 %d df2 %d\n", F, df1, df2 );  */
  return betai( df2/2.0, df1/2.0, df2/(df2+df1*F));
}

/*
int main (void)
{
	double f;
	double x;
	double df1;
	double df2;
	double prob;
	
printf(" Input F:\n");
	scanf("%f", &f);
printf (" Input df1 \n");
	scanf("%f", &df1);
printf(" Input df2 \n");
	scanf("%f", &df2);

prob = betai(0.5*df2, 0.5*df1, df2/(df2+df1*f));
 if (prob > 1.0 )  prob = 2.0-prob;
printf ("Probability is %1.10f\n", prob );
printf ("Logp is %3.4f\n", -1 * (float)log(prob));
return 0;
}
*/


double betai(double a, double b, double x)

{
double bt;

 if (x < 0.0 || x > 1.0 ) { fprintf(stderr, "Bad x=%lf in routine betai\n", x);exit(1);}
if (x == 0.0 || x == 1.0 ) bt = 0.0;

else
	bt=(double)exp(gammln(a+b) - gammln(a)-gammln(b) + a * (double) log(x) + b * (double)log(1.0 - x));

if(x < (a + 1.0)/(a + b + 2.0))
	return bt*betacf(a,b,x)/a;

else
	return 1.0- bt*betacf(b,a,1.0-x)/b;
	
}

double betacf(double a, double b, double x)
{

void nrerror(char error_text[]);
int m, m2;
double aa, c, d, del, h, qab, qam, qap;

qab = a + b;
qap = a + 1.0;
qam = a - 1.0;
c = 1.0;
 
d = 1.0 - qab * x/qap;
 if (fabs(d) < FPMIN ) d = FPMIN;
 
 d = 1.0/d;
 h = d;
 
 for (m = 1; m <= MAXIT; m++) {
 
 m2 = 2*m;
 aa = m*(b-m)*x/((qam + m2) * (a+m2));
 d = 1.0 + aa*d;
 if (fabs(d) < FPMIN) d = FPMIN;
 c = 1.0 + aa/c;
 if(fabs(c) < FPMIN) c = FPMIN;
 d = 1.0/d;
 h *= d*c;
 aa = -(a+m) * (qab + m)*x/((a + m2) * (qap + m2));
 d = 1.0 + aa * d;
  if (fabs(d) < FPMIN)  d = FPMIN;
  c = 1.0 + aa/c;
  if (fabs(c) < FPMIN)  c = FPMIN;
  d = 1.0/d;
  del = d*c;
   h *= del;
   if (fabs(del-1.0) < EPS) break;
  }
 if (m > MAXIT) {
   fprintf(stderr, "a = %e b = %e x = %e\n", a, b, x );
   nrerror("a or b too big, or MAXIT too small in betacf\n");
 }
 return h;
 }

 
 double gammln (double xx)
 
 { 

	double x,y, tmp, ser;
 
 static double cof[6] = {76.18009172947146, -86.50532032941677, 24.01409824083091,
 -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5};
 int j;
 
 y = x=xx;
 tmp = x + 5.5;
 tmp -= (x + 0.5) *(double) log(tmp);
 ser = 1.000000000190015;
 
 for (j = 0; j <=5; ++j) { 
	ser += cof[j]/++y;	
			}
 
 return -tmp + (double)log(2.5066282746310005 * ser /x);
 
 }
 
 void nrerror(char error_text[])
 { fprintf(stderr, "Numerical Recipes run-time error .... \n");
 fprintf(stderr, "%s\n", error_text);
 fprintf(stderr, "...now exiting to system...\n");
 exit(1);
 }