FALSE [IF] Laguerre algorithm for polynomial roots ---------------------------------------------------- | (c) Copyright 2001 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 include ftran201.f \ FORmula TRANslator 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 ; 7 zvariables p pp ppp zz zp R zsum 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: z -- p[z]) LOCALS| n out{ in{ | zz z! \ save guess zz" zsum = in{ n }" 0 n 1- DO \ count down from N to 1 zz" out{ I } = zsum" zz" zsum = zsum * zz + in{ I }" -1 +LOOP ( f: p[z]) zz" zsum" ; \ --------------------------------------- end synthetic division : zabs |z| ; \ aliases : zabsq |z|^2 ; : is_root? |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 **** : zfloat ( n -- ) ( f: -- n 0) \ use as zfloat(0) fdrop S>F FSWAP ; : zdiff ( n --) ( f: z -- a) LOCALS| n | zz z! \ save initial guess a{ b{ n zz" p = zsynth(zz) " \ compute p(z) zz" p" is_root? \ done? IF zz" p" EXIT THEN \ a=p(z) b{ c{ n 1- zz" pp = zsynth(zz) " \ compute p'(z) c{ d{ n 2 - zz" ppp = zsynth(zz) * 2 " \ compute p"(z) n 1- n ( n-1 n) zz" R = zsqrt( (zfloat(0) * (pp^2 - p*ppp) - pp^2 ) * zfloat(0) )" n zz" zfloat(0)* p /zmax(pp-R,pp+R)" ; : }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 ; : apart? zz" zabs( zz - zp )" f" epsilon" F> ; : ( n --) ( f: -- root) 0 LOCALS| #iter n | \ #iter = 0 n zz" zp = zz - zdiff(zz)" \ zp = zz - a BEGIN apart? #iter max_iter < AND WHILE zz" zz = zp " n zz" zp = zz - zdiff(zz)" #iter 1+ TO #iter REPEAT ; : quadroots zz" zabsq( a{ 2 } )" F0= IF zz" zabsq( a{ 1 } )" F0= IF CR ." no more roots!" \ a{1} and a{2} = 0 ELSE CR ." single root! " \ a{2} = 0 zz" -a{ 0 } / a{ 1 } " z. THEN ELSE \ a{2} <> 0 hence 2 roots zz" zz = -a{ 1 }/(2*a{ 2 })" zz" zp = zsqrt( zz^2 - a{ 0 }/a{ 2 } )" zz" zz + zp" CR z. zz" zz - zp" 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 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 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]