(PC-12xx~14xx) Weierstrass / Durand-Kerner

Code:
' WEIERSTRASS / DURAND-KERNER POLYNOMIAL ROOT FINDER

' Estimate the complex roots of a polynomial

' 2nd order convergence, worst-case linear convergence for roots with multiplicities

'

' For SHARP PC-12xx to 14xx by Robert van Engelen

' To use with SHARP PC models with a single display line:

' - remove CURSOR and CLS

' - replace PRINT with PAUSE except at line 80

'

' Optimized implementations in C and BASIC:

' - Finds complex roots z of polynomial p such that |p(z)|<T given tolerance T

' - Prevents overflow/nan/inf by bounding the root shift adjustment magnitude during polishing

' - Prevents loss of precision when possible, e.g. in reciprocals

' - Tested with a large sample of polynomials

'

' See also: https://en.wikipedia.org/wiki/Durand–Kerner_method



' Algorithm to find complex double precision roots of a polynomial with complex coefficients:

'   int WDK(int n, const complex double a[], complex double z[], int maxit, double toler) {

'     int i, j, k, m;

'     complex double r, p, d;

'     double err, max, h, nprec, dprec;

'     if (maxit <= 0)

'       maxit = 100;

'     if (toler <= 0 || toler >= 1)

'       toler = 1e-8;

'     // optimize by removing trailing zero coefficients

'     for (m = n; m > 1 && a[m] == 0; --m)

'       z[m-1] = 0;

'     // leading coefficient cannot be zero

'     if (m < 1 || a[0] == 0)

'       return -1; // error

'     // generate initial roots

'     for (i = 1; i <= m; ++i)

'       z[i-1] = cpow(CMPLX(.4, .9), i-1);

'     // max root stride modulus permitted to prevent overflow

'     max = pow(DBL_MAX, 1./m);

'     if (max < 1)

'       max = 1;

'     // Weierstrass/Durand-Kerner root polishing

'     for (k = 1; k <= maxit; ++k) {

'       err = 0;

'       for (i = 1; i <= m; ++i) {

'         // compute p(z) using Horner's rule

'         p = a[0];

'         r = z[i-1];

'         for (j = 1; j <= m; ++j)

'           p = p * r + a[j];

'         // if p(z) = 0 within acceptable tolerance, then do not update

'         if (cabs(p) <= toler)

'           continue;

'         // compute denominator d multiplied by a[0] to force monic polynomial

'         d = a[0];

'         for (j = 1; j <= m; ++j)

'           if (i != j)

'             d *= r - z[j-1];

'         // compute p/d if d is not too small to prevent overflow

'         if (cabs(d) > toler)

'           p /= d;

'         h = cabs(p);

'         err += h;

'         // polish root if its shift adjustment is within floating point bounds to prevent overflow

'         if (h < max)

'           z[i-1] = r - p;

'       }

'       // check if convergence error is below our tolerance threshold

'       if (err <= toler)

'         break;

'     }

'     // round results (optional)

'     nprec = toler;

'     dprec = 1/toler;

'     for (i = 1; i <= m; ++i) {

'       r = z[i-1];

'       z[i-1] = nprec*CMPLX(round(creal(r)*dprec), round(cimag(r)*dprec));

'     }

'     return k; // number of iterations >0

'   }

'

' Algorithm to find complex roots of a polynomial with real coefficients, using separate real and imaginary variables:

'   int WDK(int n, const double a[], double zr[], double zi[], int maxit, double toler) {

'     int i, j, k, m;

'     double rr, pr, dr;

'     double ri, pi, di;

'     double s1, s2, k1, k2, k3;

'     double err2, max2, tol2; // err, max, toler squared

'     double h, t, nprec, dprec;

'     const double zh = sqrt(.97), zt = atan(2.25);

'     if (maxit <= 0)

'       maxit = 100;

'     if (toler <= 0 || toler >= 1)

'       toler = 1e-8;

'     tol2 = toler * toler;

'     // optimize by removing trailing zero coefficients

'     for (m = n; m > 1 && a[m] == 0; --m)

'       zr[m-1] = zi[m-1] = 0;

'     // leading coefficient cannot be zero

'     if (m < 1 || a[0] == 0)

'       return -1; // error

'     // generate initial roots z[i]=(.4+.9I)^i for i=0 to m-1

'     h = zh;

'     t = zt;

'     zr[0] = 1;

'     zi[0] = 0;

'     for (i = 2; i <= m; ++i) {

'       zr[i-1] = h * cos(t);

'       zi[i-1] = h * sin(t);

'       h *= zh;

'       t += zt;

'     }

'     // square of max root stride modulus permitted to prevent overflow

'     max2 = m > 2 ? pow(DBL_MAX, 2./m) : DBL_MAX;

'     if (max2 < 1)

'       max2 = 1;

'     // Weierstrass/Durand-Kerner root polishing

'     for (k = 1; k <= maxit; ++k) {

'       err2 = 0;

'       for (i = 1; i <= m; ++i) {

'         // compute p(z) using Horner's rule

'         pr = a[0];

'         pi = 0;

'         rr = zr[i-1];

'         ri = zi[i-1];

'         for (j = 1; j <= m; ++j) {

'           // compute p = p*r + a[j]

'           k1 = pr;

'           pr = pr * rr - pi * ri + a[j];

'           pi = pi * rr + k1 * ri;

'         }

'         // if p(z) = 0 within acceptable tolerance, then do not update

'         if (pr * pr + pi * pi <= tol2)

'           continue;

'         // compute denominator d multiplied by a[0] to force monic polynomial

'         dr = a[0];

'         di = 0;

'         for (j = 1; j <= m; ++j) {

