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| ****

: <root>    ( 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 <root>        \ 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]
