import static java.lang.Math.*;
import java.io.*;
import java.util.*;
import javax.swing.*;
import javax.swing.table.*;
import java.awt.*;
@FunctionalInterface
interface if_x 
{     
public double func(double x);
//conversion method if_x to f_x(abstract class)
default public f_x f(double x)
{ class ff extends f_x{
  public double func(double x){return func(x);}}
  ff f=new ff();
return f;
}
  

default void Plot(double a,double b)
{   Plot pp=new Plot(this,a,b);
	pp.plot();
}
default void vplot(double a,double b,int n)
{   vplot.plot(this,a,b,n);}
default void vplot(double a,double b)
{   vplot.plot(this,a,b,400);}
default void Plot(double a,double b,int n)
{   Plot pp=new Plot(this,a,b,n);
	pp.plot();
}
default double dfunc(double x)
{double h=1.0e-3;
 int n=1;
 return dfunc_n(x,n,h);
}
default double dfunc(double x,double h)
{int n=1;
 return dfunc_n(x,n,h);
}
default double dfunc(double x,int n)
{double h=1.0e-3;
 return dfunc_n(x,n,h);
}

default double dfunc_n(double x,int n,double h)
{
double hh=1/h;
double df=0;
if(n==0) df=func(x);
else if(n==1)
{df=(3.0*func(x-4.0*h)-32.0*func(x-3.0*h)+168.0*func(x-2.0*h)-672.0*func(x-h)+672.0*func(x+h)-168.0*func(x+2.0*h)+32.0*func(x+3.0*h)-3.0*func(x+4.0*h))/840.0*hh;}
else if(n==2)
{df=(-14350.0*func(x)-9.0*func(x-4*h)+128*func(x-3*h)-1008*func(x-2*h)+8064*func(x-h)+8064.0*func(x+h)-1008.0*func(x+2.0*h)+128.0*func(x+3.0*h)-9.0*func(x+4.0*h))/5040.0*hh*hh;}
else if(n==3)
{df=(-7.0*func(x-4.0*h)+72.0*func(x-3.0*h)-338.0*func(x-2.0*h)+488.0*func(x-h)-488.0*func(x+h)+338.0*func(x+2.0*h)-72.0*func(x+3.0*h)+7.0*func(x+4.0*h))/240.0*hh*hh*hh;}
else if(n==4)
{df=(2730.0*func(x)+7.0*func(x-4.0*h)-96.0*func(x-3.0*h)+676.0*func(x-2*h)-1952.0*func(x-h)-1952.0*func(x+h)+676.0*func(x+2.0*h)-96.0*func(x+3.0*h)+7.0*func(x+4.0*h))/240.0*hh*hh*hh*hh;}
else if(n==5)
{df=(func(x-4.0*h)-9.0*func(x-3.0*h)+26.0*func(x-2.0*h)-29.0*func(x-h)+29.0*func(x+h)-26.0*func(x+2.0*h)+9.0*func(x+3.0*h)-func(x+4.0*h))/6.0*hh*hh*hh*hh*hh;}
else if(n==6)
{df=(-150.0*func(x)-func(x-4.0*h)+12.0*func(x-3.0*h)-52.0*func(x-2.0*h)+116.0*func(x-h)+116.0*func(x+h)-52.0*func(x+2.0*h)+12.0*func(x+3.0*h)-func(x+4.0*h))/4.0*hh*hh*hh*hh*hh*hh;}
else if(n==7)
{df=(-func(x-4.0*h)+6.0*func(x-3.0*h)-14.0*func(x-2.0*h)+14.0*func(x-h)-14.0*func(x+h)+14.0*func(x+2.0*h)-6.0*func(x+3.0*h)+func(x+4.0*h))/2.0*hh*hh*hh*hh*hh*hh*hh;}
else if(n==8)
{df=(70.0*func(x)+func(x-4.0*h)-8.0*func(x-3.0*h)+28.0*func(x-2.0*h)-56.0*func(x-h)-56.0*func(x+h)+28.0*func(x+2.0*h)-8.0*func(x+3.0*h)+func(x+4.0*h))*hh*hh*hh*hh*hh*hh*hh*hh;}
else df=0;
return df;
}
//INTEGRAL
default double integral(double a,double b)
{ //10 point Gauss-Legendre formula
//integral f(x)dx
 double r[]={-0.973906528517171,-0.865063366688984,-0.679409568299024,-0.433395394129247,-0.148874338981631,
 0.148874338981631,0.433395394129247,0.679409568299024,0.865063366688984,0.973906528517171};
 double c[]={0.066671344308684,0.149451349150580,0.219086362515982,0.269266719309996,0.295524224714752,
 0.295524224714752,0.269266719309996,0.219086362515982,0.149451349150580,0.066671344308684};
double z=0,x,y;
double k1=(b-a)/2.0;
double k2=(b+a)/2.0;
for(int i=0;i<r.length;i++)
{
x=k2+k1*r[i];
y=func(x);
z+=k1*c[i]*y;
}
return z;
}
default double[][] gauss_legendre_coefficients(double x1,double x2,int n)
{
	//calculates legendre gauss-coefficients as coefficients of the integral
	//for n terms
	double EPS=3.0e-15;
	int m,j,i;
	double z1,z,xm,xl,pp,p3,p2,p1;
    //double x[]=new double[n];
    //double w[]=new double[n];
    double a[][]=new double[2][n];//a[0][i]=x[i]   a[1][i]=w[i]
	m=(n+1)/2;
	xm=0.5*(x2+x1);
	xl=0.5*(x2-x1);
	for (i=1;i<=m;i++)  {
		z=Math.cos(Math.PI*((i-0.25)/(n+0.5)));
		do {
			p1=1.0;
			p2=0.0;
			for (j=1;j<=n;j++) {
				p3=p2;
				p2=p1;
				p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
			}
			pp=n*(z*p1-p2)/(z*z-1.0);
			z1=z;
			z=z1-p1/pp;
		} while (Math.abs(z-z1) > EPS);
		a[0][i-1]=xm-xl*z;
		a[0][n-i]=xm+xl*z;
		a[1][i-1]=2.0*xl/((1.0-z*z)*pp*pp);
		a[1][n-i]=a[1][i-1];
			}
return a;
}

default double integral(double x1,double x2,int n)
{
// n : number of integral coefficients
// this routine first generates gauss legendre coefficients
// for [x1,x2] band 
// then calculates gauss legendre integral
double a[][]=new double[2][n];
a=gauss_legendre_coefficients(x1,x2,n);
Text.print(Text.T(a));
//System.out.println(Matrix.toString(a));
double z=0;
for(int i=0;i<n;i++)
{ z+=a[1][i]*func(a[0][i]);}
return z;
}
    default double [] lngamma(double x)
    {//logaritmik hata fonksiyonu
    //xx[0] hata fonksiyonunun logaritmasinin değerini 
    //xx[1] hata fonksiyonunun kendisinin işaretini verir
	double xx[]=new double[2];    
    double result;
    double a;
    double b;
    double c;
    double p;
    double q;
    double u;
    double w;
    double z;
    int i;
    double logpi;
    double ls2pi;
    double tmp;
    double sgngam = 1;
    logpi = 1.14472988584940017414;
    ls2pi = 0.91893853320467274178;
    if( x<-34.0 )
    {
        q = -x;
        double yy[]=lngamma(q);
        w = yy[0];
        tmp=yy[1];
        p = Math.floor(q);
        i = (int)Math.round(p);
        if( i%2==0 )
        {
            sgngam = -1;
        }
        else
        {
            sgngam = 1;
        }
        z = q-p;
        if( z>0.5 )
        {
            p = p+1;
            z = p-q;
        }
        z = q*Math.sin(Math.PI*z);
        result = logpi-Math.log(z)-w;
        xx[0]=result;
        xx[1]=sgngam;
        return xx;
    }
    if( x<13 )
    {
        z = 1;
        p = 0;
        u = x;
        while(u>=3)
        {
            p = p-1;
            u = x+p;
            z = z*u;
        }
        while(u<2)
        {
            z = z/u;
            p = p+1;
            u = x+p;
        }
        if( z<0 )
        {
            sgngam = -1;
            z = -z;
        }
        else
        {
            sgngam = 1;
        }
        if( u==2 )
        {
            result = Math.log(z);
            //return result;
            xx[0]=result;
            xx[1]=sgngam;
            return xx;
        }
        p = p-2;
        x = x+p;
        b = -1378.25152569120859100;
        b = -38801.6315134637840924+x*b;
        b = -331612.992738871184744+x*b;
        b = -1162370.97492762307383+x*b;
        b = -1721737.00820839662146+x*b;
        b = -853555.664245765465627+x*b;
        c = 1;
        c = -351.815701436523470549+x*c;
        c = -17064.2106651881159223+x*c;
        c = -220528.590553854454839+x*c;
        c = -1139334.44367982507207+x*c;
        c = -2532523.07177582951285+x*c;
        c = -2018891.41433532773231+x*c;
        p = x*b/c;
        result = Math.log(z)+p;
        //return result;
        xx[0]=result;
        xx[1]=sgngam;
        return xx;
    }
    q = (x-0.5)*Math.log(x)-x+ls2pi;
    if( x>100000000 )
    {
        result = q;
        //return result;
        xx[0]=result;
        xx[1]=sgngam;
        return xx;
    }
    p = 1/(x*x);
    if( x>=1000.0 )
    {
        q = q+((7.9365079365079365079365*0.0001*p-2.7777777777777777777778*0.001)*p+0.0833333333333333333333)/x;
    }
    else
    {
        a = 8.11614167470508450300*0.0001;
        a = -5.95061904284301438324*0.0001+p*a;
        a = 7.93650340457716943945*0.0001+p*a;
        a = -2.77777777730099687205*0.001+p*a;
        a = 8.33333333333331927722*0.01+p*a;
        q = q+a/x;
    }
    xx[0]= q;
    xx[1]=sgngam;
    return xx;
    }
    	
default double[][] gauss_laguerre_coefficients(int n,double alpha)
{   //Gauss-Laguerre quadratürü katsayilarini hesaplar
    double sgngam=1;
    double a[][]=new double[2][n];
    int i;
    int j;
    double r=0; 
    double r1;
    double p1;
    double p2;
    double p3;
    double dp3;
    double tsg;
    double eps=1.0e-10;
    for(i = 0; i <= n-1; i++)
    {
        if( i==0 )
        {
            r = (1.0+alpha)*(3.0+0.92*alpha)/(1+2.4*n+1.8*alpha);
        }
        else
        {
            if( i==1 )
            {r = r+(15.0+6.25*alpha)/(1.0+0.9*alpha+2.5*n);}
            else
            {r = r+((1.0+2.55*(i-1))/(1.9*(i-1))+1.26*(i-1)*alpha/(1+3.5*(i-1)))/(1.0+0.3*alpha)*(r-a[0][i-2]);}
        }
        do
        {
            p2 = 0;
            p3 = 1;
            for(j = 0; j <= n-1; j++)
            {
                p1 = p2;
                p2 = p3;
                p3 = ((-r+2*j+alpha+1)*p2-(j+alpha)*p1)/(j+1);
            }
            dp3 = (n*p3-(n+alpha)*p2)/r;
            r1 = r;
            r = r-p3/dp3;
        }
        while(Math.abs(r-r1)>=eps*(1+Math.abs(r))*100);
        a[0][i] = r;
        double lng[]=lngamma((double)n);
        double lng1[]=lngamma((double)(alpha+n));
        a[1][i] = -Math.exp(lng1[0]-lng[0])/(dp3*n*p2);
        tsg=lng[1];
    }    
return a;
}

default double gauss_laquerre_integral(double xa,int n,double alpha)
{
//n : number of integral coefficients
// this routine first generates gauss legendre coefficients
// for [x1,x2] band 
// then calculates gauss hermite integral
double a[][]=new double[2][n];
a=gauss_laguerre_coefficients(n,alpha);
//System.out.println(Matrix.toString(a));
double z=0;
for(int i=0;i<n;i++)
{ z+=a[1][i]*func(a[0][i]+xa)*Math.exp(a[0][i])*Math.pow(a[0][i],-alpha);}
return z;
}

default double gauss_laquerre_integral(double xa,int n)
{
double alpha=1.0;	
double a[][]=new double[2][n];
a=gauss_laguerre_coefficients(n,alpha);
//System.out.println(Matrix.toString(a));
double z=0;
for(int i=0;i<n;i++)
{ z+=a[1][i]*func(a[0][i]+xa)*Math.exp(a[0][i])/a[0][i];}
return z;
}
default double adaptive_simpson_integral(double a,double b)
{return adaptive_simpson_integral(a,b,1.0e-20,30);}

default double adaptive_simpson_integral(double a,double b,double eps)
{return adaptive_simpson_integral(a,b,eps,30);}

default double adaptive_simpson_integral(double a,double b,int MAX_ITERATION)
{return adaptive_simpson_integral(a,b,1.0e-12,MAX_ITERATION);}

default double adaptive_simpson_integral(double a,double b,double eps,int MAX_ITERATION)
{   int i;
    double h=(b-a);
    double h2=h/4.0;
	double x[]=new double[5];
    double f[]=new double[5];
    for(i=0;i<5;i+=2) {x[i]=a+i*h2;f[i]=func(x[i]);}
	double sum=h*(f[0]+4.0*f[2]+f[4])/6.0; 
	double Q=0;
    Q=simpson1_3(x,f,sum,eps,0,MAX_ITERATION,Q);
    return Q;	
}
default double adaptive_gauss_legendre_integral(double x1,double x2)
{return adaptive_gauss_legendre_integral(x1,x2,1.0e-12,100);}

default double adaptive_gauss_legendre_integral(double x1,double x2,double eps)
{return adaptive_gauss_legendre_integral(x1,x2,eps,100);}

default double adaptive_gauss_legendre_integral(double x1,double x2,double eps,double MAXITER)
{
int n=1;
double ans1=integral(x1,x2,n);
double ans2;
do
{ans2=ans1;n++;ans1=integral(x1,x2,n);}
while(Math.abs(ans2-ans1)>eps && n<MAXITER);
if(n==MAXITER) System.out.println("istenilen hata seviyesi tutturulamadi n="+n );
return ans2;
}
default double gauss_konrad_integral(double a,double b,double eps,int n)
    {
	    //eps hata miktari
        int Aused=0;
	    boolean sonuc;
        double integral=0;
        double A[][];
        int Aw = 0;
        double toplamhata = 0;
        double[] c;
        double[] xg;
        double[] xk;
        int nn = 0;
        int ng = 0;
        int i = 0;
        int j = 0;
        int h = 0;
        double v = 0;
        double k1 = 0;
        double k2 = 0;
        double intg = 0;
        double intk = 0;
        double ta = 0;
        double tb = 0;
        Aw = 4;
        A = new double[n-1+1][Aw-1+1];
        nn = 61;
        ng = 15;
        c = new double[nn];
        xk = new double[nn];
        xg = new double[nn];
        //Gauss katsayilari
        xg[0] = 0.007968192496166605615465883474674;
        xg[1] = 0.018466468311090959142302131912047;
        xg[2] = 0.028784707883323369349719179611292;
        xg[3] = 0.038799192569627049596801936446348;
        xg[4] = 0.048402672830594052902938140422808;
        xg[5] = 0.057493156217619066481721689402056;
        xg[6] = 0.065974229882180495128128515115962;
        xg[7] = 0.073755974737705206268243850022191;
        xg[8] = 0.080755895229420215354694938460530;
        xg[9] = 0.086899787201082979802387530715126;
        xg[10] = 0.092122522237786128717632707087619;
        xg[11] = 0.096368737174644259639468626351810;
        xg[12] = 0.099593420586795267062780282103569;
        xg[13] = 0.101762389748405504596428952168554;
        xg[14] = 0.102852652893558840341285636705415;
        c[0] = 0.999484410050490637571325895705811;
        c[1] = 0.996893484074649540271630050918695;
        c[2] = 0.991630996870404594858628366109486;
        c[3] = 0.983668123279747209970032581605663;
        c[4] = 0.973116322501126268374693868423707;
        c[5] = 0.960021864968307512216871025581798;
        c[6] = 0.944374444748559979415831324037439;
        c[7] = 0.926200047429274325879324277080474;
        c[8] = 0.905573307699907798546522558925958;
        c[9] = 0.882560535792052681543116462530226;
        c[10] = 0.857205233546061098958658510658944;
        c[11] = 0.829565762382768397442898119732502;
        c[12] = 0.799727835821839083013668942322683;
        c[13] = 0.767777432104826194917977340974503;
        c[14] = 0.733790062453226804726171131369528;
        c[15] = 0.697850494793315796932292388026640;
        c[16] = 0.660061064126626961370053668149271;
        c[17] = 0.620526182989242861140477556431189;
        c[18] = 0.579345235826361691756024932172540;
        c[19] = 0.536624148142019899264169793311073;
        c[20] = 0.492480467861778574993693061207709;
        c[21] = 0.447033769538089176780609900322854;
        c[22] = 0.400401254830394392535476211542661;
        c[23] = 0.352704725530878113471037207089374;
        c[24] = 0.304073202273625077372677107199257;
        c[25] = 0.254636926167889846439805129817805;
        c[26] = 0.204525116682309891438957671002025;
        c[27] = 0.153869913608583546963794672743256;
        c[28] = 0.102806937966737030147096751318001;
        c[29] = 0.051471842555317695833025213166723;
        c[30] = 0.000000000000000000000000000000000;
        //Kronrod katsayilari
        xk[0] = 0.001389013698677007624551591226760;
        xk[1] = 0.003890461127099884051267201844516;
        xk[2] = 0.006630703915931292173319826369750;
        xk[3] = 0.009273279659517763428441146892024;
        xk[4] = 0.011823015253496341742232898853251;
        xk[5] = 0.014369729507045804812451432443580;
        xk[6] = 0.016920889189053272627572289420322;
        xk[7] = 0.019414141193942381173408951050128;
        xk[8] = 0.021828035821609192297167485738339;
        xk[9] = 0.024191162078080601365686370725232;
        xk[10] = 0.026509954882333101610601709335075;
        xk[11] = 0.028754048765041292843978785354334;
        xk[12] = 0.030907257562387762472884252943092;
        xk[13] = 0.032981447057483726031814191016854;
        xk[14] = 0.034979338028060024137499670731468;
        xk[15] = 0.036882364651821229223911065617136;
        xk[16] = 0.038678945624727592950348651532281;
        xk[17] = 0.040374538951535959111995279752468;
        xk[18] = 0.041969810215164246147147541285970;
        xk[19] = 0.043452539701356069316831728117073;
        xk[20] = 0.044814800133162663192355551616723;
        xk[21] = 0.046059238271006988116271735559374;
        xk[22] = 0.047185546569299153945261478181099;
        xk[23] = 0.048185861757087129140779492298305;
        xk[24] = 0.049055434555029778887528165367238;
        xk[25] = 0.049795683427074206357811569379942;
        xk[26] = 0.050405921402782346840893085653585;
        xk[27] = 0.050881795898749606492297473049805;
        xk[28] = 0.051221547849258772170656282604944;
        xk[29] = 0.051426128537459025933862879215781;
        xk[30] = 0.051494729429451567558340433647099;
        for(i=nn-1; i>=nn/2; i--)
        {c[i] = -c[nn-1-i];}
        for(i=nn-1; i>=nn/2; i--)
        { xk[i] = xk[nn-1-i]; }
        for(i=ng-1; i>=0; i--)
        { xg[nn-2-2*i] = xg[i];
          xg[1+2*i] = xg[i];
        }
        for(i=0; i<=nn/2; i++)
        { xg[2*i] = 0; }
        k1 = 0.5*(b-a);
        k2 = 0.5*(b+a);
        intg = 0;
        intk = 0;
        for(i=0; i<=nn-1; i++)
        { v = func(k1*c[i]+k2);
          intk = intk+v*xk[i];
          if( i%2==1 )
          {intg = intg+v*xg[i];}
        }
        intk = intk*(b-a)*0.5;
        intg = intg*(b-a)*0.5;
        A[0][0] = Math.abs(intg-intk);
        A[0][1] = intk;
        A[0][2] = a;
        A[0][3] = b;
        toplamhata = A[0][0];
        if( toplamhata<eps )
        {
            sonuc = true;
            integral = intk;
            Aused = 1;
            return integral;
        }
        Aused = 1;
        for(h=1; h<=n-1; h++)
        {
            Aused = h+1;
            gir(A, h, Aw);
            toplamhata = toplamhata-A[h-1][0];
            ta = A[h-1][2];
            tb = A[h-1][3];
            A[h-1][2] = ta;
            A[h-1][3] = 0.5*(ta+tb);
            A[h][2] = 0.5*(ta+tb);
            A[h][3] = tb;
            for(j=h-1; j<=h; j++)
            {
                k1 = 0.5*(A[j][3]-A[j][2]);
                k2 = 0.5*(A[j][3]+A[j][2]);
                intg = 0;
                intk = 0;
                for(i=0; i<=nn-1; i++)
                {   v = func(k1*c[i]+k2);
                    intk = intk+v*xk[i];
                    if( i%2==1 )
                    { intg = intg+v*xg[i];}
                }
                intk = intk*(A[j][3]-A[j][2])*0.5;
                intg = intg*(A[j][3]-A[j][2])*0.5;
                A[j][0] = Math.abs(intg-intk);
                A[j][1] = intk;
                toplamhata = toplamhata+A[j][0];
            }
            cik(A, h-1, Aw);
            cik(A, h, Aw);
            if( toplamhata<eps )
            {break;}
        }
        sonuc = toplamhata<eps;
        integral = 0;
        for(j=0; j<=Aused-1; j++)
        { integral = integral+A[j][1];}
        return integral;
    }

