314 lines
8.1 KiB
Plaintext
314 lines
8.1 KiB
Plaintext
pair sqrt(explicit pair z)
|
|
{
|
|
return sqrt(abs(z))*expi(angle(z,false)/2);
|
|
}
|
|
|
|
struct laguer {
|
|
pair x;
|
|
int its;
|
|
}
|
|
|
|
real EPS=10*realEpsilon;
|
|
laguer laguer(pair[] a, int m=a.length-1, pair x)
|
|
{/*
|
|
Given the degree 'm' and 'm+1' complex coefficients a[0...m] of polynomial sum(a[i]*x^i),
|
|
and given a complex value "x", this routine improves 'x' by Laguerre's method until it converges to a root of the given polynomial. The value of x found and the number of iterations taken is returned as 'laguer' structure.
|
|
Adapted from http://www.nrbook.com/a/bookcpdf/c9-5.pdf
|
|
*/
|
|
|
|
static real MR=8, MT=10, MAXIT=MT*MR;
|
|
int its;
|
|
laguer ol;
|
|
real[] ret;
|
|
int iter,j;
|
|
real abx,abp,abm,err;
|
|
pair dx,x1,b,d,f,g,h,sq,gp,gm,g2,ox=x;
|
|
static real[] frac={0, 0.5, 0.25, 0.75, 0.13, 0.38, 0.62, 0.88, 1};
|
|
|
|
laguer ret(){
|
|
ol.x=x;
|
|
ol.its=its;
|
|
return ol;
|
|
}
|
|
for (iter=1; iter <= MAXIT ; ++iter) {
|
|
its=iter;
|
|
b=a[m];
|
|
err=abs(b);
|
|
d=(0,0); f=d;
|
|
abx=abs(x);
|
|
for (int j=m-1; j >=0 ; --j) {
|
|
f=x*f+d;
|
|
d=x*d+b;
|
|
b=x*b+a[j];
|
|
err=abs(b)+abx*err;
|
|
}
|
|
err *= EPS;
|
|
if(abs(b) <= err) return ret();
|
|
g=d/b;
|
|
g2=g^2;
|
|
h=g2-2*f/b;
|
|
sq=sqrt((m-1)*(m*h-g2));
|
|
gp=g+sq;
|
|
gm=g-sq;
|
|
abp=abs(gp);
|
|
abm=abs(gm);
|
|
if(abp < abm) gp=gm;
|
|
dx=(max(abp,abm) > 0.0) ? (m,0)/gp :
|
|
(1+abx)*(cos(iter),sin(iter));
|
|
x1=x-dx;
|
|
if(x.x == x1.x && x.y == x1.y) return ret();
|
|
if(iter%MT == 0) x=x1;
|
|
else x=x-frac[floor(iter/MT)]*dx;
|
|
}
|
|
abort("Too many iteration in laguer.");
|
|
return ret();
|
|
}
|
|
|
|
pair[] zroots(pair[] a, bool polish=true)
|
|
{/*
|
|
Given the complex coefficients a[0...m] of the polynomial sum(a[i]*x^i),
|
|
this routine returns all the complex roots by Laguer's method.
|
|
If 'polish' is true, the roots are polished (also by Laguer's method');
|
|
Adapted from http://www.nrbook.com/a/bookcpdf/c9-5.pdf
|
|
*/
|
|
laguer L;
|
|
int m=a.length-1,its=L.its;
|
|
pair[] ad=copy(a);
|
|
pair x=L.x,b,c;
|
|
pair[] roots;
|
|
for (int j=m; j >= 1 ; --j) {
|
|
x=(0,0);
|
|
L=laguer(ad,j,L.x);
|
|
x=L.x; its=L.its;
|
|
if(abs(x.y) <= 2*EPS*abs(x.x)) x=(x.x,0);
|
|
roots[j-1]=x;
|
|
b=ad[j];
|
|
for (int k=j-1; k >= 0 ; --k) {
|
|
c=ad[k];
|
|
ad[k]=b;
|
|
b=x*b+c;
|
|
}
|
|
}
|
|
|
|
if(polish) {
|
|
for (int j=0; j < m ; ++j) {
|
|
L.x=roots[j];
|
|
L.its=its;
|
|
L=laguer(a,L.x);
|
|
roots[j]=L.x;
|
|
its=L.its;
|
|
}
|
|
}
|
|
for (int j=0; j < m ; ++j) {
|
|
pair x=roots[j];
|
|
int i;
|
|
for (i=j-1; i >= 0 ; --i) {
|
|
if(roots[i].y == 0) break;
|
|
roots[i+1]=roots[i];
|
|
}
|
|
roots[i+1]=x;
|
|
}
|
|
return roots;
|
|
}
|
|
|
|
real sgn(real a, real b)
|
|
{
|
|
return (b == 0) ? abs(a): sgn(b)*abs(a);
|
|
}
|
|
|
|
void hqr(real[][] a, int n, real[] wr=new real[]{}, real[] wi= new real[]{})
|
|
{/*
|
|
Finds all eigenvalues of an upper Hessenberg matrix a[1...n][1..n].
|
|
The real and imaginary parts of eigenvalues are returned in wr[1...n] and wi[1...n], respectively.
|
|
Adapted from http://www.nrbook.com/a/bookcpdf/c11-6.pdf
|
|
*/
|
|
int nn,m,l,k,j,its,i,mmin;
|
|
real z,y,x,w,v,u,t,s,r,q,p,anorm=0;
|
|
for (i=1; i <= n ; ++i)
|
|
for (j=max(i-1,1); j <= n ; ++j) anorm += abs(a[i][j]);
|
|
nn=n;t=0;
|
|
while(nn >= 1) {
|
|
its=0;
|
|
do {
|
|
for (l=nn; l >= 2 ; --l) {
|
|
s=abs(a[l-1][l-1])+abs(a[l][l]);
|
|
if(s == 0) s=anorm;
|
|
if(abs(a[l][l-1])+s == s) {
|
|
a[l][l-1]=0;
|
|
break;
|
|
}
|
|
}
|
|
x=a[nn][nn];
|
|
if(l == nn) {
|
|
wr[nn]=x+t;
|
|
wi[--nn]=0;
|
|
} 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^2+w;
|
|
z=sqrt(abs(q));
|
|
x += t;
|
|
if(q >= 0) {
|
|
z=p+sgn(z,p);
|
|
wr[nn-1]=wr[nn]=x+z;
|
|
if(z != 0) wr[nn]=x-w/z;
|
|
wi[nn-1]=wi[nn]=0;
|
|
} else {
|
|
wr[nn-1]=wr[nn]=x+p;
|
|
wi[nn-1]=-(wi[nn]=z);
|
|
}
|
|
nn -= 2;
|
|
} else {
|
|
if(its == 30) abort("Too many iterations in hqr...");
|
|
if(its == 10 || its == 20) {
|
|
t += x;
|
|
for(i=1; i <= nn; ++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^2;
|
|
}
|
|
++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 == 1) 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+v == v) break;
|
|
}
|
|
for (i=m+2; i <= nn; ++i) {
|
|
a[i][i-2]=0;
|
|
if(i != m+2) a[i][i-3]=0;
|
|
}
|
|
for (k=m; k <= nn-1; ++k) {
|
|
if(k != m) {
|
|
p=a[k][k-1];
|
|
q=a[k+1][k-1];
|
|
r=0;
|
|
if(k != nn-1) r=a[k+2][k-1];
|
|
if((x=abs(p)+abs(q)+abs(r)) != 0) {
|
|
p /= x;
|
|
q /= x;
|
|
r /=x;
|
|
}
|
|
}
|
|
if((s=sgn(sqrt(p^2+q^2+r^2),p)) != 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 <= nn; ++j) {
|
|
p=a[k][j]+q*a[k+1][j];
|
|
if(k != nn-1) {
|
|
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=l; i <= mmin; ++i) {
|
|
p=x*a[i][k]+y*a[i][k+1];
|
|
if(k != nn-1) {
|
|
p += z*a[i][k+2];
|
|
a[i][k+2] -= p*r;
|
|
}
|
|
a[i][k+1] -= p*q;
|
|
a[i][k] -= p;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
} while(1 < nn-1);
|
|
}
|
|
}
|
|
|
|
real RADIX=8;
|
|
void balanc(real[][] a, int n=a[0].length)
|
|
{/* Given a matrix 'a[1...n][1..n]', this routine replaces it by a balanced matrix with identical eigenvalues.
|
|
A symmetric matrix is already balanced and is unaffected by this procedure.
|
|
The variable RADIX should be the machine's floating-point radix.
|
|
Adapted from http://www.nrbook.com/a/bookcpdf/c11-5.pdf
|
|
*/
|
|
int last=0, j, i;
|
|
real s, r, g, f, c, sqrdx=RADIX^2;
|
|
while(last == 0) {
|
|
last=1;
|
|
for (i=1; i <= n; ++i) {
|
|
r=c=0;
|
|
for (j=1; j <= n; ++j)
|
|
if(j != i) {
|
|
c += abs(a[j][i]);
|
|
r += abs(a[i][j]);
|
|
}
|
|
if(c != 0 && r != 0) {
|
|
g=r/RADIX;
|
|
f=1;
|
|
s=c+r;
|
|
while(c < g) {
|
|
f *= RADIX;
|
|
c *= sqrdx;
|
|
}
|
|
g=r*RADIX;
|
|
while(c > g) {
|
|
f /= RADIX;
|
|
c /= sqrdx;
|
|
}
|
|
if((c+r)/f < 0.95*s) {
|
|
last=0;
|
|
g=1/f;
|
|
for (j=1; j <= n; ++j) a[i][j] *= g;
|
|
for (j=1; j <= n; ++j) a[j][i] *= j;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
pair[] zroots(real[] a, bool polish=true)
|
|
{/* Return all the roots of a polynomial with real coefficients, sum(a[i]x^i).
|
|
The method is to construct an upper Hessenberg matrix whose eigenvalues are
|
|
the desired roots, and the use the routine 'balanc' and 'hqr'.
|
|
If polish is true, root-polishing by Newton-Raphson's method is applied.
|
|
Adapted from http://www.nrbook.com/a/bookcpdf/c9-5.pdf
|
|
*/
|
|
int m=a.length-1;
|
|
int j,k;
|
|
real[][] hess=new real[m+1][m+1];
|
|
real[] rtr=sequence(new real(int n){return 0;},m+1),
|
|
rti=sequence(new real(int n){return 0;},m+1);
|
|
real xr,xi;
|
|
if(a[m] == 0) abort("Bad args in zroots.");
|
|
for (k=1; k <= m; ++k) {
|
|
hess[1][k]=-a[m-k]/a[m];
|
|
for(j=2; j <= m; ++j) hess[j][k]=0;
|
|
if(k != m) hess[k+1][k]=1;
|
|
}
|
|
balanc(hess,m);
|
|
hqr(hess,m,rtr,rti);
|
|
for (j=2; j <= m; ++j) {
|
|
xr=rtr[j];
|
|
xi=rti[j];
|
|
for (k=j-1; k >= 1; --k) {
|
|
if(rtr[k] <= xr) break;
|
|
rtr[k+1]=rtr[k];
|
|
rti[k+1]=rti[k];
|
|
}
|
|
rtr[k+1]=xr;
|
|
rti[k+1]=xi;
|
|
}
|
|
pair[] roots=sequence(new pair(int n){return (rtr[n],rti[n]);},m+1);
|
|
roots.delete(0);
|
|
return roots;
|
|
}
|