Post Reply 
(PC-12xx~14xx) Householder 3rd order method
03-10-2022, 01:24 AM
Post: #1
(PC-12xx~14xx) Householder 3rd order method
Another root finding method in my NumAlg collection. I've rewritten the algorithm to incorporate some new ideas to speed it up. This version does not include more advanced root starting point choices, such as suggested by Madsen, and stopping criteria, such as by Igarashi. However, this implementation appears to strike a good balance between computational cost (the speed of convergence) and accuracy.

HOUSEHOLDER POLYNOMIAL ROOT FINDER
Householder 3rd order method to find complex roots of a polynomial
4th order convergence, worst-case linear convergence for roots with multiplicities

Optimized implementations in C and BASIC:
- Finds complex roots z of polynomial p(z)
- Customizable floating point accuracy MACHEPS
- Uses a novel root starting point to improve performance
- Prevents overflow/nan/inf
- Prevents loss of precision when possible, e.g. in reciprocals
- Tested with a large sample of polynomials

See also: https://en.wikipedia.org/wiki/Householder%27s_method
W.H. Press et al. Numerical Recipes FORTRAN 2ed
David Goldberg. What Every Computer Scientist Should Know About Floating Point, ACM Computing Surveys, 1991
Madsen. A root-finding algorithm based on Newton Method, Bit 13 (1973) 71-75

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 60

1,774 bytes BASIC image (PC-1350)

Degree: max polynomial degree limited by time and memory only

Example:

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


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


(Solve x^4+2x^3-13x^2-14x+24=0)
[DEF-A]
N=4 A4=1 A3=2 A2=-13 A1=-14 A0=24
CALCULATING...
RE          IM
          3.          0.
         -4.          0.
         -2.          0.
          1.          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
CALCULATING...
RE          IM
       3000.          0.
      -4000.          0.
      -2000.          0.
          1.          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
CALCULATING...
RE          IM
          3.          4.
          3.         -4.
          2.          0.
         -1.          0.
          1.         -1.
          1.          1.


Code:
' HOUSEHOLDER POLYNOMIAL ROOT FINDER
' Householder 3rd order method to find complex roots of a polynomial
' 4th 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 60
'
' Optimized implementations in C and BASIC:
' - Finds complex roots z of polynomial p(z)
' - Customizable floating point accuracy MACHEPS
' - Uses a novel root starting point to improve performance
' - Prevents overflow/nan/inf
' - Prevents loss of precision when possible, e.g. in reciprocals
' - Tested with a large sample of polynomials

' See also: https://en.wikipedia.org/wiki/Householder%27s_method
' W.H. Press et al. Numerical Recipes FORTRAN 2ed
' David Goldberg. What Every Computer Scientist Should Know About Floating Point, ACM Computing Surveys, 1991
' Madsen. A root-finding algorithm based on Newton Method, Bit 13 (1973) 71-75