   default void gir(double A[][],int n,int Awidth)
    {
        int i = 0;
        int p = 0;
        double t = 0;
        int maxcp = 0;
        if( n==1 )
        { return; }
        for(i=0; i<=Awidth-1; i++)
        {   t = A[n-1][i];
            A[n-1][i] = A[0][i];
            A[0][i] = t;
        }//end of for
        p = 0;
        while( 2*p+1<n-1 )
        {   maxcp = 2*p+1;
            if( 2*p+2<n-1 )
            {   if( A[2*p+2][0]>A[2*p+1][0] )
                { maxcp = 2*p+2; }
            }
            if( A[p][0]<A[maxcp][0] )
            {   for(i=0; i<=Awidth-1; i++)
                {   t = A[p][i];
                    A[p][i] = A[maxcp][i];
                    A[maxcp][i] = t;
                }
                p = maxcp;
            }
            else
            { break; }
        } //end of while
    }

    default void cik(double A[][],int n,int Awidth)
    {
        int i = 0;
        int p = 0;
        double t = 0;
        int kk = 0;
        if( n==0 )
        {return;}
        p = n;
        while( p!=0 )
        {   kk = (p-1)/2;
            if( A[p][0]>A[kk][0] )
            { for(i=0; i<=Awidth-1; i++)
                { t = A[p][i];
                  A[p][i] = A[kk][i];
                  A[kk][i] = t;
                }
                p = kk;
            }
            else
            {break;}
        }
    }
default double adaptive_simpson_integral2(double a,double b)
{   int i;
	Plot(a,b);
    double h=(b-a);
    double h2=h/4.0;
    double h3=h/8.0;
    double m=(a+b)/2.0;
    double x[]=new double[9];
    double y[]=new double[9];
    for(i=0;i<9;i++) {x[i]=a+i*h3;y[i]=func(x[i]);}
    double is=h/28350.0*(989.0*y[0]+5888.0*y[1]-928.0*y[2]+10496.0*y[3]-4540.0*y[4]+
           10496.0*y[5]-928.0*y[6]+5888.0*y[7]+989.0*y[8]); 
    double fa = y[0]; 
    double fm = y[4]; 
    double fb = y[8];
    if(is==0) is = b-a;
    double Q = adaptsimstp(a,b,fa,fm,fb,is);
    return Q;	
}
	
default double adaptsimstp(double a,double b,double fa,double fm,double fb,double is)
{
  double m = (a + b)/2; 
  double h = (b - a)/4;
  double x[] = {a + h, b - h};
  int n=x.length;
  double y[] = new double[n];
  for(int i=0;i<n;i++) {y[i]=func(x[i]);}
  double fml = y[0];
  double fmr = y[1];
  double i1 = h/1.5 * (fa + 4*fm + fb);
  double i2 = h/3 * (fa + 4*(fml + fmr) + 2*fm + fb);
  double Q=0;
  i1 = (16.0*i2 - i1)/15.0; //romberg integrasyonu ile iyileştirme
  if (((is + (i1-i2)) == is) || (m <= a) || (b<=m))
  {
    if ((m <= a) || (b<=m))
       { System.out.println("required tolerance is not  reached");}
    Q = i1;
  }  
  else
  {
    Q = adaptsimstp(a,m,fa,fml,fm,is)+
        adaptsimstp(m,b,fm,fmr,fb,is);
  }
  return Q;
}
default double adaptive_Lobatto_integral(double a,double b)
{   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;
    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]=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;
    //is=s*Math.abs(is)*tol/eps; 
    if(is==0) is=b-a;
  double Q=0;  
  Q=Lobatto(a,b,fa,fb,is);
  return Q;	
}
	
