#include <stdio.h>
#include <math.h>
#include "nrutil.h"
#include "nr.h"
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <time.h>
#include <string.h>
#include <ctype.h>
#include <malloc.h>

#define MAXIT 100
#define EPS 3.0e-7
#define FPMIN 1.0e-30
//#define DIM 1000
#define MAXLENGTH 50000

void fit(float x[], float y[], int ndata, float sig[], int mwt, float *a,
	float *b, float *siga, float *sigb, float *chi2, float *q);
double erfcc(double x);
double betai(double a, double b, double x);
double gammln(double xx);
double betacf(double a, double b, double x);



void usage(){
	printf("Usage: stat_sign1 [-R] [-L] -O originLog(L)File -M mutationLog(L)File [-A alphaFile] [-P prob(normal distribution)File] [-F funcFile]\n");
	printf("-R not use regression mode\n");
	printf("-L use maximal likelihood regression, default: linear regression\n");
	printf("-B use beta distribution to calculate the probability, default use normal distribution\n");
	exit(1);
}

int main(int argc, char **argv)
{
	double *Cn, **Cd, *avgCd, *sigmaCd, alfa, beta,**Cexp,*avgCexp,*sigmaCexp;
	double **deltaC, *AVGdeltaC, *STDdeltaC, *prob1, *prob2, *Z1, *Z2;
	double *alpha, *normalProb;
	int DIM;
	int i,j,k;
	char alnName[50], dataFile[50], alphaFile[50], probFile[50], funcName[50];
	int seqLength;
	int *func;
	double regressionFlag=true;
	bool reciprocalFlag=false;
	bool mlFlag=false;
	bool betaFlag=false;

	if(argc<=1) usage();
	for(i=1;i<argc;++i) {
		char option, optidx;
		if(argv[i][0]=='-') option=argv[i][1];
		else usage();

		switch(option) {
			case 'R':
				regressionFlag=false;
				break;
			case 'B':
				betaFlag=true;
				break;
			case 'L':
				mlFlag=true;
				break;
			case 'r':
				reciprocalFlag=true;
				break;
			case 'O':
				++i;
				optidx=0;
				if(i<argc) strcpy(alnName,argv[i]);
				else usage();
				break;
			case 'M':
				++i;
				optidx=0;
				if(i<argc) strcpy(dataFile,argv[i]);
				else usage();
				break;
			case 'A':
				++i;
                                optidx=0;
                                if(i<argc) strcpy(alphaFile,argv[i]);
                                else usage();
                                break;
			case 'P':
                                ++i;
                                optidx=0;
                                if(i<argc) strcpy(probFile,argv[i]);
                                else usage();
                                break;
			case 'F':
                                ++i;
                                optidx=0;
                                if(i<argc) strcpy(funcName,argv[i]);
                                else usage();
                                break;
			default:
				usage();
		}
	}
/*
	if(argc<3) usage();
	else {
		strcpy(alnName,argv[1]);
		strcpy(dataFile,argv[2]);
		if(argc>3) strcpy(alphaFile,argv[3]);
		if(argc>4) strcpy(probFile,argv[4]);
		if(argc>5) strcpy(funcName,argv[5]);
	}
*/
	//read original log(L) file and get seqLength
	FILE *input;
	char line[MAXLENGTH];
	input=fopen(alnName, "r");
	if(!input) {
		printf("can not open file: %s to read\n",alnName);
		exit(1);
	}
	seqLength=0;
	while(fgets(line,MAXLENGTH,input)!=NULL){
		bool flag=false;
		for(i=0;i<strlen(line);i++) if(isalnum(line[i])) { flag=true; break; }
		if(strlen(line)>1 && flag) {
			++seqLength;
		}
	}
	rewind(input);

	//define data array
	func=ivector(1,seqLength);
	//assign vector and matrix
	Cn=dvector(1,seqLength);
	avgCd=dvector(1,seqLength);
	sigmaCd=dvector(1,seqLength);
	avgCexp=dvector(1,seqLength);
	sigmaCexp=dvector(1,seqLength);
	AVGdeltaC=dvector(1,seqLength);
	STDdeltaC=dvector(1,seqLength);
	prob1=dvector(1,seqLength);
	prob2=dvector(1,seqLength);
	Z1=dvector(1,seqLength);
	Z2=dvector(1,seqLength);
	alpha=dvector(1,seqLength);
	normalProb=dvector(1,seqLength);

	//continue read original log(L) file
	seqLength=0;
	while(fgets(line,MAXLENGTH,input)!=NULL){
		bool flag=false;
		for(i=0;i<strlen(line);i++) if(isalnum(line[i])) { flag=true; break; }
		if(strlen(line)>1 && flag) {
			Cn[++seqLength]=atof(line);
		}
	}
	fclose(input);
	// output seqLength
	printf("After reading original Z-score file, seqlength: %d\n",seqLength);

	//read prob(normal distribution)
	for(i=1;i<=seqLength;i++) normalProb[i]=1.0;
	int tmpRow=0;
	if((input=fopen(probFile,"r"))!=NULL) {
		while(fgets(line,MAXLENGTH,input)!=NULL) {
			bool flag=false;
			for(i=0;i<strlen(line);i++) if(isalnum(line[i])) { flag=true; break; }
			if(strlen(line)>1 && flag) {
				++tmpRow;
				if(tmpRow>seqLength) {printf("row number is wrong in file: %s\n",probFile); exit(1);}
				normalProb[tmpRow]=atof(line);
			}
		}
	}

	//open functional file to compare
	int *cutoff5,*cutoff8,index5=1,index8=1;
	cutoff5=ivector(1,seqLength);
	cutoff8=ivector(1,seqLength);
	for(i=1;i<=seqLength;i++) func[i]=0;
	FILE *cutoffFile;
	char line1[1000];
	char seps[]=" \n";
	char *token;
	if((cutoffFile=fopen(funcName,"r"))!=NULL){
		if(fgets(line1,1000,cutoffFile)!=NULL){
			token=strtok(line1,seps);
			while(token!=NULL){
				if(strcmp(token,"cutoff5:")!=0){
					cutoff5[index5++]=atoi(token);
					func[cutoff5[index5-1]]=1;
				}
				token=strtok(NULL,seps);
			}
		}
		if(fgets(line1,1000,cutoffFile)!=NULL){
			token=strtok(line1,seps);
			while(token!=NULL){
				if(strcmp(token,"cutoff8:")!=0){
					cutoff8[index8++]=atoi(token);
				}
				token=strtok(NULL,seps);
			}
		}
	}
/*	else { printf("can not open file: %s\n",funcName); exit(1); }
	for(i=1;i<=seqLength;i++) func[i]=0;
	for(i=1;i<index5;i++) func[cutoff5[i]]=1;
	printf("cutoff5: ");
	for(i=1;i<index5;i++) printf("%d ",cutoff5[i]);
	printf("\n");
*/

/*	printf("Cn vector\n");
	for(i=1;i<=seqLength;i++) printf("%12.6f\n",Cn[i]);
*/	
	//read alpha
	input=fopen(alphaFile,"r");
	if(!input) { printf("can not open file: %s\n", alphaFile); for(i=1;i<=seqLength;i++) alpha[i]=1.01; }
	else {
		for(i=1;i<=seqLength;i++) {
			fgets(line,MAXLENGTH,input);
			alpha[i]=atof(line);
		}
		fclose(input);
	}
/*
	printf("Alpha vector\n");
	for(i=1;i<=seqLength;i++) printf("%12.6f\n",alpha[i]);
*/
	printf("site origin log<L>    alpha     prob<normal distribution>\n");
	for(i=1;i<=seqLength;i++) printf("%4d %12.6f   %12.6f    %12.6e\n",i,Cn[i],alpha[i],normalProb[i]);

	//read Cd
	input=fopen(dataFile,"r");
	if(!input) { printf("can not open file: %s\n", dataFile); exit(1); }

	//decide the DIM value
	DIM=0;
	while(fgets(line,MAXLENGTH,input)!=NULL){
		bool flag=false;
		for(i=0;i<strlen(line);i++) if(isalnum(line[i])) {flag=true; break;}
		if(strlen(line)>1 && flag) ++DIM;
	}
	rewind(input);

	Cd=dmatrix(1,seqLength,1,DIM);
	Cexp=dmatrix(1,seqLength,1,DIM);
	deltaC=dmatrix(1,seqLength,1,DIM);

	int tmpRowIndex=0, tmpColIndex=0;
	while(fgets(line,MAXLENGTH,input)!=NULL){
		bool flag=false;
		for(i=0;i<strlen(line);i++) if(isalnum(line[i])) {flag=true; break;}
		if(strlen(line)>1 && flag) {
			++tmpRowIndex;
//			printf("tmpRowIndex:%d\n",tmpRowIndex);
			if(tmpRowIndex>DIM) {
				printf("row number is wrong, tmpRowIndex:%d and DIM: %d\n",tmpRowIndex,DIM);
				exit(1);
			}
			tmpColIndex=0;
			token=strtok(line,seps);
			while(token!=NULL) {
				++tmpColIndex;
				if(tmpColIndex>seqLength) {
//					printf("two row have different number of data\n");
//					exit(1);
					break;
				}
				Cd[tmpColIndex][tmpRowIndex]=atof(token);
				token=strtok(NULL,seps);
			}
		}
	}
/*	if(tmpRowIndex!=DIM) {
		printf("the row number is wrong, tmpRowIndex:%d and DIM: %d\n",tmpRowIndex,DIM);
		exit(1);
	}
*/
	fclose(input);
/*
	printf("Cd matrix\n");
	for(i=1;i<=seqLength;i++) {
		for(j=1;j<=DIM;j++) printf("%12.4f ",Cd[i][j]);
		printf("\n");
	}
*/

	//calculate the reciprocal of Cn and Cd matrix
	if(reciprocalFlag) {
/*		for(i=1;i<=seqLength;i++) {
			if(fabs(Cn[i])>1.0e-8) Cn[i]=-1.0/Cn[i];
			for(j=1;j<=DIM;j++) if(fabs(Cd[i][j])>1.0e-8) Cd[i][j]=-1.0/Cd[i][j];
		}
*/
		
	}

	//calculate avgCd and sigmaCd
	for(i=1;i<=seqLength;i++){
		avgCd[i]=0.0;
		for(j=1;j<=DIM;j++) avgCd[i]+=Cd[i][j]/(1.0*DIM);
	}
	for(i=1;i<=seqLength;i++){
		sigmaCd[i]=0.0;
		for(j=1;j<=DIM;j++) sigmaCd[i]+=(Cd[i][j]-avgCd[i])*(Cd[i][j]-avgCd[i])/(DIM-1);
		sigmaCd[i]=sqrt(sigmaCd[i]);
	}

	//reciprocal flag, remove those sites which avgCd/Cn< avg(avgCd/Cn)-sigma(avgCd/Cn)
	if(reciprocalFlag) {
		double *ratio; ratio=dvector(1,seqLength);
		double avgRatio=0.0, sigmaRatio=0.0;
		int ratioLength=0;
		for(i=1;i<=seqLength;i++) {
			if(fabs(Cn[i])>1.0e-8) ratio[i]=avgCd[i]/Cn[i];
			else ratio[i]=0.0;
			printf("i:%4d ratio:%f avgCd:%f Cn:%f\n",i,ratio[i],avgCd[i],Cn[i]);
		}
		for(i=1;i<=seqLength;i++) {
			if(ratio[i]>1.0e-8 && fabs(ratio[i]-1.0)>1.0e-8) {
				avgRatio+=ratio[i];
				ratioLength++;
			}
		}
		avgRatio/=(1.0*ratioLength);
		for(i=1;i<=seqLength;i++) {
			if(ratio[i]>1.0e-8 && fabs(ratio[i]-1.0)>1.0e-8) sigmaRatio+=(ratio[i]-avgRatio)*(ratio[i]-avgRatio);
		}
		sigmaRatio=sqrt(sigmaRatio/(ratioLength-1));
		printf("avgRatio:%f  sigmaRatio:%f ratioLength:%d\n",avgRatio,sigmaRatio,ratioLength);
		for(i=1;i<=seqLength;i++) {
			if(ratio[i]<avgRatio-sigmaRatio) alpha[i]=0.0;
		}
	}

	//find the new seqLength, discard those alpha<0.2 alpha>2.0 and alpha=1.0 site
//	int tmpseqLength=0;
//	for(i=1;i<=seqLength;i++) if(alpha[i]>0.2 && alpha[i]<2.0 && fabs(alpha[i]-1.0)>1.0e-8) tmpseqLength++;

	int *remainSite;
	remainSite=ivector(1,seqLength);
	for(i=1;i<=seqLength;i++) remainSite[i]=i;

	for(i=1;i<=seqLength;i++){
		if(fabs(alpha[i]-1.0)<1.0e-8 || alpha[i]<1.0e-8 || sigmaCd[i]<1.0e-8 || fabs(Cn[i]-avgCd[i])<1.0e-10){
//		if(fabs(alpha[i]-1.0)<1.0e-8 || alpha[i]<1.0e-8 || sigmaCd[i]<1.0e-8 || Cn[i]<avgCd[i]){
//		if(alpha[i]<0.5 || alpha[i]>2.0 || fabs(alpha[i]-1.0)<1.0e-8 || normalProb[i]<0.05) {
			for(j=i;j<=seqLength-1;j++){
				remainSite[j]=remainSite[j+1]; alpha[j]=alpha[j+1]; normalProb[j]=normalProb[j+1];
				Cn[j]=Cn[j+1]; avgCd[j]=avgCd[j+1]; sigmaCd[j]=sigmaCd[j+1];func[j]=func[j+1];
				for(k=1;k<=DIM;k++) Cd[j][k]=Cd[j+1][k];
			}
			seqLength--;
			i--;
		}
	}
/*
	for(i=1;i<=DIM;i++) {
		Cd[1][i]/=seqLength;
		for(j=2;j<=seqLength;j++) Cd[1][i]+=Cd[j][i]/seqLength;
		for(j=2;j<=seqLength;j++) Cd[j][i]=Cd[1][i];
	}
	for(j=1;j<=seqLength;j++) {
		avgCd[j]=0.0;
		for(i=1;i<=DIM;i++) avgCd[j]+=Cd[j][i]/DIM;
		sigmaCd[j]=0.0;
		for(i=1;i<=DIM;i++) sigmaCd[j]+=(Cd[j][i]-avgCd[j])*(Cd[j][i]-avgCd[j]);
		sigmaCd[j]=sqrt(sigmaCd[j]/(DIM-1));
	}
*/
		
	//reciprocal flag
	if(reciprocalFlag) {
		//find the average of Cn and avg(avgCd)
/*		double avgCn=0.0, avg_avgCd=0.0;
		for(i=1;i<=seqLength;i++) {
			avgCn+=Cn[i]/(1.0*seqLength);
			avg_avgCd+=avgCd[i]/(1.0*seqLength);
		}
		if(avgCn<0.0) avgCn*=-1.0; if(avg_avgCd<0.0) avg_avgCd*=-1.0;
		for(i=1;i<=seqLength;i++) {
			Cn[i]/=avgCn;
			avgCd[i]/=avg_avgCd;
			sigmaCd[i]/=avg_avgCd;
			for(j=1;j<=DIM;j++) Cd[i][j]/=avg_avgCd;
		}
*/	}

	double **sort;
	sort=dmatrix(1,seqLength,1,3);

	printf("seqLength: %d\n",seqLength);
	//output site Cn Cd F(Y/N) with Cn(sc)!=-1 or Cn(alpha)!=1
	printf("site      Cn      avgCd    alpha     F\n");
	for(i=1;i<=seqLength;i++){
		printf("%3d %12.4f %12.4f %12.4f %1d\n",remainSite[i],Cn[i],avgCd[i],alpha[i],func[i]);
	}
	printf("\nFunctional site\n");
	for(i=1;i<=seqLength;i++){
		if(func[i]==1) printf("%3d %12.4f %12.4f %12.4f %1d\n",remainSite[i],Cn[i],avgCd[i],alpha[i],func[i]);
	}
/*
	printf("average Cd and sigma Cd\n");
	for(i=1;i<=seqLength;i++) printf("%12.6f %12.6f\n",avgCd[i],sigmaCd[i]);
*/
	//calculate alpha and beta
	float a,b,chi2,q,siga,sigb,*sig,*x,*y;
	int mwt=0;
	x=vector(1,seqLength);
	y=vector(1,seqLength);
	sig=vector(1,seqLength);
	for(i=1;i<=seqLength;i++) {
		if(mlFlag){
			mwt=1;
			x[i]=(float) Cn[i];
			y[i]=(float) avgCd[i];
			sig[i]=sigmaCd[i];
		}
		else {
			x[i]=(float) avgCd[i];
			y[i]=(float) Cn[i];
		}
	}
	fit(x,y,seqLength,sig,mwt,&a,&b,&siga,&sigb,&chi2,&q);
	printf("b:%f a:%f siga: %f sigb: %f chi2: %f q: %f\n",b,a,siga,sigb,chi2,q);
	if(mlFlag) {
		alfa=1.0/b; beta=-b/a;
	}
	else
		alfa=b; beta=a;

	//test
	if(!regressionFlag) {
//		double avg_cn=0.0,avg_cd=0.0;
//		for(i=1;i<=seqLength;i++) { avg_cn+=Cn[i]/seqLength; avg_cd+=avgCd[i]/seqLength; }
//		alfa=avg_cn/avg_cd; beta=0.0;
		alfa=1.0; beta=0.0;
	}

	printf("\nalpha: %12.6f  and beta: %12.6f\n",alfa,beta);

	//calculate Cexp, avgCexp and sigmaCexp
	for(i=1;i<=seqLength;i++){
		avgCexp[i]=0.0;
		for(j=1;j<=DIM;j++) {
			Cexp[i][j]=alfa*Cd[i][j]+beta;
			avgCexp[i]+=Cexp[i][j]/(1.0*DIM);
		}
	}
	for(i=1;i<=seqLength;i++){
		sigmaCexp[i]=0.0;
		for(j=1;j<=DIM;j++) sigmaCexp[i]+=(Cexp[i][j]-avgCexp[i])*(Cexp[i][j]-avgCexp[i])/(DIM-1);
		sigmaCexp[i]=sqrt(sigmaCexp[i]);
	}
	//add experiment test pvalue
	double *pvalue_origin, *pvalue_fit;
	pvalue_origin=dvector(1,seqLength);
	pvalue_fit=dvector(1,seqLength);
	for(i=1;i<=seqLength;i++) {
		pvalue_origin[i]=0.0; pvalue_fit[i]=0.0;
		for(j=1;j<=DIM;j++) {
			if(Cd[i][j]>Cn[i]) pvalue_origin[i]+=1.0;
			if(Cexp[i][j]>Cn[i]) pvalue_fit[i]+=1.0;
		}
		pvalue_origin[i]/=(1.0*DIM);
		pvalue_fit[i]/=(1.0*DIM);
	}
	printf("pvalue from experiment\n");
	printf("site pvalue_origin pvalue_fit\n");
	for(i=1;i<=seqLength;i++) {
		printf("%4d      %8.6f     %8.6f\n", remainSite[i],pvalue_origin[i],pvalue_fit[i]);
	}
/*
	printf("Cexp\n");
	for(i=1;i<=seqLength;i++) {
		for(j=1;j<=DIM;j++) printf("%8.4f ",deltaC[i][j]);
		printf("\n");
	}
	printf("average Cexp and sigma Cexp\n");
	for(i=1;i<=seqLength;i++) printf("%12.6f %12.6f\n",avgCexp[i],sigmaCexp[i]);
*/

	//calculate Z1 and prob1 (Model1)
	if(!betaFlag) {
		for(i=1;i<=seqLength;i++) {
			Z1[i]=(Cn[i]-avgCexp[i])/sigmaCexp[i];
			prob1[i]=0.5*erfcc(Z1[i]/sqrt(2.0));
		}
	}
	else {
		double a1, b1, variance;
		for(i=1;i<=seqLength;i++) {
			variance=sigmaCexp[i]*sigmaCexp[i];
			a1=((avgCexp[i]-avgCexp[i]*avgCexp[i])/variance-1.0)*avgCexp[i];
			b1=a1*(1.0/avgCexp[i]-1.0);
			if(remainSite[i]==15) printf("site 15, its Cn:%10.6f avgCd:%10.6f sigmaCd:%10.6f alpha: %10.6f beta:%10.6f\n",Cn[i],avgCexp[i],sigmaCexp[i],a1,b1);
			prob1[i]=1-betai(a1,b1,Cn[i]);
		}
	}

	//sort prob1 and output
	for(i=1;i<=seqLength;i++) {
		sort[i][1]=i; sort[i][2]=prob1[i]; sort[i][3]=0;
	}
	//for(i=1;i<index5;i++) { for(j=1;j<=seqLength;j++) if(cutoff5[i]==remainSite[j]) sort[cutoff5[i]][3]=1;}
	for(i=1;i<=seqLength;i++) sort[i][3]=func[i];

	for(i=1;i<seqLength;i++) {
		for(j=i+1;j<=seqLength;j++) {
			if(sort[i][2]>sort[j][2]) {
				double sort1, sort2, sort3;
				sort1=sort[i][1]; sort2=sort[i][2]; sort3=sort[i][3];
				sort[i][1]=sort[j][1]; sort[i][2]=sort[j][2]; sort[i][3]=sort[j][3];
				sort[j][1]=sort1; sort[j][2]=sort2; sort[j][3]=sort3;
			}
		}
	}
	printf("\n\nsite    prob1  Likelihood   avg_likelihood avg_exp_logL func:Y/N  alpha   sigmaCd   sigmaCexp ndProb\n");
	for(i=1;i<=seqLength;i++) {
		int site= (int) (sort[i][1]+1.0e-8);
		printf("%4d %12.6e %12.6f %12.6f %12.6f %6.0f %8.5f %10.6f %10.6f %14.6e\n",remainSite[site],sort[i][2], Cn[site], avgCd[site], avgCexp[site],sort[i][3], alpha[site],sigmaCd[site],sigmaCexp[site],normalProb[site]);
		if(i==10) printf("\n");
	}

	//print result to another file
	char tmp_alnName[100];
	strcpy(tmp_alnName, alnName);
	strcat(tmp_alnName, ".rank");
	FILE *tmpOutput;
	tmpOutput=fopen(tmp_alnName,"w");
	for(i=1;i<=seqLength;i++) {
		int site= (int) (sort[i][1]+1.0e-8);
		fprintf(tmpOutput,"%d\n",remainSite[site]);
	}
	fclose(tmpOutput);
	strcpy(tmp_alnName, alnName);
	strcat(tmp_alnName, ".pvalue");
	tmpOutput=fopen(tmp_alnName,"w");
        for(i=1;i<=seqLength;i++) {
        	int site= (int) (sort[i][1]+1.0e-8);
		fprintf(tmpOutput,"%12.6e\n",sort[i][2]);
	}
	fclose(tmpOutput);

	//calculate Z2 and prob2 (Model2)
	for(i=1;i<=seqLength;i++) {
		Z2[i]=(Cn[i]-avgCexp[i])/sigmaCexp[i];
//		Z2[i]=pvalue_origin[i];
	}

	//sort prob2 and output
	for(i=1;i<=seqLength;i++) {
		sort[i][1]=i; sort[i][2]=Z2[i]; sort[i][3]=0;
	}
	for(i=1;i<=seqLength;i++) sort[i][3]=func[i];

	for(i=1;i<seqLength;i++) {
		for(j=i+1;j<=seqLength;j++) {
			if(sort[i][2]<sort[j][2]) {
				double sort1, sort2, sort3;
				sort1=sort[i][1]; sort2=sort[i][2]; sort3=sort[i][3];
				sort[i][1]=sort[j][1]; sort[i][2]=sort[j][2]; sort[i][3]=sort[j][3];
				sort[j][1]=sort1; sort[j][2]=sort2; sort[j][3]=sort3;
			}
		}
	}
	printf("\n\nsite    Z1  Likelihood   avg_likelihood func:Y/N\n");
	for(i=1;i<=seqLength;i++) {
		printf("%4d %12.6e %12.6f %12.6f %6.0f\n",remainSite[(int) sort[i][1]],sort[i][2], Cn[(int) sort[i][1]], avgCd[(int) sort[i][1]], sort[i][3]);
//		if(i==10) printf("\n");
	}


/*	printf("prob1,prob2 and Cn\n");
	for(i=1;i<=seqLength;i++) printf("%12.6e %12.6e %12.6f\n",prob1[i],prob2[i],Cn[i]);
*/
}


double erfcc(double x) {
	double t,z,ans;

	z=fabs(x);
	t=1.0/(1.0+0.5*z);
	ans=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
		t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
		t*(-0.82215223+t*0.17087277)))))))));
	return x >= 0.0 ? ans : 2.0-ans;
}

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 bt;

	if (x < 0.0 || x > 1.0) nrerror("Bad x in routine betai");
	if (x == 0.0 || x == 1.0) bt=0.0;
	else
		bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*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 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)*log(tmp);
	ser=1.000000000190015;
	for (j=0;j<=5;j++) ser += cof[j]/++y;
	return -tmp+log(2.5066282746310005*ser/x);
}


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) nrerror("a or b too big, or MAXIT too small in betacf");
	return h;
}
#undef MAXIT
#undef EPS
#undef FPMIN

