FALSE [IF] Laguerre algorithm for polynomial roots ---------------------------------------------------- | (c) Copyright 2000 Julian V. Noble. | | Permission is granted by the author to | | use this software for any application pro- | | vided this copyright notice is preserved. | ---------------------------------------------------- This is an ANS Forth program requiring the FLOAT, FLOAT EXT, FILE and TOOLS EXT wordsets. Environmental dependences: Assumes PARSE can be used interpretively. Assumes "address units" in MOVE are bytes. (dpANS94: 6.1.1900) Assumes independent floating point stack Complex numbers reside on the fp stack as ( f: x y) where z = x + iy (Im above Re). The complex sqrt function, ZSQRT, is assumed to map (0, 2*pi) into (0, pi). That is, its branch cut is the positive real axis. Non-Standard words Most of these are in the files complex.f and arrays.f However as there is as yet no agreement as to the names of certain complex functions, I have chosen names that seemed sensible. Thus instead of ZABS I have defined |z| (more telegraphic). I also use the function |z|^2 which computes x^2 + y^2 = conjg(z) * z . The lexicon for arrays (arrays.f) builds both the number of elements and the data size into the header of an array. This information is required by the word }mov that moves data from one array to another array, since }mov is generic (that is, it works for any 2 arrays). }mov could easily be replaced by a definition that works only for complex arrays, should another array definition be preferred. Definitions using |z|, |z|^2 and }mov have been marked with **** for easy reference. Reference: F.S. Acton, "Numerical Methods that (Usually) Work" (Mathematical Ass'n of America, Washington, DC, 1990) Algorithm: For a given z, assumes z - z1 = a, and for all other roots, z - zn = b ; then G = p'(z)/p(z) = 1/a + (n-1)/b and H = G^2 - p"(z)/p(z) = 1/a^2 + (n-1)/b^2 . Eliminate b to get a = n/( G +- sqrt((nH-G^2)*(n-1)) ) . The next guess is z' = z - a . Iterate until converged, then deflate polynomial by the factor (z - root) and repeat. [THEN] MARKER -laguer \ define "[undefined]" if it does not already exist BL PARSE [undefined] DUP PAD C! PAD CHAR+ SWAP CHARS MOVE PAD FIND NIP 0= [IF] : [undefined] BL WORD FIND NIP 0= ; [THEN] include complex.f \ lexicon for complex arithmetic include arrays.f \ lexicon for arrays \ these are complex fp arrays 20 long 2 FLOATS 1array a{ \ (complex) coefficients of input polynomial 20 long 2 FLOATS 1array b{ \ coefficients of quotient polynomial 20 long 2 FLOATS 1array c{ \ coefficients of 1st derivative 20 long 2 FLOATS 1array d{ \ coefficients of 2nd derivative \ complex variables : zvariables ( n --) 0 DO CREATE 2 FLOATS ALLOT LOOP ; 3 zvariables p' zz zp fvariable epsilon 6 VALUE max_iter \ ------------------------------------------- synthetic division \ p[z] = (z-s) * q[z] + p[s] \ adr1 is address of coeff array of input polynomial p[z] \ adr2 is address of coeff array of quotient polynomial q[z] \ n is degree of polynomial : }zsynth ( adr1 adr2 n -- ) ( f: s -- p[s]) LOCALS| n chummy{ dummy{ | dummy{ n } z@ ( f: -- z sum) 0 n 1- DO \ count down from N to 1 zdup chummy{ I } z! \ b{ I } = sum zover z* dummy{ I } z@ z+ \ sum = sum * x + a{ I } -1 +LOOP ( f: s p[s]) znip ; FALSE [IF] Test case for synthetic division: : }. ( adr n -- ) \ display complex 1array 0 SWAP DO DUP I } z@ CR I . z. -1 +LOOP DROP ; 7.e0 0e0 a{ 0 } z! -5.e0 0e0 a{ 1 } z! 1.e0 0e0 a{ 2 } z! -14.e0 0e0 a{ 3 } z! 0.e0 0e0 a{ 4 } z! 3.e0 0e0 a{ 5 } z! a{ b{ 5 5.e0 0e0 )zsynth CR z. b{ 4 }. answers should be 7.63200E3 + i 0.00000E-1 ok 4 3.00000E0 + i 0.00000E-1 3 1.50000E1 + i 0.00000E-1 2 6.10000E1 + i 0.00000E-1 1 3.06000E2 + i 0.00000E-1 0 1.52500E3 + i 0.00000E-1 ok [THEN] \ --------------------------------------- end synthetic division : guessed? |z| epsilon F@ F< ; \ uses |z| **** : zmax ( f: z1 z2 -- z1 | z2) \ leave value with larger |z| zover zover ( f: z1 z2 z1 z2) |z|^2 f-rot |z|^2 ( f: -- z1 z2 |z2|^2 |z1|^2) F< IF zdrop ELSE znip THEN ; \ uses |z|^2 **** : new_diff ( n --) ( f: z -- a) LOCALS| n | zdup zz z! \ save initial guess a{ b{ n }zsynth ( f: -- p[z]) zdup guessed? IF EXIT THEN zz z@ b{ c{ n 1- }zsynth ( f: -- p[z] p'[z] ) p' z! zz z@ c{ d{ n 2 - }zsynth ( f: -- p[z] p"[z]/2 ) z2* zover z* znegate ( f: -- p[z] -p[z]*p"[z] ) p' z@ z^2 z+ ( f: -- p[z] H*p^2) n s>f z*f p' z@ z^2 z- ( f: -- p[z] n*H*p^2-p'^2) n 1- s>f z*f zsqrt ( f: -- p[z] R) zdup znegate ( f: -- p[z] R -R) p' z@ z+ ( f: -- p[z] R p'-R) zswap p' z@ z+ ( f: -- p[z] p'-R p'+R) zmax z/ ( f: -- p[z]/zmax[p'-R,p'+R]) n s>f z*f \ a = n*p/( p' +|- sqrt((nH*p^2-p'^2)*(n-1)) ) ; : }mov ( src dst n --) \ array_dst <-- array_src LOCALS| n dst{ src{ | dst{ @ \ get size of source data src{ @ \ get size of dst'n data <> ABORT" Inconsistent data types" 0 n DO src{ I } dst{ I } ( src[I] dst[I]) src{ @ ( #bytes) MOVE -1 +LOOP ; : new_z ( f: a --) znegate zz z@ z+ zp z! ; : apart? zz z@ zp z@ z- |z| epsilon f@ F> ; \ uses |z| **** : ( n --) ( f: -- root) 0 \ #iter = 0 LOCALS| #iter n | zz z@ n new_diff new_z \ compute zp = zz - a BEGIN apart? #iter max_iter < AND WHILE zp z@ zdup zz z! \ zz = zp n new_diff new_z \ zp = zz - a #iter 1+ TO #iter REPEAT ; : quadroots a{ 2 } z@ |z|^2 F0= \ uses |z|^2 **** IF a{ 1 } z@ |z|^2 F0= IF CR ." no roots!" ELSE CR ." only one root! " a{ 0 } z@ a{ 1 } z@ z/ znegate z. THEN ELSE a{ 1 } z@ znegate ( f: -b) zdup z^2 a{ 0 } z@ a{ 2 } z@ z* z2* z2* z- zsqrt ( f: -b d) zover zover ( f: -b d -b d) z+ z2/ a{ 2 } z@ z/ CR z. z- z2/ a{ 2 } z@ z/ CR z. THEN ; : roots ( n --) ( f: z0 epsilon --) LOCALS| n | epsilon F! zz z! BEGIN n 2 > WHILE n \ get a root CR zp z@ z. \ display it n 1- TO n \ n = n-1 b{ a{ n }mov \ create q_{n-1} [uses }mov]**** REPEAT quadroots ; FALSE [IF] Test case for roots: p(z) = z^6 + 4*z^5 - 6*z^4 - 4*z^3 - 7*z^2 - 48*z + 60 1.e0 0e0 a{ 6 } z! 4.e0 0e0 a{ 5 } z! -6.e0 0e0 a{ 4 } z! -4.e0 0e0 a{ 3 } z! -7.e0 0e0 a{ 2 } z! -48.e0 0e0 a{ 1 } z! 60.e0 0e0 a{ 0 } z! Say: 6 -10e0 0e0 1e-9 roots Should get -5.00000E0 + i 0.00000E-1 -2.00000E0 + i 0.00000E-1 7.09585E-11 + i 1.73205E0 1.00000E0 + i 6.14840E-12 2.00000E0 + i -3.36886E-12 2.93416E-11 + i -1.73205E0 ok [THEN]