default double Lobatto(double a,double b,double fa,double fb,double is)
{
  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]=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);
  boolean c1,c2,c3;
  double w1=i1-i2;;
  double w2=is;
  double w3=w2+w1;  
  c1=(w3==w2);
  c2=(mll<=a);
  c3=(b<=mrr);
  if( c1|| c2 ||c3 )
  {
   if (((m <= a) || (b<=m)))
   {
       System.out.println("istenilen toleransa ulaşilamadi");
   }
   Q=i1;
  }  
  else
  {
    Q=Lobatto(a,mll,fa,fmll,is) + Lobatto(mll,ml,fmll,fml,is)+
      Lobatto(ml,m,fml,fm,is)   + Lobatto(m,mr,fm,fmr,is)+
      Lobatto(mr,mrr,fmr,fmrr,is)+Lobatto(mrr,b,fmrr,fb,is);
  }
  return Q;
}
	
default double simpson1_3(double x[],double f[],double sum,double eps,int k,int MAX_ITERATION,double Q)
{       //Newton-cotes integral 3 points
        //adaptive Simpson 1/3 rule
        double h=(x[4]-x[0]);
        double h1=h/2.0;
        double h2=h/4.0;        
        int p,i;
        double x1[]=new double[5];
        double f1[]=new double[5];
        double x2[]=new double[5];
        double f2[]=new double[5];
        k++;
        for(i=1;i<5;i+=2) {x[i]=x[0]+i*h2;f[i]=func(x[i]);}
        double sum1=h1*(f[0]+4.0*f[1]+f[2])/6.0;
        double sum2=h1*(f[2]+4.0*f[3]+f[4])/6.0;     
        x1[0]=x[0];f1[0]=f[0];
        x1[2]=x[1];f1[2]=f[1];
        x1[4]=x[2];f1[4]=f[2];
        x2[0]=x[2];f2[0]=f[2];
        x2[2]=x[3];f2[2]=f[3];
        x2[4]=x[4];f2[4]=f[4];      
        if(1.0/15.0*Math.abs(sum1+sum2-sum)<eps || k>MAX_ITERATION) Q+=sum1+sum2;
        else Q+=simpson1_3(x1,f1,sum1,eps,k,MAX_ITERATION,Q)+simpson1_3(x2,f2,sum2,eps,k,MAX_ITERATION,Q);
        return Q;
}

