public class sayisal
{
public static double gammln(double xx)
{
	//gamma function
	double x,tmp,ser;
	double cof[]={76.18009173,-86.50532033,24.01409822,
		-1.231739516,0.120858003e-2,-0.536382e-5};
	int j;

	x=xx-1.0;
	tmp=x+5.5;
	tmp -= (x+0.5)*Math.log(tmp);
	ser=1.0;
	for (j=0;j<=5;j++) {
		x += 1.0;
		ser += cof[j]/x;
	}
	return -tmp+Math.log(2.50662827465*ser);
}

public static double factrl(int n)
{
	// n! of an integer number
	int ntop=4;
	double a[]=new double[33];
	a[0]=1.0;
	a[1]=1.0;
	a[2]=2.0;
	a[3]=6.0;
	int j;
	if (n < 0) System.out.println("Negative factorial in routine FACTRL");
	if (n > 32) return Math.exp(gammln(n+1.0));
	while (ntop<=n) {
		j=ntop-1;
		a[ntop]=a[j]*ntop;
		ntop++;
	}
	return a[n];
}

public static double factln(int n)
{
	//Math.log(factorial) of an integer number
	double a[]=new double[33];
	int ntop=4;
	a[0]=1.0;
	a[1]=1.0;
	a[2]=2.0;
	a[3]=6.0;
	int j;
	if (n < 0) System.out.println("Negative factorial in routine FACTLN");
	if (n <= 1) return 0.0;
	if (n > 32) return gammln(n+1.0);
	while (ntop<=n) {
		j=ntop-1;
		a[ntop]=a[j]*ntop;
		ntop++;
	}
	return Math.log(a[n]);
}

public static double  bico(int n,int k)
{
	//binomial coefficients
	return Math.floor(0.5+Math.exp(factln(n)-factln(k)-factln(n-k)));
}

public static double beta(double z,double w)
{
	//gamma(z)*gamma(w)/gamma(z+w);
	return Math.exp(gammln(z)+gammln(w)-gammln(z+w));
}

public static double gammp(double a,double x)
{
	//incomplete gamma function P(a,x)
	double gamser,gammcf;
	if (x < 0.0 || a <= 0.0) System.out.println("Invalid arguments in routine GAMMP");
	if (x < (a+1.0)) {
		gamser=gser(a,x);
		return gamser;
	} else {
		gammcf=gcf(a,x);
		return 1.0-gammcf;
	}
}

public static double gammq(double a,double x)
{
	//incomplete gamma function Q(a,x)=1-P(a,x);
	double gamser,gammcf;
	if (x < 0.0 || a <= 0.0) System.out.println("Invalid arguments in routine GAMMQ");
	if (x < (a+1.0)) {
		gamser=gser(a,x);
		return 1.0-gamser;
	} else {
		gammcf=gcf(a,x);
		return gammcf;
	}
}

public static double gser(double a,double x)
{
	int ITMAX=100;
	double EPS=3.0e-7;
	int n;
	double sum,del,ap;
	double gln=gammln(a);
	double gamser=0;
	if (x <= 0.0) {
		if (x < 0.0) System.out.println("x less than 0 in routine GSER");
		gamser=0.0;
		return gamser;
	} else {
		ap=a;
		del=sum=1.0/a;
		for (n=1;n<=ITMAX;n++) {
			ap += 1.0;
			del *= x/ap;
			sum += del;
			if (Math.abs(del) < Math.abs(sum)*EPS) {
				gamser=sum*Math.exp(-x+a*Math.log(x)-(gln));
				return gamser;
			}
		}
		System.out.println("a too large, ITMAX too small in routine GSER");
		return gamser;
	}
}

public static double gcf(double a,double x)
{
	int ITMAX=100;
	double EPS=3.0e-7;
	int n;
	double gln;
	double gammcf=0;
	double gold=0.0,g,fac=1.0,b1=1.0;
	double b0=0.0,anf,ana,an,a1,a0=1.0;
	gln=gammln(a);
	a1=x;
	for (n=1;n<=ITMAX;n++) {
		an=(double) n;
		ana=an-a;
		a0=(a1+a0*ana)*fac;
		b0=(b1+b0*ana)*fac;
		anf=an*fac;
		a1=x*a0+anf*a1;
		b1=x*b0+anf*b1;
		if (a1!=0) {
			fac=1.0/a1;
			g=b1*fac;
			if (Math.abs((g-gold)/g) < EPS) {
				gammcf=Math.exp(-x+a*Math.log(x)-gln)*g;
				return gammcf;
			}
			gold=g;
		}
	}
	System.out.println("a too large, ITMAX too small in routine GCF");
    return gammcf;
}

public static double  erf(double x)
{
   //error function
   return x < 0.0 ? -gammp(0.5,x*x) : gammp(0.5,x*x);
}

public static double  erfc(double x)
{
	//complimentary error function
	return x < 0.0 ? 1.0+gammp(0.5,x*x) : gammq(0.5,x*x);
}

public static double  erfcc(double x)
{
	double t,z,ans;
	z=Math.abs(x);
	t=1.0/(1.0+0.5*z);
	ans=t*Math.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 betacf(double a,double b,double x)
{
	int ITMAX=100;
	double EPS=3.0e-7;
	double qap,qam,qab,em,tem,d;
	double bz,bm=1.0,bp,bpp;
	double az=1.0,am=1.0,ap,app,aold;
	int m;
	qab=a+b;
	qap=a+1.0;
	qam=a-1.0;
	bz=1.0-qab*x/qap;
	for (m=1;m<=ITMAX;m++) {
		em=(double) m;
		tem=em+em;
		d=em*(b-em)*x/((qam+tem)*(a+tem));
		ap=az+d*am;
		bp=bz+d*bm;
		d = -(a+em)*(qab+em)*x/((qap+tem)*(a+tem));
		app=ap+d*az;
		bpp=bp+d*bz;
		aold=az;
		am=ap/bpp;
		bm=bp/bpp;
		az=app/bpp;
		bz=1.0;
		if (Math.abs(az-aold) < (EPS*Math.abs(az))) return az;
	}
	System.out.println("a or b too big, or ITMAX too small in BETACF");
    return az;
}

double betai(double a,double b,double x)
{
	//incomplete beta function
	double bt;
	if (x < 0.0 || x > 1.0) System.out.println("Bad x in routine BETAI");
	if (x == 0.0 || x == 1.0) bt=0.0;
	else
		bt=Math.exp(gammln(a+b)-gammln(a)-gammln(b)+a*Math.log(x)+b*Math.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;
}
}