02-10-2022, 09:52 PM
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.
Edit: minor changes to reduce program size.
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: RETURNEdit: minor changes to reduce program size.