default double bisection(double xl,double xu)
{
//bisection root finding method
double test;
double xr=0;
double es,ea;
double fxl,fxr;
int maxit=100,iter=0;
es=0.0000001;
ea=1.1*es;
while((ea>es)&&(iter<maxit))
{
    xr=(xl+xu)/2.0;
    iter++;
    if((xl+xu)!=0)
        { ea=Math.abs((xu-xl)/(xu+xl))*100;}
    fxl= func(xl);
    fxr= func(xr);
    test= fxl*fxr;
    if(test==0.0)      ea=0;
    else if(test<0.0)  xu=xr;
    else              {xl=xr;}
} 
if(iter>=maxit) JOptionPane.showMessageDialog(null,"Maximum number of iteration is exceeded \n"+
   " result might not be valid","MAKSİMUM NUMBER OF ITERATION WARNING",JOptionPane.WARNING_MESSAGE); 
return xr;
}
default double illinois(double xa,double xu)
{
//Modified False position root finding: Illinois method
double test;
double p=0;
double es,ea;
double fa,fp,fu;
int maxit=100,iter=0;
es=0.000001;
ea=1.1*es;
double xold;
int ia=0,iu=0;
fa=func(xa);
fu=func(xu);
p=xu;
double xpold;
while(iter<maxit && ea>es)
{
xpold=p;	
p=xu-fu*(xa-xu)/(fa-fu);
fp=func(p);
if(Math.signum(fp)!=Math.signum(fu))
{xa=xu;fa=fu;}
else
{fa=fa/2;}
xu=p;fu=fp;
iter++;
if(p!=0) ea=Math.abs((p-xpold)/p)*100.0;
}
if(iter>=maxit) JOptionPane.showMessageDialog(null,"Maximum number of iteration is exceeded \n"+
   " result might not be valid","MAKSİMUM NUMBER OF ITERATION WARNING",JOptionPane.WARNING_MESSAGE); 
return p;
}
default double newton(double x)
{
int nmax=100;
double tolerance=1.0e-10;
double fx,dfx;
for(int i=0;i<nmax;i++)
{
  fx=func(x);
  dfx=dfunc(x);
  x-=fx/dfx;
  if(Math.abs(fx)<tolerance) { return x;}
}
JOptionPane.showMessageDialog(null,"Warning : Maximum number of iteration is exceeded \n"+
   " çözüm geçerli olmiyabilir","MAXİMUM ITERATION NUMBER WARNING",JOptionPane.WARNING_MESSAGE);
return x;
}
default double muller(f_x f,double x0,double x1,double x2)
{
//Finding real roots by using Muller method
//orta xr noktasindan +/- h kadar giderek diğer 2 noktaya ulaşiyoruz.
int maxit=5;
double xr=x0;
int iter=0;
double es1=0.001;
double es=es1;
double ea=1.1*es;
double h0,h1;
double d0,d1;
double a,b,c;
double fx0,fx1,fx2;
double determinant;
double den;
double dxr;
while(iter<maxit && ea>es)
{
    fx0=f.func(x0);
    fx1=f.func(x1);
    fx2=f.func(x2);
    h0=x1-x0;
    h1=x2-x1;
    d0=(fx1-fx0)/h0;
    d1=(fx2-fx1)/h1;
    a=(d1-d0)/(h1+h0);
    b=a*h1+d1;
    c=fx2;
    determinant=b*b-4.0*a*c;
    if(determinant<=0) determinant=0;//we do not want to calculate complex roots here
    else determinant=Math.sqrt(determinant);
    if(Math.abs(b+determinant)>Math.abs(b-determinant)) {den=b+determinant;;}
    else {den=b-determinant;}
    dxr=-2*c/den;
    xr=x2+dxr;
    System.out.println("a="+a+"b="+b+"c="+c+"disc="+determinant+"xr="+xr+"iter="+iter);
    ea=Math.abs((xr-x2)/xr)*100;
    iter++;
    x0=x1;
    x1=x2;
    x2=xr;
} 
if(iter>=maxit) JOptionPane.showMessageDialog(null,"Warning : MAximum number of iteration is exceeded \n"+
   " çözüm geçerli olmiyabilir","MAXİMUM ITERATION NUMBER WARNING",JOptionPane.WARNING_MESSAGE);
return xr;
}
default double newton_bisection(double x1, double x2)
{
double xacc=1.0e-10; //accuracy
double h=0.0001; //derivative function step
int MAXIT=100;  //Maximumnumber of iterations
int j;
double dfr,dx,dxold,fr,fh,fl;
double temp,b,a,r;
fl=func(x1);
fh=func(x2);
//Check if a root is existed in the given region
if ((fl > 0.0 && fh > 0.0) || (fl < 0.0 && fh < 0.0))
System.out.println("There are no root in the given region or double root?");
if (fl == 0.0) return x1;
if (fh == 0.0) return x2;
if (fl < 0.0)   { a=x1;b=x2;} 
else            { b=x1;a=x2;}
r=0.5*(x1+x2); //midpoint value
dxold=Math.abs(x2-x1); 
dx=dxold; 
fr=func(r);   //function value at midpoint
dfr=dfunc(r,h); //derivative of function value at midpoint 
for (j=1;j<=MAXIT;j++) 
{ 
if ((((r-b)*dfr-fr)*((r-a)*dfr-fr) > 0.0) || (Math.abs(2.0*fr) > Math.abs(dxold*dfr))) 
    { //bisection step
	  dxold=dx;
      dx=0.5*(b-a);
      r=a+dx;
      if (a == r) return r; //çözüm!!!!
    } 
else{ //Newton-Raphson (Secant) step
	  dxold=dx;
      dx=fr/dfr;
      temp=r;
      r -= dx;
      if (temp == r) return r;  //solution!!!!!
    }
if (Math.abs(dx) < xacc) return r; //solution!!!!!
fr=func(r);
dfr=dfunc(r,h);
if (fr < 0.0) 
a=r;
else
b=r;
}
System.out.println("Maximum number of iterations are exceeded");
return 0.0; //program should never reach here
}
default double brent(double x1,double x2)
{
    // roots of the nonlinear equation: Brent method
    // reference : Numerical Recipes in C, second edition, William H. Press, 
    // Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery
    // cambridge university press	
	double ITMAX=100;
               double EPS=2.22044604925031e-16;
	int iter;
	double a=x1,b=x2,c,d,e,min1,min2;
	double tol=1.0e-5;
	c=0;
	d=0;
	e=0;
	double fa=func(a),fb=func(b),fc,p,q,r,s,tol1,xm;
	if (fb*fa > 0.0) JOptionPane.showMessageDialog(null,"Root limit selection is uncorrect \n result may be wrong"
	,"ROOT LIMITS WARNING",JOptionPane.WARNING_MESSAGE);
	fc=fb;
	for (iter=1;iter<=ITMAX;iter++) {
		if (fb*fc > 0.0) {
			c=a;
			fc=fa;
			e=d=b-a;
		}
		if (Math.abs(fc) < Math.abs(fb)) {
			a=b;
			b=c;
			c=a;
			fa=fb;
			fb=fc;
			fc=fa;
		}
		tol1=2.0*EPS*Math.abs(b)+0.5*tol;
		xm=0.5*(c-b);
		if (Math.abs(xm) <= tol1 || fb == 0.0) return b;
		if (Math.abs(e) >= tol1 && Math.abs(fa) > Math.abs(fb)) {
			s=fb/fa;
			if (a == c) {
				p=2.0*xm*s;
				q=1.0-s;
			} else {
				q=fa/fc;
				r=fb/fc;
				p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
				q=(q-1.0)*(r-1.0)*(s-1.0);
			}
			if (p > 0.0)  q = -q;
			p=Math.abs(p);
			min1=3.0*xm*q-Math.abs(tol1*q);
			min2=Math.abs(e*q);
			if (2.0*p < (min1 < min2 ? min1 : min2)) {
				e=d;
				d=p/q;
			} else {
				d=xm;
				e=d;
			}
		} else {
			d=xm;
			e=d;
		}
		a=b;
		fa=fb;
		if (Math.abs(d) > tol1)
			b += d;
		else
			b += (xm > 0.0 ? Math.abs(tol1) : -Math.abs(tol1));
		fb=func(b);
	}
JOptionPane.showMessageDialog(null,"Maximum number of iterations are exceeded \n"+
   " Result mAy not be valid","MAXIMUM ITERATION NUMBER WARNING",JOptionPane.WARNING_MESSAGE);
return b;
}
default double halley(double x)
{
int nmax=500;
double tolerance=1.0e-8;
double fx,dfx,fxm,dfxm,d2fx;
for(int i=0;i<nmax;i++)
{  fx=func(x);
   dfx=dfunc(x);
   d2fx=dfunc(x,2);
   x-=2.0*fx*dfx/(2.0*dfx*dfx-fx*d2fx);
  if(Math.abs(fx)<tolerance) { return x;}
}
JOptionPane.showMessageDialog(null,"Maximum number of iteration is exceeded in Halley method\n"+
   " results may not be valid","MAXIMUM ITERATION WARNING",JOptionPane.WARNING_MESSAGE);
return x;
}
default double abbasbandy(double x)
{
int nmax=500;
double tolerance=1.0e-8;
double fx,dfx,d2fx,d3fx;
for(int i=0;i<nmax;i++)
{  fx=func(x);
   dfx=dfunc(x);
   d2fx=dfunc(x,2);
   d3fx=dfunc(x,3);
   x-=(fx/dfx+fx*fx*d2fx/(2.0*dfx*dfx*dfx)+fx*fx*fx*d3fx/(6.0*dfx*dfx*dfx*dfx));
  if(Math.abs(fx)<tolerance) { return x;}
}
  JOptionPane.showMessageDialog(null,"Maximum number of iteration is exceeded in Abbasbandy\n"+
   " results may not be valid","MAXIMUM ITERATION WARNING",JOptionPane.WARNING_MESSAGE);
   return x;
}   

