import static java.lang.Math.*;

abstract  class f_xxi extends Mathd
{ 
  // single function single independent variable
  // example f=x*x
  // Reference :"Generation of Finite Difference Formulas on Arbitrary Spaced Grids", 
  // Bength Fornberg, Mathematics of Computation, Volume 51, Number 184, October 1988
  // pages 699-706
double df;
double c[][][];
double alpha[];
double h;
int M;
public double a[];
//constructors
   public f_xxi(double ai[])
   {setA(ai);
   derivative_coef();
   df=0;
   }
   
   public f_xxi()
   { double ai[]=new double[0];
	 setA(ai);
	 derivative_coef();
	 df=0;
   }
   
  public void setA(double ai[])
  {int n=ai.length;
   a=new double[n];
   for(int i=0;i<n;i++)
   {a[i]=ai[i];}
  }
// abstract function func  
  abstract public double func(double x,int i);

// 1st order derivative 

public double dfunc(double x,int i)
{double h=1.0e-3;
 int n=1;
 return dfunc(x,i,n,h);
}
// nth order derivative  
public double dfunc(double x,int i,int n)
{double h=1.0e-2;
 return dfunc(x,i,n,h);
}

double dfunc(double x,int k,int N,double hi)
{// order of the maximum derivative
 // N order of derivative
  setH(hi);
  double deriv=0;
  double h1=1/h;
  double h2=1;
  for(int j=0;j<N;j++)
  {h2*=h1;}
  for(int i=0;i<c[0].length;i++)
  {deriv+=c[M][i][N]*func(x+alpha[i]*h,k);}
  return deriv*h2; 
}


public void setH(double hi)
{h=hi;}
public void coef_out(int N)
{  double a[][]=new double[alpha.length][2];
   for(int i=0;i<alpha.length;i++)
   {a[i][0]=alpha[i];
    a[i][1]=c[M][i][N];
   }
   String s1="Derivative coefficients for order of the highest derivative M = "+M;
   String s[]=new String[2]; 
   s[0]=""+'\u0394'+"x";
   if(N==1){s[1]="df(x)/dx";}
   else if(N==2) {s[1]="d"+'\u00B2'+"f(x)/dx"+'\u00B2';}
   else if(N==3) {s[1]="d"+'\u00B3'+"f(x)/dx"+'\u00B3';}
   else  {s[1]="d("+N+")f(x)/dx"+N;}
	Text.print(a,s,s1);
}
public void coef_out()
{  double a[][]=new double[alpha.length][M+1];
   String s[]=new String[M+1];
   s[0]=""+'\u0394'+"x";
   for(int j=1;j<=M;j++)
   {   
   if(j==1){s[1]="df(x)/dx";}
   else if(j==2) {s[j]="d"+'\u00B2'+"f(x)/dx"+'\u00B2';}
   else if(j==3) {s[j]="d"+'\u00B3'+"f(x)/dx"+'\u00B3';}
   else  {s[j]="d("+j+")f(x)/dx"+j;}
   }
   for(int i=0;i<alpha.length;i++)
   {a[i][0]=alpha[i];
    for(int j=1;j<=M;j++)
    a[i][j]=c[M][i][j];
   }
    String s1="Derivative coefficients for order of the highest derivative M = "+M;
	Text.print(a,s,s1);
}
public void setOrder(int Mi)
{M=Mi;}

public  void derivative_coef()
{
h=0.0001;
double x0=0;
M=60;
double alphai[]={0,1,-1,2,-2,3,-3,4,-4,5,-5,6,-6,7,-7,8,-8,9,-9,10,-10,11,-11,12,-12,13,-13,14,-14,15,
-15,16,-16,17,-17,18,-18,19,-19,20,-20,21,-21,22,-22,23,-23,24,-24,25,-25,26,-26,27,-27,28,-28,29,-29,30,-30,31,-31,32,-32,33,-33,34,-34,35,-35,36,-36,37,-37,38,-38,39,-39,30,-40,41,-41};
alpha=alphai;
int N=alpha.length-1;
// M   degree of highest derivative
// N+1 number of coefficients 
double delta[][][]=new double[N+1][N+1][M+1];
double c1,c2,c3;
delta[0][0][0]=1.0;
c1=1.0;
//for(int i=0;i<alpha.length;i++)
//{System.out.println(i+" "+alpha[i]);}
for(int n=1;n<=N;n++)
{  c2=1;
   for(int nu=0;nu<=(n-1);nu++)
   {c3=alpha[n]-alpha[nu];
    c2=c2*c3;
    if(n<=M) delta[n-1][nu][n]=0.0;
    for(int m=0;m<=Math.min(n,M);m++)
    {   //System.out.println("n="+n+"m="+m+"nu="+nu);
        if(m==0)
        { delta[n][nu][m]=((alpha[n]-x0)*delta[n-1][nu][m])/c3;}
        else
        {delta[n][nu][m]=((alpha[n]-x0)*delta[n-1][nu][m]-m*delta[n-1][nu][m-1])/c3;}
    }//next m
   }//next nu
   for(int m=0;m<=Math.min(n,M);m++)
   {   if(m==0)
	   {delta[n][n][m]=c1/c2*(-(alpha[n-1]-x0)*delta[n-1][n-1][m]);}
	   else
	   {delta[n][n][m]=c1/c2*(m*delta[n-1][n-1][m-1]-(alpha[n-1]-x0)*delta[n-1][n-1][m]);}   
   }//next m
   c1=c2;
}//next n
c=delta;
}

public double dfunc_n(double x,int i,int n,double h)
{
double hh=1/h;
double df=0;
if(n==0) df=func(x,i);
else if(n==1)
{df=(3.0*func(x-4.0*h,i)-32.0*func(x-3.0*h,i)+168.0*func(x-2.0*h,i)-672.0*func(x-h,i)+672.0*func(x+h,i)-168.0*func(x+2.0*h,i)+32.0*func(x+3.0*h,i)-3.0*func(x+4.0*h,i))/840.0*hh;}
else if(n==2)
{df=(-14350.0*func(x,i)-9.0*func(x-4*h,i)+128*func(x-3*h,i)-1008*func(x-2*h,i)+8064*func(x-h,i)+8064.0*func(x+h,i)-1008.0*func(x+2.0*h,i)+128.0*func(x+3.0*h,i)-9.0*func(x+4.0*h,i))/5040.0*hh*hh;}
else if(n==3)
{df=(-7.0*func(x-4.0*h,i)+72.0*func(x-3.0*h,i)-338.0*func(x-2.0*h,i)+488.0*func(x-h,i)-488.0*func(x+h,i)+338.0*func(x+2.0*h,i)-72.0*func(x+3.0*h,i)+7.0*func(x+4.0*h,i))/240.0*hh*hh*hh;}
else if(n==4)
{df=(2730.0*func(x,i)+7.0*func(x-4.0*h,i)-96.0*func(x-3.0*h,i)+676.0*func(x-2*h,i)-1952.0*func(x-h,i)-1952.0*func(x+h,i)+676.0*func(x+2.0*h,i)-96.0*func(x+3.0*h,i)+7.0*func(x+4.0*h,i))/240.0*hh*hh*hh*hh;}
else if(n==5)
{df=(func(x-4.0*h,i)-9.0*func(x-3.0*h,i)+26.0*func(x-2.0*h,i)-29.0*func(x-h,i)+29.0*func(x+h,i)-26.0*func(x+2.0*h,i)+9.0*func(x+3.0*h,i)-func(x+4.0*h,i))/6.0*hh*hh*hh*hh*hh;}
else if(n==6)
{df=(-150.0*func(x,i)-func(x-4.0*h,i)+12.0*func(x-3.0*h,i)-52.0*func(x-2.0*h,i)+116.0*func(x-h,i)+116.0*func(x+h,i)-52.0*func(x+2.0*h,i)+12.0*func(x+3.0*h,i)-func(x+4.0*h,i))/4.0*hh*hh*hh*hh*hh*hh;}
else if(n==7)
{df=(-func(x-4.0*h,i)+6.0*func(x-3.0*h,i)-14.0*func(x-2.0*h,i)+14.0*func(x-h,i)-14.0*func(x+h,i)+14.0*func(x+2.0*h,i)-6.0*func(x+3.0*h,i)+func(x+4.0*h,i))/2.0*hh*hh*hh*hh*hh*hh*hh;}
else if(n==8)
{df=(70.0*func(x,i)+func(x-4.0*h,i)-8.0*func(x-3.0*h,i)+28.0*func(x-2.0*h,i)-56.0*func(x-h,i)-56.0*func(x+h,i)+28.0*func(x+2.0*h,i)-8.0*func(x+3.0*h,i)+func(x+4.0*h,i))*hh*hh*hh*hh*hh*hh*hh*hh;}
else df=0;
return df;
}
 
}