'           if (i != j) {

'             // compute d = d * (r - z[j])

'             s1 = rr - zr[j-1];

'             s2 = ri - zi[j-1];

'             k1 = dr;

'             dr = dr * s1 - di * s2;

'             di = di * s1 + k1 * s2;

'           }

'         }

'         // compute p/d if d is not too small to prevent overflow

'         if (dr * dr + di * di > tol2) {

'           // compute p/d without loss of precision

'           if (fabs(dr) >= fabs(di)) {

'             k1 = di/dr;

'             k2 = k1*di + dr;

'             k3 = pr;

'             pr = (k1*pi + pr)/k2;

'             pi = (pi - k1*k3)/k2;

'           }

'           else {

'             k1 = dr/di;

'             k2 = k1*dr + di;

'             k3 = pr;

'             pr = (k1*pr + pi)/k2;

'             pi = (k1*pi - k3)/k2;

'           }

'         }

'         // polish root if its shift adjustment is within floating point bounds to prevent overflow

'         h = pr * pr + pi * pi;

'         err2 += h;

'         if (h < max2) {

'           zr[i-1] = rr - pr;

'           zi[i-1] = ri - pi;

'         }

'       }

'       // check if convergence error is below our tolerance threshold

'       if (err2 <= tol2)

'         break;

'     }

'     // round results (optional)

'     nprec = toler;

'     dprec = 1/toler;

'     for (i = 1; i <= m; ++i) {

'       zr[i-1] = nprec*round(zr[i-1]*dprec);

'       zi[i-1] = nprec*round(zi[i-1]*dprec);

'     }

'     return k; // number of iterations >0

'   }



' VARIABLES

'  N           degree

'  T           tolerance (1E-6)

'  L           max iterations (20), updated to the iterations executed

'  A(27..27+N) N+1 coefficients indexed from highest to lowest

'  A(O..O+N-1) N roots real parts

'  A(Z..Z+N-1) N roots imag parts

'  E           estimated error of the roots squared

'  O           array offset in A() with the real parts of the roots, O=28+N

'  Z           array offset in A() with the imag parts of the roots, Z=O+N

'  S           tolerance T squared

'  R           max root stride modulus squared

'  I,J,K       loop iterators

'  P,Q         real and imag parts of p(z)

'  A,B,C,D     scratch

'  U,V,W,X,Y   scratch



' driver program

10 "A" CLS: WAIT 0: INPUT "N=";N: IF N<2 GOTO 10

20 FOR I=0 TO N: CLS: PRINT "A";STR$(N-I);"=";A(27+I): CURSOR 24: INPUT A(27+I)

30 NEXT I

40 CURSOR: T=1E-6,L=20: INPUT "T=";T

50 INPUT "L=";L

60 CLS: PRINT "CALCULATING...": O=28+N: GOSUB "WDK"

70 CURSOR: PRINT "RE","IM": WAIT

80 FOR I=0 TO N-1: PRINT A(O+I),A(Z+I): NEXT I

90 END



' polynomial root finder: N degree, T tolerance, L max iters, A(27..27+N) coefficients

100 "WDK" S=T*T,O=28+N,Z=O+N,M=N

' optimize: remove leading zero coefficients (trailing A)

110 IF M>1 IF A(27+M)=0 LET M=M-1,A(O+M)=0,A(Z+M)=0: GOTO 110

120 IF M<1 OR A(27)=0 RETURN

' generate initial roots, square of max root shift distance bounded to prevent overflow

130 RADIAN: X=SQR .97,Y=ATN 2.25,U=X,V=Y,A(O)=1,A(Z)=0,R=1E99: IF M>2 LET R=R^(2/M): IF R<1 LET R=1

140 IF M>1 FOR I=1 TO M-1: A(O+I)=U*COS V,A(Z+I)=U*SIN V,U=U*X,V=V+Y: NEXT I

' iterative root polishing until convergence or max iterations reached

150 FOR K=1 TO L: E=0

160 FOR I=0 TO M-1

' compute p(z) using Horner's rule

170 A=A(O+I),B=A(Z+I),C=A(27),D=0,P=C,Q=0

180 FOR J=0 TO M-1: U=P,P=P*A-Q*B+A(28+J),Q=Q*A+U*B: NEXT J

' optimized: if |p(z)|<tolerance then do not update root z

190 IF P*P+Q*Q<=S GOTO 270

' compute denominator d multiplied by a[0] to force monic polynomial

200 FOR J=0 TO M-1

210 IF I<>J LET X=A-A(O+J),Y=B-A(Z+J),U=C,C=C*X-D*Y,D=D*X+U*Y

220 NEXT J

' compute p/d if d is not too small to prevent overflow

230 X=C*C+D*D: IF X<=S GOTO 260

240 IF ABS C>=ABS D LET U=D/C,V=U*D+C,W=P,P=(U*Q+P)/V,Q=(Q-U*W)/V: GOTO 260

250 U=C/D,V=U*C+D,W=P,P=(U*P+Q)/V,Q=(U*Q-W)/V

' polish root if its shift adjustment is within max bounds to prevent overflow

260 X=P*P+Q*Q,E=E+X: IF X<R LET A(O+I)=A-P,A(Z+I)=B-Q

270 NEXT I

280 CURSOR 24: PRINT "K=";STR$ K;":E=";SQR E: IF E<=S LET U=L,L=K,K=U

290 NEXT K

' round results (optional)

300 U=1/T

310 FOR I=0 TO M-1: A(O+I)=T*INT(U*A(O+I)+.5),A(Z+I)=T*INT(U*A(Z+I)+.5): NEXT I

' return error estimate

320 E=(E=0)*T+SQR E: RETURN