default double TM1(double x)
{ //Thiele 1 method
int nmax=500;
double tolerance=1.0e-8;
double fx,dfx,d2fx,d3fx;
double P1,Q1;
for(int i=0;i<nmax;i++)
{  fx=func(x);
   dfx=dfunc(x);
   d2fx=dfunc(x,2);
   d3fx=dfunc(x,3);
   P1=4.0*fx*dfx*(3.0*dfx*dfx*d2fx-3.0*fx*d2fx*d2fx+fx*dfx*d3fx);
   Q1=12.0*dfx*dfx*dfx*dfx*d2fx-18.0*fx*dfx*dfx*d2fx*d2fx+3.0*fx*fx*d2fx*d2fx*d2fx+4.0*fx*dfx*dfx*dfx*d3fx;
   x-=P1/Q1;
  if(Math.abs(fx)<tolerance) { return x;}
}
JOptionPane.showMessageDialog(null,"Maximum number of iteration is exceeded in Thiele 1\n"+
   " results may not be valid","MAXIMUM ITERATION WARNING",JOptionPane.WARNING_MESSAGE);
return x;
}

default double TM2(double x)
{ //Thiele 2 method
int nmax=500;
double tolerance=1.0e-8;
double fx,dfx,d2fx,d3fx;
double P2,Q2;
for(int i=0;i<nmax;i++)
{  fx=func(x);
   dfx=dfunc(x);
   d2fx=dfunc(x,2);
   d3fx=dfunc(x,3);
   P2=fx*(6.0*dfx*dfx*d2fx-3.0*fx*d2fx*d2fx+2*fx*dfx*d3fx);
   Q2=2.0*dfx*(3.0*dfx*dfx*d2fx-3.0*fx*d2fx*d2fx+fx*dfx*d3fx);
   x-=P2/Q2;
  if(Math.abs(fx)<tolerance) { return x;}
}
JOptionPane.showMessageDialog(null,"Maximum number of iteration is exceeded in Thiele 2 method\n"+
   " results may not be valid","MAXIMUM ITERATION WARNING",JOptionPane.WARNING_MESSAGE);
return x;
}
default double hasan47a(double x)
{ // third order
int nmax=500;
double tolerance=1.0e-8;
double fx,dfx;
double x1=0,dfx1=0;
for(int i=0;i<nmax;i++)
{  fx=func(x);
   dfx=dfunc(x);
   x1=x-fx/(2.0*dfx);
   dfx1=dfunc(x1);
   x-=fx/dfx1;
  if(Math.abs(fx)<tolerance) { return x;}
}
JOptionPane.showMessageDialog(null,"Maximum number of iteration is exceeded in hasan47a Method\n"+
   " results may not be valid","MAXIMUM ITERATION WARNING",JOptionPane.WARNING_MESSAGE);
return x;
}

