import java.io.*;
import javax.swing.*;
class ff extends f_x
{
public double func(double x)
{return Math.exp(x);}
}

// Legendre  rational approximation
class Legendre
{ 
public double q[];
public double p[];
public double a[];

public static double P(int n, double x)
{
//Legendre polynomials
if(n==0) return 1.0;
else if(n==1) return x;
else
{return (P((n-1),x)*x*(double)(2*n-1)-P((n-2),x)*(double)(n-1))/(double)n;}
// end legendre polynomials
}
	    
    public  double fk(f_x f,int k,double x)
    {   double a=0;
        double mult=2.0/Math.PI;
	    a=mult*f.func(Math.cos(x))*Math.cos(k*x);
	    return a;
	}

public double adaptive_Lobatto_integral(f_x f,double a,double b,int k)
{   
	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]=fk(f,k,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(f,a,b,fa,fb,is,k);
  return Q;	
}
	
public  double Lobatto(f_x f,double a,double b,double fa,double fb,double is,int k)
{
  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]=fk(f,k,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şılamadı");
   }
   Q=i1;
  }  
  else
  {
    Q=Lobatto(f,a,mll,fa,fmll,is,k) + Lobatto(f,mll,ml,fmll,fml,is,k)+
      Lobatto(f,ml,m,fml,fm,is,k)   + Lobatto(f,m,mr,fm,fmr,is,k)+
      Lobatto(f,mr,mrr,fmr,fmrr,is,k)+Lobatto(f,mrr,b,fmrr,fb,is,k);
  }
  return Q;
}

//=============   


public Chebyshev(f_x f,int n,int m)
{
int N=m+n;
int Npm=N+m;
double ai[]=new double[Npm+1];
int k;
for(k=0;k<=(N+m);k++)
{ai[k]=adaptive_Lobatto_integral(f,0.0,Math.PI,k);  
}
initilize(ai,n,m);
}

public Chebyshev(double ai[],int n,int m)
{initilize(ai,n,m);}

public void initilize(double ai[],int n,int m)
{
int N=m+n;
int NP1=N+1;
int NPm=N+m;
int mp1=m+1;
int np1=n+1;
a=new double[NPm+1];
for(int i=0;i<=NPm;i++)
{a[i]=ai[i];}
q=new double[mp1];
p=new double[np1];
q[0]=1.0;
int i,j;
double b[][]=new double[NP1][NP1];
double c[]=new double[NP1];
for(i=0;i<=N;i++)
{  for(j=0;j<=i;j++)
   {if(j<=n) b[i][j]=0;}
if(i<=n) b[i][i]=1;
   for(j=i+1;j<=n;j++)
   {b[i][j]=0;}
   int k1,k2;
   for(j=n+1;j<=N;j++)
   {k1=i+j-n;
    k2=Math.abs(i-j+n);
   if(i!=0) b[i][j]=-0.5*(a[k1]+a[k2]);
   else     b[i][j]=-0.5*a[j-n];
  }
  if(i!=0) c[i]=a[i];
  else     c[i]=0.5*a[i];
}
double x[]=pivotlugauss(b,c);
for(i=0;i<=n;i++)
{p[i]=x[i];}
for(i=n+1;i<=N;i++)
{q[i-n]=x[i];}
}

  public static double[] pivotlugauss(double a[][],double b[])
  { //gauss elimination with partial pivoting
  int n=b.length;	  
  double x[]=new double[n];
  double carpan=0;
  double toplam=0;
  double buyuk;
  double dummy=0;
  //gauss elimination
  int i,j,k,p,ii,jj;
  for(k=0;k<(n-1);k++)
  { //pivoting
	p=k;
	buyuk=Math.abs(a[k][k]);
	for(ii=k+1;ii<n;ii++)
	{  dummy=Math.abs(a[ii][k]);
	   if(dummy > buyuk) {buyuk=dummy;p=ii;}
	}
	if(p!=k)
	{ for(jj=k;jj<n;jj++)
	  { dummy=a[p][jj];
	    a[p][jj]=a[k][jj];
	    a[k][jj]=dummy;
	  }
	  dummy=b[p];
	  b[p]=b[k];
	  b[k]=dummy;
	} 
    for(i=k+1;i<n;i++)
    {  carpan=a[i][k]/a[k][k];
       a[i][k]=0;
       for(j=k+1;j<n;j++)
       {   a[i][j]-=carpan*a[k][j]; }
           b[i]   =b[i]   -carpan*b[k];
	}
  }
  x[n-1]=b[n-1]/a[n-1][n-1];
  for(i=n-2;i>=0;i--)
  {
    toplam=0;
    for(j=i+1;j<n;j++)
    {  toplam+=a[i][j]*x[j];}
    x[i]=(b[i]-toplam)/a[i][i];
  }
  return x;
  }
  
public double func(double x)
{
return poly_func(p,x)/poly_func(q,x);
}
  
public static double poly_func(double a[],double x)
{
// this function calculates the value of
// least square curve fitting function
int n=a.length;
double b=0;
for(int k=0;k<n;k++)
{b+=a[k]*Ti(x,k);}
return b;
}

public void plot_func(f_x f,double a,double b)
{	double xx[]=new double[101];
    double yy[]=new double[101];
	for(int i=0;i<=100;i++)
	{xx[i]=a+i*(b-a)/100.0;yy[i]=func(xx[i]);}
	Plot pp=new Plot(xx,yy);
    pp.setPlabel("Chebyshev rational approximation function");
	pp.plot();}

	public void plot_error(f_x f,double a,double b)
    { double xx[]=new double[101];
      double yy[]=new double[101];
	for(int i=0;i<=100;i++)
	{xx[i]=a+i*(b-a)/100.0;yy[i]=f.func(xx[i])-func(xx[i]);}
	Plot pp=new Plot(xx,yy);
    pp.setPlabel("Chebyshev rational approximation error");
    pp.setYlabel("error");
	pp.plot();
}

public String toString()
{
	String s="p : \n"+Matrix.toStringT(p)+"\n";
	      s+="q : \n"+Matrix.toStringT(q)+"\n";
	return s;
}

public String toString(f_x f,double a,double b,int n)
{
	
return Matrix.toStringT(result(f,a,b,n));
}

public double[][] result(f_x f,double a,double b,int n)
{   double xx[][]=new double[4][n+1];
    for(int i=0;i<=n;i++)
	{ xx[0][i]=a+i*(b-a)/n;
	  xx[1][i]=func(xx[0][i]); 
	  xx[2][i]=f.func(xx[0][i]);
	  xx[3][i]=xx[2][i]-xx[1][i];
	}
	return xx;
}


public static void print(String s,String s1)
{JOptionPane.showMessageDialog(null,s,s1,JOptionPane.PLAIN_MESSAGE);}
public static void main(String arg[])
{
	ff f=new ff();
	int n=2;
	int m=2;
	Chebyshev Pa=new Chebyshev(f,n,m);
	print(Pa.toString(),"Chebyshev rational polynomials coefficients");
	String s1[]={"x","Chebyschev f(x)=p(x)/q(x)","f(x)","error"};
	Text.print(Text.T(Pa.result(f,-1.0,1.0,20)),s1,"Chebyshev rational polynomials n = "+n+" m = "+m);
	Pa.plot_func(f,-1.0,1.0);
	Pa.plot_error(f,-1.0,1.0);
}

}
