package stability;

/*
 * @author：fy  2004-11-25 15:35:51
 */
public class NiuDun {
/**
 * NiuDun 構造子注解。
 */
public NiuDun() {
	super();
}

//求實系數代數方程全部根的牛頓下山法
public static boolean ssrt(int n,double[] a ,double[] x ,double[] y){
	double[] s=new double[4];
	double x0,y0,x1,y1,p,q,p1,q1,t=1,c,d,u,dc;
	int i =0;
	boolean falg;
	x0=y0=x1=y1=p=q=p1=q1=0;
	if(a[0]== 0 )return false;
	if(n==1){
		x[0] = -a[1]/a[0];
		y[0]=0;
		return true;
		}
	if(n==2){
		p=a[1]*a[1]-4*a[0]*a[2];
		if(p>0){
			x[0]=(-a[1]+Math.sqrt(p))/2/a[0];
			x[1]=(-a[1]-Math.sqrt(p))/2/a[0];
			y[0]=y[1]=0;
			}
		else
		{
			x[0]=x[1]=-a[1]/2/a[0];
			y[0]=Math.sqrt(-p)/2/a[0];
			y[1]=-y[0];
			}
		return true;
	}
	
	while(Math.abs(a[n])<1e-12){
		x[n-1]=y[n-1]=0;
		n--;
		}
	u=Math.exp(Math.log(Math.abs(a[n]/a[0]))/n);
	p=u;
	for(i=n-1;i>=0;i--)
	   { a[i]*=p;  p*=u; }
	for(i=n;i>0;i--) a[i]/=a[0];
	a[0]=1;
	while(n>0){
		F(n,a,x0,y0,s);
		p=s[0]*s[0]+s[1]*s[1];
		q=s[2]*s[2]+s[3]*s[3];
		falg =true;
		while(true){
			if (p<1e-16 ) break;
			if (q>1e-16&& falg){
				t=1;
				while(true){
					x1=x0-t*(s[0]*s[2]+s[1]*s[3])/q;
					y1=y0-t*(s[1]*s[2]-s[0]*s[3])/q;
					F(n,a,x1,y1,s);
					p1=s[0]*s[0]+s[1]*s[1];
					q1=s[2]*s[2]+s[3]*s[3];
					if(p1<p){
						p=p1; q=q1; x0=x1;y0=y1;t=1;break;}
					t/=2;
					if(t<1e-10){ falg=false; break;}
					}
				
				}
			else
			{   d=0.01;
				while(true){
					c=0;
					dc=Math.PI/50;
					for(i=0;i<100;i++){
						x1=x0+d*Math.cos(c);
						y1=y0+d*Math.sin(c);
						c+=dc;
						if(x1*x1+y1*y1>1) continue;
						F(n,a,x1,y1,s);
						p1=s[0]*s[0]+s[1]*s[1];
				    	q1=s[2]*s[2]+s[3]*s[3];
				    	if(p1<p) break;
						}
					if(i==100) d+=0.01;
					else
					{
						p=p1; q=q1; t=1; d=0.1; x0=x1; y0=y1; falg=true; break;}
					}
				
				}
			}
		if(Math.abs(y0)<1e-7){
			x[n-1]=x1*u;
			y[n-1]=0;
			for(i=1;i<n;i++) a[i]+=x0*a[i-1];
			a[n]=0;
			n--;
			}
		else{
			x[n-1]=x[n-2]=x0*u;
			y[n-1]=y0*u;
			y[n-2]=-y[n-1];
			p=-2*x0; 
			q=x0*x0+y0*y0;
			for(i=1;i<n-1;i++){
				a[i]-=a[i-1]*p;
				a[i+1]-=a[i-1]*q;
				}
			a[n-1]=a[n]=0;
			n-=2;
			}
		x0=y0=0;
		}
	return true;
	}


   /*功能：使用秦九韶算法計算給定復數點處的函數值和導數值
     參數：n--整數，方程的次數
           a--雙精度實數數組，長度為(n+1)，記系數其中a[i]為(n-i)次方系數
           x--雙精度實數，要計算點的實部
           y--雙精度實數，要計算點的虛部
           s--雙精度實數數組，長度為4,s[0]、s[1]分別記函數值的實部和虛部。s[2]、s[3]記導數值實部和虛部。*/
   public static void F(int n,double[] a,double x,double y,double[] s){
	   int i;
	   double temp;
	   s[0]=a[0];
	   s[1]=0;
	   s[2]=a[0];
	   s[3]=0;
	   for(i=1;i<n;i++){
		   temp =s[0];
		   s[0]=x*temp-y*s[1]+a[i];
		   s[1]=x*s[1]+y*temp;
		   temp=s[2];
		   s[2]=x*temp-y*s[3]+s[0];
		   s[3]=x*temp+y*s[3]+s[1];
		   }
	   temp=s[0];
	   s[0]=x*temp-y*s[1]+a[n];
	   s[1]=x*s[1]+y*temp;
	   }

   public static void main(java.lang.String[] args) {
  
    double[] aa ={1,16,4,56};
    double[] r1=new double[3];
    double[] r2=new double[3];
    int i;
    ssrt(3,aa,r1,r2);

    for (i = 0; i < 3; i++) {
	    r1[i]=Math.round(r1[i]*1e8)/1e8;
	    r2[i]=Math.round(r2[i]*1e8)/1e8;
            if (r2[i] == 0)
                System.out.println( i+1+" "+r1[i]);
            else if(r2[i]<0) System.out.println(i+1+" "+ r1[i]+" "+ r2[i]+"i");
               else System.out.println(i+1+" "+r1[i]+"+"+r2[i]+"i");
        }
    }
}