default double hasan47b(double x)
{ //third order
int nmax=500;
double tolerance=1.0e-8;
double fx,dfx;
double x1=0,dfx1=0;
for(int i=0;i<nmax;i++)
{  fx=func(x);
   dfx=dfunc(x);
   x1=x-fx/dfx;
   dfx1=dfunc(x1);
   x-=2.0*fx/(dfx+dfx1);
  if(Math.abs(fx)<tolerance) { return x;}
}
JOptionPane.showMessageDialog(null,"Maximum number of iteration is exceeded in hasan47b Method\n"+
   " results may not be valid","MAXIMUM ITERATION WARNING",JOptionPane.WARNING_MESSAGE);
return x;
}
default double hasan47c(double x,double alfa)
{ //third order
int nmax=500;
double beta=1-alfa;
double tolerance=1.0e-8;
double fx,dfx;
double x1=0,dfx1=0;
double x2=0,dfx2=0;
for(int i=0;i<nmax;i++)
{  fx=func(x);
   dfx=dfunc(x);
   x1=x-alfa*fx/dfx;
   dfx1=dfunc(x1);
   x2=x-beta*fx/dfx;
   dfx2=dfunc(x2);
   x-=2.0*fx/(dfx1+dfx2);
  if(Math.abs(fx)<tolerance) { return x;}
}

JOptionPane.showMessageDialog(null,"Maximum number of iteration is exceeded in hasan47c Method\n"+
   " results may not be valid","MAXIMUM ITERATION WARNING",JOptionPane.WARNING_MESSAGE);
return x;
}

