break;
}
}
x=a[nn][nn];
if (l == nn) {
wri[nn]=new Complex(x+t,0);
a[nn][nn]=x+t;
nn--;
} else {
y=a[nn-1][nn-1];
w=a[nn][nn-1]*a[nn-1][nn];
if (l == nn-1) {
p=0.5*(y-x);
q=p*p+w;
z=sqrt(abs(q));
x += t;
a[nn][nn]=x;
a[nn-1][nn-1]=y+t;
if (q >= 0.0) {
z=p+SIGN(z,p);
wri[nn-1]=new Complex(x+z,0);
wri[nn]=new Complex(x+z,0);
if (z != 0.0) wri[nn]=new Complex(x-w/z,0);
x=a[nn][nn-1];
s=abs(x)+abs(z);
p=x/s;
q=z/s;
r=sqrt(p*p+q*q);
p /= r;
q /= r;
for (j=nn-1;j<n;j++) {
z=a[nn-1][j];
a[nn-1][j]=q*z+p*a[nn][j];
a[nn][j]=q*a[nn][j]-p*z;
}
for (i=0;i<=nn;i++) {
z=a[i][nn-1];
a[i][nn-1]=q*z+p*a[i][nn];
a[i][nn]=q*a[i][nn]-p*z;
}
for (i=0;i<n;i++) {
z=zz[i][nn-1];
zz[i][nn-1]=q*z+p*zz[i][nn];
zz[i][nn]=q*zz[i][nn]-p*z;
}
} else {
wri[nn]=new Complex(x+p,-z);
wri[nn-1]=wri[nn].conj();
}
nn -= 2;
} else {
if (its == 30) throw new IllegalArgumentException("Too many iterations in hqr");
if (its == 10 || its == 20) {
t += x;
for (i=0;i<nn+1;i++) a[i][i] -= x;
s=abs(a[nn][nn-1])+abs(a[nn-1][nn-2]);
y=x=0.75*s;
w = -0.4375*s*s;
}
++its;
for (m=nn-2;m>=l;m--) {
z=a[m][m];
r=x-z;
s=y-z;
p=(r*s-w)/a[m+1][m]+a[m][m+1];
q=a[m+1][m+1]-z-r-s;
r=a[m+2][m+1];
s=abs(p)+abs(q)+abs(r);
p /= s;
q /= s;
r /= s;
if (m == l) break;
u=abs(a[m][m-1])*(abs(q)+abs(r));
v=abs(p)*(abs(a[m-1][m-1])+abs(z)+abs(a[m+1][m+1]));
if (u <= EPS*v) break;
}
for (i=m;i<nn-1;i++) {
a[i+2][i]=0.0;
if (i != m) a[i+2][i-1]=0.0;
}
for (k=m;k<nn;k++) {
if (k != m) {
p=a[k][k-1];
q=a[k+1][k-1];
r=0.0;
if (k+1 != nn) r=a[k+2][k-1];
if ((x=abs(p)+abs(q)+abs(r)) != 0.0) {
p /= x;
q /= x;
r /= x;
}
}
if ((s=SIGN(sqrt(p*p+q*q+r*r),p)) != 0.0) {
if (k == m) {
if (l != m)
a[k][k-1] = -a[k][k-1];
} else
a[k][k-1] = -s*x;
p += s;
x=p/s;
y=q/s;
z=r/s;
q /= p;
r /= p;
for (j=k;j<n;j++) {
p=a[k][j]+q*a[k+1][j];
if (k+1 != nn) {
p += r*a[k+2][j];
a[k+2][j] -= p*z;
}
a[k+1][j] -= p*y;
a[k][j] -= p*x;
}
mmin = nn < k+3 ? nn : k+3;
for (i=0;i<mmin+1;i++) {
p=x*a[i][k]+y*a[i][k+1];
if (k+1 != nn) {
p += z*a[i][k+2];
a[i][k+2] -= p*r;
}
a[i][k+1] -= p*q;
a[i][k] -= p;
}
for (i=0; i<n; i++) {
p=x*zz[i][k]+y*zz[i][k+1];
if (k+1 != nn) {
p += z*zz[i][k+2];
zz[i][k+2] -= p*r;
}
zz[i][k+1] -= p*q;
zz[i][k] -= p;
}
}
}
}
}
} while (l+1 < nn);
}
if (anorm != 0.0) {
for (nn=n-1;nn>=0;nn--) {
p=wri[nn].re();
q=wri[nn].im();
na=nn-1;
if (q == 0.0) {
m=nn;
a[nn][nn]=1.0;
for (i=nn-1;i>=0;i--) {
w=a[i][i]-p;
r=0.0;
for (j=m;j<=nn;j++)
r += a[i][j]*a[j][nn];
if (wri[i].im() < 0.0) {
z=w;
s=r;
} else {
m=i;
if (wri[i].im() == 0.0) {
t=w;
if (t == 0.0)
t=EPS*anorm;
a[i][nn]=-r/t;
} else {
x=a[i][i+1];
y=a[i+1][i];
q=SQR(wri[i].re()-p)+SQR(wri[i].im());
t=(x*s-z*r)/q;
a[i][nn]=t;
if (abs(x) > abs(z))
a[i+1][nn]=(-r-w*t)/x;
else
a[i+1][nn]=(-s-y*t)/z;
}
t=abs(a[i][nn]);
if (EPS*t*t > 1)
for (j=i;j<=nn;j++)
a[j][nn] /= t;
}
}
} else if (q < 0.0) {
m=na;
if (abs(a[nn][na]) > abs(a[na][nn])) {
a[na][na]=q/a[nn][na];
a[na][nn]=-(a[nn][nn]-p)/a[nn][na];
} else {
Complex temp=new Complex(0.0,-a[na][nn]).div(new Complex(a[na][na]-p,q));
a[na][na]=temp.re();
a[na][nn]=temp.im();
}
a[nn][na]=0.0;
a[nn][nn]=1.0;
for (i=nn-2;i>=0;i--) {
w=a[i][i]-p;
ra=sa=0.0;
for (j=m;j<=nn;j++) {
ra += a[i][j]*a[j][na];
sa += a[i][j]*a[j][nn];
}
if (wri[i].im() < 0.0) {
z=w;
r=ra;
s=sa;
} else {
m=i;
if (wri[i].im() == 0.0) {
Complex temp = new Complex(-ra,-sa).div(new Complex(w,q));
a[i][na]=temp.re();
a[i][nn]=temp.im();
} else {
x=a[i][i+1];
y=a[i+1][i];
vr=SQR(wri[i].re()-p)+SQR(wri[i].im())-q*q;
vi=2.0*q*(wri[i].re()-p);
if (vr == 0.0 && vi == 0.0)
vr=EPS*anorm*(abs(w)+abs(q)+abs(x)+abs(y)+abs(z));
Complex temp= new Complex(x*r-z*ra+q*sa,x*s-z*sa-q*ra).div(new Complex(vr,vi));
a[i][na]=temp.re();
a[i][nn]=temp.im();
if (abs(x) > abs(z)+abs(q)) {
a[i+1][na]=(-ra-w*a[i][na]+q*a[i][nn])/x;
a[i+1][nn]=(-sa-w*a[i][nn]-q*a[i][na])/x;
} else {
temp= new Complex(-r-y*a[i][na],-s-y*a[i][nn]).div(new Complex(z,q));
a[i+1][na]=temp.re();
a[i+1][nn]=temp.im();
}
}
}
t=max(abs(a[i][na]),abs(a[i][nn]));
if (EPS*t*t > 1)