' Algorithm to find complex roots of a polynomial with complex coefficients:
'   // define MachEps machine epsilon precision, pow(FLT_RADIX, -DBL_MANT_DIG + 1) for double precision
'   #define MACHEPS 1E-16
'   // squared magnitude of complex z
'   #define cnorm(z) (creal(z) * creal(z) + cimag(z) * cimag(z))
'   void HH(int n, const complex double a[], complex double z[]) {
'     int i, m;
'     double e, t;
'     complex double *b, p, p1, p2, p3, r, h, u, v;
'     // find the trailing zero coefficients and remove them
'     for (m = n; m > 1 && a[m] == 0; --m)
'       z[m-1] = 0;
'     // leading coefficient cannot be zero
'     if (m < 1 || a[0] == 0)
'       abort();
'     // copy coefficients to deflate the polynomial
'     b = (complex double*)malloc((m + 1) * sizeof(complex double));
'     memcpy(b, a, (m + 1) * sizeof(complex double));
'     // loop to find the complex roots z[m-1] to z[2]
'     while (m > 2) {
'       // pick a suitable stating point for the root search
'       r = sqrt(0.5);
'       if (b[m-1] != 0)
'         r *= -b[m] / b[m-1];
'       // my novel idea (invention) to improve the complex root starting point for faster convergence
'       if (cimag(r) == 0)
'         r = CMPLX(creal(r), creal(r)/m);
'       // root modulus error tolerance
'       e = 6 * m * cabs(b[m]) * MACHEPS;
'       // loop to compute root z[m-1]
'       do {
'         // compute p(z), p'(z), p"(z), p'"(Z) using Horner's rule
'         p = b[0];
'         p1 = 0;
'         p2 = 0;
'         p3 = 0;
'         for (i = 1; i <= m; ++i) {
'           p3 = p3 * r + p2;
'           p2 = p2 * r + p1;
'           p1 = p1 * r + p;
'           p = p * r + b[i];
'         }
'         // if |p(z)| is (close to) zero then stop searching; using triangle inequality to prevent overflow
'         if (fabs(creal(p)) < e && fabs(cimag(p)) < e && cnorm(p) < e*e)
'           break;
'         // if stuck at a saddle point then randomize new starting point
'         if (fabs(creal(p1)) + fabs(cimag(p1)) < MACHEPS) {
'           r = -b[m] * CMPLX(1.0, 1.0/m) / (rand() + 1.0) * RAND_MAX;
'           continue;
'         }
'         h = -p/p1;
'         u = p2/p1;
'         v = p3/p1;
'         h *= (1 + h * u/2) / (1 + h * (u + h * v/6));
'         r += h;
'       }
'       while (cnorm(h) > (cnorm(r) + 1) * MACHEPS*MACHEPS);
'       if (fabs(cimag(r)) < MACHEPS) {
'         // |z.im| < MachEps, consider root is real
'         r = creal(r);
'       }
'       else {
'         // compute p(z.re) using Horner's rule
'         t = cnorm(p);
'         p = b[0];
'         for (i = 1; i <= m; ++i)
'           p = p * creal(r) + b[i];
'         // if |p(z.re)| <= |p(z)| then root is real
'         if (cnorm(p) <= t)
'           r = creal(r);
'       }
'       // save the complex root
'       z[m-1] = r;
'       // forward deflate the polynomial
'       v = b[0];
'       for (i = 1; i < m; ++i)
'         v = b[i] += v * r;
'       --m;
'     }
'     // compute the roots for the remaining m=1 and m=2 cases
'     if (m == 1) {
'       // linear case
'       z[0] = -b[1] / b[0];
'     }
'     else if (b[1] == 0) {
'       // quadratic case when b=0
'       z[0] = csqrt(-b[2] / b[0]);
'       z[1] = -z[0];
'     }
'     else {
'       // quadratic formula without potential catastrophic cancellation
'       v = b[1] * b[1] - 4 * b[0] * b[2];
'       if (creal(v) <= 0 || creal(b[1]) < 0)
'         z[0] = (-b[1] + csqrt(v)) / (2 * b[0]);
'       else
'         z[0] = (2 * b[2]) / (-b[1] - csqrt(v));
'       z[1] = b[2] / (b[0] * z[0]);
'     }
'     free(b);
'   }
'
' Algorithm to find complex roots of a polynomial with real coefficients, using separate real and imaginary variables:
'   // define MachEps machine epsilon precision, pow(FLT_RADIX, -DBL_MANT_DIG + 1) for double precision
'   #define MACHEPS 1E-16
'   // squared magnitude of complex (zr,zi)
'   #define norm(zr, zi) ((zr) * (zr) + (zi) * (zi))
'   void hh(int n, double a[], double zr[], double zi[]) {
'     int i, m;
'     double *b, e, k1, k2, pr, pi, pr1, pi1, pr2, pi2, pr3, pi3, rr, ri, hr, hi, ur, ui, vr, vi;
'     // find the trailing zero coefficients and remove them
'     for (m = n; m > 1 && a[m] == 0; --m)
'       zr[m-1] = zi[n-1] = 0;
'     // leading coefficient cannot be zero
'     if (m < 1 || a[0] == 0)
'       abort();
'     // copy coefficients to deflate the polynomial
'     b = (double*)malloc((m + 1) * sizeof(double));
'     memcpy(b, a, (m + 1) * sizeof(double));
'     // loop to find the complex roots z[m-1] to z[2]
'     while (m > 2) {
'       // pick a suitable stating point for the root search
'       rr = sqrt(0.5);
'       if (b[m-1] != 0)
'         rr *= -b[m] / b[m-1];
'       // my novel idea (invention) to improve the complex root starting point for faster convergence
'       ri = rr/m;
'       // root modulus error tolerance
'       e = 6 * m * fabs(b[m]) * MACHEPS;
'       // loop to compute root z[m-1]
'       do {
'         // compute p(z), p'(z), p"(z), p'"(Z) using Horner's rule
'         pr = b[0];
'         pi = 0;
'         pr1 = 0;
'         pi1 = 0;
'         pr2 = 0;
'         pi2 = 0;
'         pr3 = 0;
'         pi3 = 0;
'         for (i = 1; i <= m; ++i) {
'           // compute p3 = p3 * r + p2
'           k1 = pr3;
'           pr3 = pr3 * rr - pi3 * ri + pr2;
'           pi3 = pi3 * rr + k1  * ri + pi2;
'           // compute p2 = p2 * r + p1
'           k1 = pr2;
'           pr2 = pr2 * rr - pi2 * ri + pr1;
'           pi2 = pi2 * rr + k1  * ri + pi1;
'           // compute p1 = p1 * r + p
'           k1 = pr1;
'           pr1 = pr1 * rr - pi1 * ri + pr;
'           pi1 = pi1 * rr + k1  * ri + pi;
'           // compute p = p * r + b[i]
'           k1 = pr;
'           pr = pr * rr - pi * ri + b[i];
'           pi = pi * rr + k1 * ri;
'         }
'         // if |p(z)| is (close to) zero then stop searching; using triangle inequality to prevent overflow
'         if (fabs(pr) < e && fabs(pi) < e && norm(pr, pi) < e*e)
'           break;
'         // if stuck at a saddle point then randomize new starting point
'         if (fabs(pr1) + fabs(pi1) < MACHEPS) {
'           rr = -b[m] / (rand() + 1.0) * RAND_MAX;
'           ri = rr/m;
'           continue;
'         }
'         // compute h = -p/p1, u = p2/p1, v = p3/p1 without loss of precision
'         if (fabs(pr1) >= fabs(pi1)) {
'           k1 = pi1/pr1;
'           k2 = k1*pi1 + pr1;
'           hr = -(k1*pi + pr) / k2;
'           hi = (k1*pr - pi) / k2;
'           ur = (k1*pi2 + pr2) / k2;
'           ui = (pi2 - k1*pr2) / k2;
'           vr = (k1*pi3 + pr3) / k2;
'           vi = (pi3 - k1*pr3) / k2;
'         }
'         else {
'           k1 = pr1/pi1;
'           k2 = k1*pr1 + pi1;
'           hr = -(k1*pr + pi) / k2;
'           hi = (pr - k1*pi) / k2;
'           ur = (k1*pr2 + pi2) / k2;
'           ui = (k1*pi2 - pr2) / k2;
'           vr = (k1*pr3 + pi3) / k2;
'           vi = (k1*pi3 - pr3) / k2;
'         }
'         // compute denominator v' = h * (h * v/6 + u) + 1
'         k1 = (hr * vr - hi * vi) / 6 + ur;
'         k2 = (hi * vr + hr * vi) / 6 + ui;
'         vr = hr * k1 - hi * k2 + 1;
'         vi = hi * k1 + hr * k2;
'         // compute numerator u' = h * h * u/2 + h
'         k1 = (hr * ur - hi * ui) / 2;
'         k2 = (hi * ur + hr * ui) / 2;
'         ur = hr * k1 - hi * k2 + hr;
'         ui = hi * k1 + hr * k2 + hi;
'         // compute h = u'/v' = h * (1 + h * u/2) / (1 + h * (u + h * v/6))
'         if (fabs(vr) >= fabs(vi)) {
'           k1 = vi/vr;
'           k2 = k1*vi + vr;
'           hr = (k1*ui + ur) / k2;
'           hi = (ui - k1*ur) / k2;
'         }
'         else {
'           k1 = vr/vi;
'           k2 = k1*vr + vi;
'           hr = (k1*ur + ui) / k2;
'           hi = (k1*ui - ur) / k2;
'         }
'         rr += hr;
'         ri += hi;
'       }
'       while (norm(hr, hi) > (norm(rr, ri) + 1) * MACHEPS*MACHEPS);
'       // if |z.im| < MachEps then root is real
'       if (fabs(ri) < MACHEPS) {
'         ri = 0;
'       }
'       else {
'         // compute p(z.re) using Horner's rule
'         k1 = b[0];
'         for (i = 1; i <= m; ++i)
'           k1 = k1 * rr + b[i];
'         // if |p(z.re)| <= |p(z)| then root is real
'         if (k1*k1 <= norm(pr, pi))
'           ri = 0;
'       }
'       // save the complex root
'       zr[m-1] = rr;
'       zi[m-1] = ri;
'       if (ri == 0) {
'         // root is real; deflate the polynomial
'         k1 = b[0];
'         for (i = 1; i < m; ++i)
'           k1 = b[i] += k1 * rr;
'       }
'       else {
'         // root is complex; add conjugate root and deflate by quadratic
'         --m;
'         zr[m-1] = rr;
'         zi[m-1] = -ri;
'         k1 = -2 * rr;
'         k2 = rr*rr + ri*ri;
'         b[1] -= k1 * b[0];
'         for (i = 2; i < m; ++i)
'           b[i] -= k1 * b[i-1] + k2 * b[i-2];
'       }
'       --m;
'     }
'     // compute the roots for the remaining m=1 and m=2 cases
'     if (m == 1) {
'       // linear case
'       zr[0] = -b[1] / b[0];
'       zi[0] = 0;
'     }
'     else if (b[1] == 0) {
'       // quadratic case when b=0
'       k1 = -b[2] / b[0];
'       if (k1 > 0) {
'         zr[0] = sqrt(k1);
'         zi[0] = 0;
'         zr[1] = -zr[0];
'         zi[1] = 0;
'       }
'       else {
'         zr[0] = 0;
'         zi[0] = sqrt(-k1);
'         zr[1] = 0;
'         zi[1] = -zi[0];
'       }
'     }
'     else {
'       // quadratic formula without potential catastrophic cancellation
'       k1 = b[1] * b[1] - 4 * b[0] * b[2];
'       if (k1 >= 0) {
'         if (b[1] > 0) {
'           k2 = -b[1] - sqrt(k1);
'           zr[0] = 2 * b[2] / k2;
'           zi[0] = 0;
'           zr[1] = 0.5 * k2 / b[0];
'           zi[1] = 0;
'         }
'         else {
'           k2 = -b[1] + sqrt(k1);
'           zr[0] = 0.5 * k2 / b[0];
'           zi[0] = 0;
'           zr[1] = 2 * b[2] / k2;
'           zi[1] = 0;
'         }
'       }
'       else {
'         zr[0] = -0.5 * b[1] / b[0];
'         zi[0] = 0.5 * sqrt(-k1) / fabs(b[0]);
'         zr[1] = zr[0];
'         zi[1] = -zi[0];
'       }
'     }
'     free(b);
'   }