default double hasan43(double x)
{ //fourth order
int r=3;
int nmax=500;
double tolerance=1.0e-8;
double fx,dfx,d2fx;
for(int i=0;i<nmax;i++)
{  fx=func(x);
   dfx=dfunc(x);
   d2fx=dfunc(x,2);
   x-=fx/(dfx+1.0/r*dfx*Math.log((dfx*dfx-fx*d2fx)/(dfx*dfx)));
  if(Math.abs(fx)<tolerance) { return x;}
}
JOptionPane.showMessageDialog(null,"Maximum number of iteration is exceeded in hasan43 Method\n"+
   " results may not be valid","MAXIMUM ITERATION WARNING",JOptionPane.WARNING_MESSAGE);
return x;
}
default double SIGN(double a,double b)
{
return Math.abs(a)*b/Math.abs(b);
}
default double ridder(double x1, double x2)
{double xacc=1.0e-9;
 return ridder(x1,x2,xacc);
}
default double ridder(double x1, double x2, double xacc)
{  
	int NR;
	int MAXITER=200;
	double ABC=-1.0e30;
	int j;
	double answer,f_high,f_low,f_mid,fnew,s,x_high,x_low,x_mid,xnew;
	f_low=func(x1);
	f_high=func(x2);
	if ((f_low > 0.0 && f_high < 0.0) || (f_low < 0.0 && f_high > 0.0)) {
		x_low=x1;
		x_high=x2;
		answer=ABC;
		for (j=0;j<MAXITER;j++) {
			x_mid=0.5*(x_low+x_high);
			f_mid=func(x_mid);
			s=Math.sqrt(f_mid*f_mid-f_low*f_high);
			if (s == 0.0) return answer;
			xnew=x_mid+(x_mid-x_low)*((f_low >= f_high ? 1.0 : -1.0)*f_mid/s);
			if (Math.abs(xnew-answer) <= xacc) return answer;
			answer=xnew;
			fnew=func(answer);
			if (fnew == 0.0) return answer;
			if (SIGN(f_mid,fnew) != f_mid) {
				x_low=x_mid;
				f_low=f_mid;
				x_high=answer;
				f_high=fnew;
			} else if (SIGN(f_low,fnew) != f_low) {
				x_high=answer;
				f_high=fnew;
			} else if (SIGN(f_high,fnew) != f_high) {
				x_low=answer;
				f_low=fnew;
			} else System.out.println("we should not reach to this point");
			if (Math.abs(x_high-x_low) <= xacc) return answer;
		}
		JOptionPane.showMessageDialog(null,"Warning : Maximum number of iteration is exceeded \n"+
        " result may not be valid","MAXİMUM ITERATION NUMBER WARNING",JOptionPane.WARNING_MESSAGE);
	}
	else {
		if (f_low == 0.0) return x1;
		if (f_high == 0.0) return x2;
		JOptionPane.showMessageDialog(null,"Warning : No roots existed between the given limits ","NO ROOTS WARNING",JOptionPane.WARNING_MESSAGE);
	}
	return 0.0;
}
//OPTIMISATION ONE DİMENSİONAL FUNCTION
default double golden(double a,double b,double epsilon)
{return golden(a,b,epsilon,1.0e-10,0);}
default double golden(f_x f,double a,double b)
{return golden(a,b,1.0e-10,1.0e-10,0);}

default double golden(double a,double b,double epsilon,double delta,int print)
{
//  find the minimum of the function
// note  maximum f(x) = minimum (-f(x))
    double r1 = (Math.sqrt(5.0)-1.0)/2.0; // golden ratio
    double r2 = r1*r1;
    double h = b - a;
    double ya = func(a);
    double yb = func(b);
    double c = a + r2*h;
    double d = a + r1*h;
    double yc = func(c);
    double yd = func(d);
    int k = 1;
    double dp,dy,p,yp;
    while ((Math.abs(yb-ya)>epsilon) || (h>delta))
    {
      k++;
      if (yc<yd)
        {
        b = d;
        yb = yd;
        d = c;
        yd = yc;
        h = b - a;
        c = a + r2 * h;
        yc = func(c);
        }
      else
       {
         a = c;
         ya = yc;
         c = d;
         yc = yd;
         h = b - a;
         d = a + r1 * h;
         yd = func(d);
       }//end of if
    }//end of while
    dp = Math.abs(b-a);
    dy = Math.abs(yb-ya);
    p = a;
    yp = ya;
    if (yb<ya)
    {
      p = b;
      yp = yb;
    }
    if(print==1)
    {System.out.println("x min = "+p+"ymin = "+yp+"error of x ="+dp+"error of y"+dy); }
    return p;
}
default double quadratic_search(double x0,double x2)
{double x1=(x0+x2)/2.0;
 return quadratic_search(x0,x1,x2);
}
default double quadratic_search(double x0,double x1,double x2)
{
double epsilon=1.0e-10;
double delta=1.0e-5;
int print=0;
return quadratic_search(x0,x1,x2,epsilon,delta,print);
}
default double quadratic_search(double x0,double x1,double x2,double epsilon,double delta,int print)
{  
int maxit=100;    
double f0 = func(x0);
    double f1 = func(x1);
    double f2 = func(x2);
    double f3 = f2;
    double x3 = 0;
    double h = x1 - x0;
    double k=1;
    double dd=Math.abs(f1-f0);
    while ((dd >epsilon) || (h>delta))
    {
      k++;
      x3=(f0*(x1*x1-x2*x2)+f1*(x2*x2-x0*x0)+f2*(x0*x0-x1*x1))/
             (2*f0*(x1-x2)+2.0*f1*(x2-x0)+2.0*f2*(x0-x1));
      f3=func(x3);      
      if(x3 >= x0 && x3<x1)
        {x2=x1;f2=f1;x1=x3;f1=f3;}
      else if(x3 >= x1 && x3 <x2) 
        {x0=x1;f0=f1;x1=x3;f1=f3; }
      else if(x3 >x2) 
        {x0=x1;f0=f1;x1=x2;f1=f2;x2=x3;f2=f3; }
      else if(x3 < x0) 
        {x0=x3;f0=f3;x1=x0;f1=f0;x2=x1;f2=f1; }
      dd=Math.abs(f1-f0);
      h=Math.abs(x1-x0);
      if(k>maxit) break;          
    }//end of while
    if(print==1)
    {System.out.println("x  = "+x3+"f = "+f3+"hata x ="+h+"hata f(x)=y "+dd); }
    return x3;
}
default double cubic_search_opt(double x0,double x1)
{
int nmax=500;
double tolerance=1.0e-15;
double f0,f1,df0,df1,fxs,dfxs;
double h=0.00001;
double xs;
int k=0;
double X,F,G,H,a,b,c;
  f0=func(x0);
  f1=func(x1);
  df0=dfunc(x0,h);
  df1=dfunc(x1,h);
while(Math.abs(x0-x1)>tolerance && k<nmax)
{ k++;	  
  X=1.0/(x1-x0);
  F=X*(f1-f0);
  G=X*(df1-df0);
  H=X*(F-df0);
  a=X*(G-2.0*H);
  b=3.0*H-G;
  c=df0;
  xs=x0+(Math.sqrt(b*b-3.0*a*c)-b)/(3.0*a);
  fxs=func(xs);
  dfxs=dfunc(xs,h);
  if(dfxs > 0 ) {x1=xs;f1=fxs;df1=dfxs;}
  else          x0=xs;f0=fxs;df0=dfxs; 
}
if(k>=nmax)  JOptionPane.showMessageDialog(null,"Maximum itaration number is exceeded \n"+
   " results may not be valid","MAKSIMUM ITERATION NUMBUR WARNING",JOptionPane.WARNING_MESSAGE);
return (x0+x1)/2.0;
}
default double newton_search(double x)
{
// root finding newton-raphson program

int nmax=500;
double tolerance=1.0e-15;
double fx,dfx;
for(int i=0;i<nmax;i++)
{
fx=dfunc(x);
dfx=dfunc(x,2);
x-=fx/dfx;
if(Math.abs(fx)<tolerance)
{
return x;
}
}
return x;
}
default double MAX(double a,double b) {return (a > b ? a : b);}

