public class Leverrier
{
double D[][]; //inverse of matri
public static double trace(double B[][])
{
double sum=0;
int n=B.length;
for(int i=0;i<n;i++)
{sum+=B[i][i];
}
return sum;	
}

public static double[][] topla(double A[][],double a)
{
int n=A.length;
int m=A[0].length;
double B[][]=new double[n][m];
for(int i=0;i<n;i++)
{for(int j=0;j<m;j++)
  if(i!=j) B[i][j]=A[i][j];
  else     B[i][j]=A[i][j]+a;
}
return B;	
}
public static double[][] carp(double left,double[][] right)
{
//multiplying a matrix with a constant
int i,j;
int n=right.length;
int m=right[0].length;
double b[][];
b=new double[n][m];
for(i=0;i<n;i++)
  {
  for(j=0;j<m;j++)
    b[i][j]=right[i][j]*left;
  }
return b;
//end of multiplying a matrix with a constant double
}
public static double[][] carp(double[][] left,double[][] right)
{
//iki matrisin çarpımı
int ii,jj,i,j,k;
int m1=left[0].length;
int n1=left.length;
int m2=right[0].length;
int n2=right.length;
double[][] b;
b=new double[n1][m2];
 if(m1 != n2)
 {
 System.out.println("uyarı: matris boyutları aynı olmalıdır");
 for(ii=0;ii<n1;ii++)
   {
   for(jj=0;jj<m2;jj++)
     b[ii][jj]=0;
   }
   return b;
 }
 for(i=0;i<n1;i++)
   {
   for(j=0;j<m2;j++)
     {
     for(k=0;k<m1;k++){     
        //System.out.println(i+"  "+j+"   "+k);
	    b[i][j]+=left[i][k]*right[k][j];
	    //System.out.println(b[i][j]);
    }
     }
   }
return b;
}
public  double[] polynom_katsayi(double A[][])
{
int n=A.length;
int m=A[0].length;
double B[][]=new double[n][m];
for(int i=0;i<n;i++)
{for(int j=0;j<m;j++)
  B[i][j]=A[i][j];
}
System.out.println("i="+(n-1)+"\n"+Matrix.toString(B));
double a[]=new double[n+1];
double C[][]=new double[n][m];
a[n]=1.0;
a[n-1]=-trace(B);
for(int k=n-2;k>=0;k--)
{
C=topla(B,a[k+1]);	
B=carp(A,C);
System.out.println("i="+k+"\n"+Matrix.toString(C));
a[k]=-trace(B)/(double)(n-k);
}
D=carp((1.0/a[0]),C);
return a;
}
//polynom fonksiyonu
public static double func(double e[],double x)
{
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+e[i]; }
  }
else
  ff=0;
return ff;
}
public static double[][] findrootlimits(double e[],double x1,double x2,int n,int nbb)
{
	// find root limits
	// n : number of sub regions
	// x1,x2 : limit values
	// nbb = number of roots in the given region
	//  
	int nb;
	int i;
	double x,fp,fc,dx;
	double xb[][]=new double[2][nbb];
	nb=0;
	dx=(x2-x1)/n;
	x=x1;
	fp=func(e,x1);
	for (i=1;i<=n;i++) 
	{
		x+=dx;
		fc=func(e,x);
		// if the region with a root is found.....
		if (fc*fp < 0.0 || fp==0) {
			xb[0][nb]=x-dx;
			xb[1][nb]=x;
			nb++;
		}
		fp=fc;
		if (nbb == nb) return xb;
	}	
if( nb == 0)
	System.out.println("search is completed no roots is found");	
else if(nb<nbb)
  {	System.out.println("search is completed only "+nb+" roots are found \n"+
    "you search for"+nbb+" roots");	
   double xc[][]=new double[2][nb];
   for (i=0;i<nb;i++) {xc[0][i]=xb[0][i];xc[1][i]=xb[1][i];}
   return xc;
  }
return xb;
}
public static double brent(double e1[], double x1,double x2)
{
	//referans : Numerical Recipes in C, second edition, William H. Press, 
    //Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery
    //cambridge university press	
	double ITMAX=200;
    double EPS=3.0e-10;
	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(e1,a),fb=func(e1,b),fc,p,q,r,s,tol1,xm;
	if (fb*fa > 0.0) System.out.println("Root limits is not selected correctly \n results may not be correct");
	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(e1,b);
	}
	System.out.println("Warning : Maximum number of ,teration is exceeded \n"+
   " result might not be valid");
return b;
}
public  double[] Leverrier(double A[][])
{
int n=A.length;	
double a[]=polynom_katsayi(A);
System.out.println("coefficient of polynomial a[0]..a[n] \n"+Matrix.toStringT(a));
double lambda[]=new double[n];
lambda[n-1]=SCO13C.eigen_power(A);
lambda[0]=SCO13D.eigen_inversepower(A);
double xx[][]=findrootlimits(a,lambda[0]-0.1,lambda[n-1]+0.1,40,n);
System.out.println("root values : \n"+Matrix.toString(xx));
for(int k=1;k<(n-1);k++)
{lambda[k]=brent(a, xx[0][k],xx[1][k]);}
return lambda;
}

public static void main(String arg[])
{
//double A[][]={{1,2},{2,5}};
double A[][]={{12,6,-6},{6,16,2},{-6,2,16}};
Leverrier lv=new Leverrier();
double a[]=lv.Leverrier(A);
System.out.println("A-1 = \n"+Matrix.toString(lv.D));
System.out.println("lambda = \n"+Matrix.toStringT(a));
}
}