' VARIABLES
'  N              degree
'  A(27..27+N)    N+1 coefficients indexed from highest to lowest
'  A(28+N..28+2N) N+1 copied coefficients used in deflation
'  A(O..O+N-1)    N roots real parts
'  A(Z..Z+N-1)    N roots imag parts
'  O              array offset in A() with the real parts of the roots
'  Z              array offset in A() with the imag parts of the roots
'  I              loop iterator
'  L              array offset in A() with copied coefficients used in deflation
'  M              internal adjusted degree M<=N
'  A,B            real and imag parts of root z
'  P,Q            real and imag parts of p(z)
'  C,D            real and imag parts of p'(z)
'  F,G            real and imag parts of p"(z)
'  R,S            real and imag parts of p'"(z)
'  E              eps = 6*m*|a[m]|*MachEps
'  H,K            real and imag parts of h
'  U,X            real and imag parts of u
'  V,Y            real and imag parts of v
'  T,W            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 CLS: PRINT "CALCULATING...": GOSUB "HH"
50 PRINT "RE","IM": WAIT
60 FOR I=0 TO N-1: PRINT A(O+I),A(Z+I): NEXT I
70 END

' polynomial root finder: N degree, A(27..27+N) coefficients
100 "HH" L=28+N,O=L+N+1,Z=O+N,M=N
' optimized: 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
' copy coefficients to deflate the polynomial
130 FOR I=0 TO M: A(L+I)=A(27+I): NEXT I
' loop to find the complex roots z[m-1] to z[2]
140 IF M<3 GOTO 440
150 FOR M=M TO 3 STEP -1: CURSOR 24: PRINT "ROOT ";M
160 A=SQR .5: IF A(L+M-1)<>0 LET A=-A*A(L+M)/A(L+M-1)
170 B=A/M,E=6E-10*M*ABS A(L+M)
' loop to compute root z[m-1]
' compute p(z), p'(z), p"(z), p'"(Z) using Horner's rule
180 P=A(L),Q=0,C=0,D=0,F=0,G=0,R=0,S=0
190 FOR I=1 TO M
200 T=R,R=R*A-S*B+F,S=S*A+T*B+G,T=F,F=F*A-G*B+C,G=G*A+T*B+D
210 T=C,C=C*A-D*B+P,D=D*A+T*B+Q,T=P,P=P*A-Q*B+A(L+I),Q=Q*A+T*B
220 NEXT I
' if |p(z)| is (close to) zero then stop searching; using triangle inequality to prevent overflow
230 IF ABS P<E IF ABS Q<E IF P*P+Q*Q < E*E GOTO 360
' if stuck at a saddle point then randomize new starting point
240 IF ABS C+ABS D<1E-10 LET A=-A(L+M)/RND 1,B=A/M: GOTO 180
' compute h = -p/p1, u = p2/p1, v = p3/p1 without loss of precision
250 IF ABS C<ABS D GOTO 280
260 T=D/C,W=T*D+C,H=-(T*Q+P)/W,K=(T*P-Q)/W,U=(T*G+F)/W,X=(G-T*F)/W
270 V=(T*S+R)/W,Y=(S-T*R)/W: GOTO 300
280 T=C/D,W=T*C+D,H=-(T*P+Q)/W,K=(P-T*Q)/W,U=(T*F+G)/W,X=(T*G-F)/W
290 V=(T*R+S)/W,Y=(T*S-R)/W
' compute denominator v' = h * (h * v/6 + u) + 1
300 T=(H*V-K*Y)/6+U,W=(K*V+H*Y)/6+X,V=H*T-K*W+1,Y=K*T+H*W
' compute numerator u' = h * h * u/2 + h
310 T=(H*U-K*X)/2,W=(K*U+H*X)/2,U=H*T-K*W+H,X=K*T+H*W+K
' compute h = u'/v' = h * (1 + h * u/2) / (1 + h * (u + h * v/6))
320 IF ABS V>=ABS Y LET T=Y/V,W=T*Y+V,H=(T*X+U)/W,K=(X-T*U)/W: GOTO 340
330 T=V/Y,W=T*V+Y,H=(T*U+X)/W,K=(T*X-U)/W
340 A=A+H,B=B+K
' loop while |h|^2 > (|r|^2 + 1) * MachEps
350 IF H*H+K*K>(A*A+B*B+1)*1E-20 GOTO 180
' if |z.im| < MachEps then root is real
360 IF ABS B<1E-10 LET B=0: GOTO 390
' compute p(z.re) using Horner's rule
370 T=A(L): FOR I=1 TO M: T=T*A+A(L+I): NEXT I
' if |p(z.re)| <= |p(z)| then root is real
380 IF T*T<=P*P+Q*Q LET B=0
390 A(O+M-1)=A,A(Z+M-1)=B
' if root is real then deflate the polynomial
400 IF B=0 LET T=A(L): FOR I=1 TO M-1: T=A(L+I)+T*A,A(L+I)=T: NEXT I: GOTO 430
' root is complex; add conjugate root and deflate by quadratic
410 M=M-1:A(O+M-1)=A,A(Z+M-1)=-B,T=-2*A,W=A*A+B*B,A(L+1)=A(L+1)-T*A(L)
420 FOR I=2 TO M-1: A(L+I)=A(L+I)-T*A(L+I-1)-W*A(L+I-2): NEXT I
430 NEXT M
' compute the roots for the remaining m=1 and m=2 cases
440 IF M=1 LET A(O)=-A(L+1)/A(L),A(Z)=0: RETURN
450 IF A(L+1)<>0 GOTO 480
460 T=-A(L+2)/A(L): IF T>0 LET A(O)=SQR T,A(Z)=0,A(O+1)=-A(O),A(Z+1)=0: RETURN
470 A(O)=0,A(Z)=SQR -T,A(O+1)=0,A(Z+1)=-A(Z): RETURN
' quadratic formula without potential catastrophic cancellation
480 T=A(L+1),W=T*T-4*A(L)*A(L+2)
490 IF W<0 LET A(O)=-.5*T/A(L),A(Z)=.5*SQR(-W)/ABS A(L),A(O+1)=A(O),A(Z+1)=-A(Z): RETURN
500 IF T>0 LET W=-T-SQR W,A(O)=2*A(L+2)/W,A(Z)=0,A(O+1)=.5*W/A(L),A(Z+1)=0: RETURN
510 W=-T+SQR W,A(O)=.5*W/A(L),A(Z)=0,A(O+1)=2*A(L+2)/W,A(Z+1)=0: RETURN