default double brent(double ax,double bx,double cx,double tol)
{
 double aa[]=brentf(ax,bx,cx,tol);
 return aa[0];
}

default double[] brentf(double ax,double bx,double cx,double tol)
{
// ax,bx,cx is three initial guesses (bx should be located between ax and cx)
// f : function (defined with f_x abstract class)
// tol tolerance	
int ITMAX=100;
double CGOLD=(3.0-Math.sqrt(5))/2.0; //golden ratio
double ZEPS=1.0e-10;
double xmin;
double aa[]=new double[2];
	int iter;
	double a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
	double e=0.0;
	d=0.0;
	a=((ax < cx) ? ax : cx);
	b=((ax > cx) ? ax : cx);
	x=w=v=bx;
	fw=fv=fx=func(x);
	for (iter=1;iter<=ITMAX;iter++) {
		xm=0.5*(a+b);
		tol2=2.0*(tol1=tol*Math.abs(x)+ZEPS);
		if (Math.abs(x-xm) <= (tol2-0.5*(b-a))) {
			xmin=x;
	        aa[0]=xmin;
	        aa[1]=fx;
	        return aa;
		}
		if (Math.abs(e) > tol1) {
			r=(x-w)*(fx-fv);
			q=(x-v)*(fx-fw);
			p=(x-v)*q-(x-w)*r;
			q=2.0*(q-r);
			if (q > 0.0) p = -p;
			q=Math.abs(q);
			etemp=e;
			e=d;
			if (Math.abs(p) >= Math.abs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
				d=CGOLD*(e=(x >= xm ? a-x : b-x));
			else {
				d=p/q;
				u=x+d;
				if (u-a < tol2 || b-u < tol2)
					d=SIGN(tol1,xm-x);
			}
		} else {
			d=CGOLD*(e=(x >= xm ? a-x : b-x));
		}
		u=(Math.abs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
		fu=func(u);
		if (fu <= fx) {
			if (u >= x) a=x; else b=x;
			{v=w;w=x;x=u;}
			{fv=fw;fw=fx;fx=fu;}
		} else {
			if (u < x) a=u; else b=u;
			if (fu <= fw || w == x) {
				v=w;
				w=u;
				fv=fw;
				fw=fu;
			} else if (fu <= fv || v == x || v == w) {
				v=u;
				fv=fu;
			}
		}
	}
	System.out.println("maximum number of iterations are exceeded in BRENT optimization");
	xmin=x;
	aa[0]=xmin;
	aa[1]=fx;
	return aa;
}
default double dbrent(double ax,double bx,double cx,double tol)
{
// ax,bx,cx is three initial guesses (bx should be located between ax and cx)
// f : function (defined with f_x abstract class)
// tol tolerance
int ITMAX=100;
double CGOLD=(3.0-Math.sqrt(5))/2.0; //altin oran
double ZEPS=1.0e-10;
	int iter;
	boolean ok1,ok2;
	double a,b,d,d1,d2,du,dv,dw,dx,e=0.0;
	d=0.0;
	double fu,fv,fw,fx,olde,tol1,tol2,u,u1,u2,v,w,x,xm;
	double xmin=bx;
	a=(ax < cx ? ax : cx);
	b=(ax > cx ? ax : cx);
	x=w=v=bx;
	fw=fv=fx=func(x);
	dw=dv=dx=dfunc(x);
	for (iter=1;iter<=ITMAX;iter++) {
		xm=0.5*(a+b);
		tol1=tol*Math.abs(x)+ZEPS;
		tol2=2.0*tol1;
		if (Math.abs(x-xm) <= (tol2-0.5*(b-a))) {
			xmin=x;
			return xmin;
		}
		if (Math.abs(e) > tol1) {
			d1=2.0*(b-a);
			d2=d1;
			if (dw != dx)  d1=(w-x)*dx/(dx-dw);
			if (dv != dx)  d2=(v-x)*dx/(dx-dv);
			u1=x+d1;
			u2=x+d2;
			ok1 = (a-u1)*(u1-b) > 0.0 && dx*d1 <= 0.0;
			ok2 = (a-u2)*(u2-b) > 0.0 && dx*d2 <= 0.0;
			olde=e;
			e=d;
			if (ok1 || ok2) {
				if (ok1 && ok2)
					d=(Math.abs(d1) < Math.abs(d2) ? d1 : d2);
				else if (ok1)
					d=d1;
				else
					d=d2;
				if (Math.abs(d) <= Math.abs(0.5*olde)) {
					u=x+d;
					if (u-a < tol2 || b-u < tol2)
						d=SIGN(tol1,xm-x);
				} else {
					d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
				}
			} else {
				d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
			}
		} else {
			d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
		}
		if (Math.abs(d) >= tol1) {
			u=x+d;
			fu=func(u);
		} else {
			u=x+SIGN(tol1,d);
			fu=func(u);
			if (fu > fx) {
				xmin=x;
				return xmin;
			}
		}
		du=dfunc(u);
		if (fu <= fx) {
			if (u >= x) a=x; else b=x;
			//MOV3(v,fv,dv, w,fw,dw)
			{v=w;fv=fw;dv=dw;}
			//MOV3(w,fw,dw, x,fx,dx)
			{w=x;fw=fx;dw=dx;}
			//MOV3(x,fx,dx, u,fu,du)
			{x=u;fx=fu;dx=du;}
		} else {
			if (u < x) a=u; else b=u;
			if (fu <= fw || w == x) {
				//MOV3(v,fv,dv, w,fw,dw)
				{v=w;fv=fw;dv=dw;}
				//MOV3(w,fw,dw, u,fu,du)
				{w=u;fw=fu;dw=du;}
			} else if (fu < fv || v == x || v == w) {
				//MOV3(v,fv,dv, u,fu,du)
				{v=u;fv=fu;dv=du;}
			}
		}
	}
	System.out.println("maximum number of iterations are exceeded in DBRENT optimization");
	return xmin;
}
}

