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

//Pade approximation
class pade1
{ 
public double q[];
public double p[];
public double a[];
double x0;
public double[] Maclaurin(if_x f,double x0i,int n)
{x0=x0i;
double a[]=new double[n+1];
double fact=1.0;
for(int i=0;i<=n;i++)
{a[i]=f.dfunc(x0,i)/fact;fact*=(i+1);}
return a;
}

public pade1(if_x f,int n,int m)
{
int N=m+n;
double a1[]=Maclaurin(f,0.0,N);
initilize(a1,n,m);
}

public pade1(double a[],int n,int m)
{initilize(a,n,m);}

public void initilize(double ai[],int n,int m)
{
int N=m+n;
int NP1=N+1;
a=new double[NP1];
for(int i=0;i<NP1;i++)
{a[i]=ai[i];}
q=new double[m+1];
p=new double[n+1];
q[0]=1.0;
p[0]=a[0];
int i,j;
double b[][]=new double[N+1][N+1];
double bb[][]=new double[N][N];
double c[]=new double[N];
for(i=1;i<=N;i++)
{ for(j=1;j<(i-1);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;}
  for(j=1;j<=i;j++)
  {if(j<=m)b[i][n+j]=-a[i-j];}
  for(j=n+i+1;j<=N;j++)
  {b[i][j]=0;}
  c[i-1]=a[i];
}
for(i=0;i<N;i++)
{for(j=0;j<N;j++)
 {bb[i][j]=b[i+1][j+1];}
}
double x[]=pivotlugauss(bb,c);
for(i=0;i<n;i++)
{p[i+1]=x[i];}
for(i=n;i<N;i++)
{q[i+1-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 double poly_func(double e[],double x)
{
// this function calculates the value of
// least square curve fitting function
int n=e.length;
double ff;
if(n!=0.0)
  { ff=e[n-1];
  for(int i=n-2;i>=0;i--)
   { ff=ff*(x-x0)+e[i]; }
  }
else
  ff=0;
return ff;
}

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

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

public double[][] result(if_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 void plot_func(if_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("Pad"+'\u00E9'+" rational approximation function");
	pp.plot();}

	public void plot_error(if_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("Pad"+'\u00E9'+" rational approximation error");
    pp.setYlabel("error");
	pp.plot();
}

public static void print(String s,String s1)
{JOptionPane.showMessageDialog(null,s,s1,JOptionPane.PLAIN_MESSAGE);}

public static void main(String arg[])
{ if_x f=(x)->Math.exp(x);
    int n=10;
	int m=0;
	pade1 Pa=new pade1(f,n,m);
	print(Pa.toString(),"Pad"+'\u00E9'+" rational polynomials coefficients");
	String s1[]={"x","Pad"+'\u00E9'+" f(x)=p(x)/q(x)","f(x)","error"};
	Text.print(Text.T(Pa.result(f,-1.0,1.0,20)),s1,"Pad"+'\u00E9'+" rational polynomials n = "+n+" m = "+m);
	Pa.plot_func(f,-1.0,1.0);
	Pa.plot_error(f,-1.0,1.0);
}

}
