//======================================================
// Numerical Analysis package in java
// error function calculation by using 
// Adaptive Gauss-Lobatto integral and serial solution
// Referans : Adapted Quadrature - Revisited
// Walter Gander - Walter Gautchi
// ftp.inf.ethz.ch
// Dr. Turhan Coban
// =====================================================

import java.io.*;
import javax.swing.*;

class f1 extends f_x
{ double func(double x)
  {return 2.0/Math.sqrt(Math.PI)*Math.exp(-x*x);}
}

class erf
{

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 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 gammln(double xx)
{
	//gamma function
	double x,y,tmp,ser;
	double cof[]={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)*Math.log(tmp);
	ser=1.000000000190015;
	for (j=0;j<=5;j++) ser += cof[j]/++y;
	return -tmp+Math.log(2.5066282746310005*ser/x);
}
public static double  erf(double x)
{
   //error function from series solution
   return x < 0.0 ? -gammp(0.5,x*x) : gammp(0.5,x*x);
}

public static double  erf_int(double x)
{
   //error function from series solution
   f1 b1=new f1();
   return adaptive_Lobatto_integral(b1,0.0,x);
}


public static double adaptive_Lobatto_integral(f_x ff,double a,double b)
{double tol=Double.MIN_VALUE;
 return adaptive_Lobatto_integral(ff,a,b,tol);
}

public static double adaptive_Lobatto_integral(f_x ff,double a,double b,double tol)
{   int i;
	double m=(a+b)/2.0; double h=(b-a)/2.0;
    double alpha=Math.sqrt(2.0/3.0); 
    double beta=1.0/Math.sqrt(5.0);
    double x1,x2,x3;
    double eps =Double.MIN_VALUE;
    //System.out.println("eps="+eps);  
    x1=0.94288241569547971905635175843185720232; 
    x2=0.64185334234578130578123554132903188354;
    x3=0.23638319966214988028222377349205292599;
    double A,B,C,D,E,F,G;
    A=0.015827191973480183087169986733305510591;
    B=0.094273840218850045531282505077108171960;
    C=0.15507198733658539625363597980210298680;
    D=0.18882157396018245442000533937297167125;
    E=0.19977340522685852679206802206648840246;
    F=0.22492646533333952701601768799639508076;
    G=0.24261107190140773379964095790325635233;
	double x[]={a,(m-x1*h),(m-alpha*h),(m-x2*h),(m-beta*h),(m-x3*h),m,(m+x3*h),(m+beta*h),(m+x2*h),(m+alpha*h),(m+x1*h),b};
	int n=x.length;
	double y[]=new double[n];
	for(i=0;i<n;i++) y[i]=ff.func(x[i]);
	double fa=y[0]; 
	double fb=y[12];
    double i2=(h/6.0)*(y[0]+y[12]+5.0*(y[4]+y[8]));
    double i1=(h/1470.0)*(77.0*(y[0]+y[12])+432.0*(y[2]+y[10])+
     625.0*(y[4]+y[8])+672.0*y[6]);
    double is=h*(A*(y[0]+y[12])+B*(y[1]+y[11])+C*(y[2]+y[10])+D*(y[3]+y[9])+E*(y[4]+y[8])+F*(y[5]+y[7])+G*y[6]);    
    int s=(int)(is/Math.abs(is)); 
    if(s==0) s=1;
    double erri1=Math.abs(i1-is);
    double erri2=Math.abs(i2-is);
    double R; 
    R=erri1/erri2;
    if(R>0 && R<1) tol=tol/R;
    //System.out.println("is="+is+"s="+s+"tol="+tol+"eps="+eps);
    if(is==0) is=b-a;
  double Q=0;  
  Q=Lobatto(ff,a,b,fa,fb,is,tol);
  return Q;	
}
	
public static double Lobatto(f_x ff,double a,double b,double fa,double fb,double is,double tol)
{
  double h=(b-a)/2; 
  double m=(a+b)/2;
  double alpha=Math.sqrt(2.0/3.0); 
  double beta=1.0/Math.sqrt(5.0);
  double mll=m-alpha*h; 
  double ml=m-beta*h; 
  double mr=m+beta*h; 
  double mrr=m+alpha*h;
  double x[]={mll,ml,m,mr,mrr};
  int n=x.length;
  double y[]=new double[n];
  for(int i=0;i<n;i++) y[i]=ff.func(x[i]);
  double fmll=y[0]; 
  double fml=y[1]; 
  double fm=y[2]; 
  double fmr=y[3]; 
  double fmrr=y[4];
  double Q;
  double i2=(h/6)*(fa+fb+5*(fml+fmr));
  double i1=(h/1470)*(77*(fa+fb)+432*(fmll+fmrr)+625*(fml+fmr)+672*fm);
  double di=Math.abs(i1-i2);
  if (((is + di) == is) || (mll <= a) || (b<=mrr) || di<tol)
  {
   if (((m <= a) || (b<=m)))
   { System.out.println("istenilen toleransa ulaşılamadı delta i = "+di);}
   Q=i1;
  }  
  else
  {
    Q=Lobatto(ff,a,mll,fa,fmll,is,tol) + Lobatto(ff,mll,ml,fmll,fml,is,tol)+
      Lobatto(ff,ml,m,fml,fm,is,tol)   + Lobatto(ff,m,mr,fm,fmr,is,tol)+
      Lobatto(ff,mr,mrr,fmr,fmrr,is,tol)+Lobatto(ff,mrr,b,fmrr,fb,is,tol);
  }
  return Q;
}

  public static void main(String args[]) throws IOException
  {
  double x=Double.parseDouble(JOptionPane.showInputDialog(" x = "));	  
  System.out.println("erf("+x+") = "+erf(x)+"integral erf="+erf_int(x));
  }
}
m