"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
03-10-2022, 02:01 AM (This post was last modified: 03-10-2022 02:03 AM by robve.)
Post: #2
RE: (PC-12xx~14xx) Householder 3rd order method
Shown below is a comparison of Householder 3rd order method for polynomials with complex coefficients (algorithm "HH") and the method for polynomials with real coefficients (algorithm "hh"). For each of the four given MachEps values, the total number of iterations iter to converge the roots is shown, together with the sum of the residual absolute errors. Note that the last root or last two roots are computed directly, hence the polynomials shown in red font require no iterations. Note that the residual errors are expected to typically larger for polynomials of higher order with large coefficients, since a small change in the root z may cause a large change in p(z).

[Image: hh.png]

For the algorithms, I came up with the following idea to use a simple starting point choice for the complex root \( z \) of the (deflated) polynomial of order \( m \):
\( z.re = \frac{a_m}{\sqrt{2}\,a_{m-1}} \) if \( a_{m-1} \neq 0 \) otherwise \( z.re = \frac{1}{\sqrt{2}} \)
\( z.im = \frac{z.re}{m} \).
Note: \( a_m \) is the last coefficient, i.e. \( p(z)=...+a_{m-2}z^2+a_{m-1}z+a_m \)
Finding the roots in increasing order of magnitude is combined with forward deflation to improve numerical stability. This choice of starting point appears to perform very well, especially for polynomials with root multiplicities, but I have no analytic proof that explains why this is the case.

Choosing a root modulus error tolerance is key to converge to an acceptable accuracy without "overdoing" the iterations:
\( \epsilon = 6\,m\,|a_m|\,\mbox{MACHEPS} \).
The stopping criterion is simply:
\( |p(z)|<\epsilon \) and \( |h|\le(|r|+1)\,\mbox{MACHEPS} \).
This is more efficient to compute with the squared magnitudes (allowing a modest numerical difference in the formula):
\( |p(z)|^2<\epsilon^2 \) and \( |h|^2\le(|r|^2+1)\,\mbox{MACHEPS}^2 \)
We pick \( |r|^2+1 \) to bound the error to ~MACHEPS if \( |r|<1 \).

The algorithm checks if the root \( z \) is real-valued instead of complex if the imaginary part is negligible or if \( |p(z.re)|^2\le|p(z)|^2 \).

The "hh" version is optimized for polynomials with real coefficients. Therefore, two roots (the root and its conjugate) are found at a time if the root has a nonzero imaginary part.

- Rob

"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)