Post Reply 
(PC-12xx~14xx) Weierstrass / Durand-Kerner
02-10-2022, 09:52 PM (This post was last modified: 03-10-2022 02:36 PM by robve.)
Post: #1
(PC-12xx~14xx) Weierstrass / Durand-Kerner
This is another one of my NumAlg programs sitting in my collection for years. I've added my tested and verified C version for comparison and commented the code. Perhaps someone finds this useful.

WEIERSTRASS / DURAND-KERNER POLYNOMIAL ROOT FINDER
Estimate the complex roots of a polynomial
Note: Weierstrass/Durand-Kerner is not guaranteed to converge
2nd order convergence, worst-case linear convergence for roots with multiplicities

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

For SHARP PC-12xx to 14xx

To use with SHARP PC models with a single display line:
- remove CURSOR and CLS
- replace PRINT with PAUSE except at line 80

954 bytes BASIC image (PC-1350)

Degree: max polynomial degree limited by time and memory only
Tolerance: T=1E-6 (adjustable)
Levels: up to L=20 iterations (adjustable)

Example:

(Solve x^2+2x+3=0)
[DEF-A]
N=2 A2=1 A1=2 A0=3
T= L=
CALCULATING...
RE          IM
         -1.  -1.414214.
         -1.   1.414214.

(Solve x^3-2x^2-x+2=0)
[DEF-A]
N=3 A3=1 A2=-2 A1=-1 A0=2
T= L=
CALCULATING...
RE          IM
          1.          0.
          2.          0.
         -1.          0.

(Solve x^3+3x^2+x+3=0)
[DEF-A]
N=3 A3=1 A2=3 A1=1 A0=3
T= L=
CALCULATING...
RE          IM
          0.         -1.
          0.          1.
         -3.          0.

(Solve x^4+2999x^3-10003e3x^2-2399e7x+24e9=0)
[DEF-A]
N=4 A4=1 A3=2999 A2=-10003E3 A1=-2399E7 A0=24E9
T= L=
CALCULATING...
RE          IM
          1.          0.
      -4000.          0.
       3000.          0.
      -2000.          0.

(Solve 5x^6-45x^5+225x^4-425x^3+170x^2+370x-500=0)
[DEF-A]
N=6 A6=5 A5=-45 A4=225 A3=-425 A2=170 A1=370 A0=-500
T= L=
CALCULATING...
RE          IM
          2.          0.
          1.          1.
          3.          4.
         -1.          0.
          3.         -4.
          1.         -1.


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

Edit: minor changes to reduce program size.

"I count on old friends" -- HP 71B,Prime|Ti VOY200,Nspire CXII CAS|Casio fx-CG50...|Sharp PC-G850,E500,2500,1500,14xx,13xx,12xx...
Visit this user's website Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)