(   Title:  Noble Laguerre algorithm tests
   Author:  David N. Williams
     File:  lagroots-test.fs
  License:  Public Domain
  Version:  1.0.0
  Revised:  November 30, 2010
)

\ debugging
[UNDEFINED] \\ [IF]
: \\  ( -- )  -1 parse 2drop BEGIN refill 0= UNTIL ;
[THEN]
: .#  ( n -- )  cr ." #" . cr ;

s" complex.fs"        INCLUDED
s" fsl-util.fs"       INCLUDED
s" lagroots.fs"       INCLUDED
s" ttester-xf.fs"     INCLUDED
s" tester-display.fs" INCLUDED

decimal
5 SET-PRECISION
true VERBOSE !
?.cr  \ before the first displayed line


\ ***BEGIN ERROR DISPLAY

VARIABLE XT-#ERRORS    0 XT-#ERRORS !

: XT-ERROR1  ( c-addr u -- )
\ Display an error message followed by the line that had the
\ error.
  red-text 1 XT-#ERRORS +! XT-ERROR-DEFAULT normal-text ;
\ ' XT-ERROR1 XT-ERROR-XT !

: xt-a.  ( &array c-addr u -- )
  blue-text type
  ( &array) dup >r @ dup
  IF   1- r> cell+ swap cells over +
       DO i @ . -1 cells +LOOP
  ELSE r> 2drop ."  none" THEN cr normal-text ;

: XT-ERROR2  ( c-addr u -- )
  XT-ERROR1
  XT-RESULTS s"   xresults: " xt-a.
  XT-GOALS   s"     xgoals: " xt-a. ;
' XT-ERROR2 XT-ERROR-XT !

VARIABLE FT-#ERRORS    0 FT-#ERRORS !

: FT-ERROR1  ( c-addr u -- )
\ Display an error message followed by the line that had the
\ error.
  red-text 1 FT-#ERRORS +! FT-ERROR-DEFAULT normal-text ;
\ ' FT-ERROR1 FT-ERROR-XT !

: ft-a.  ( &array c-addr u -- )
  blue-text type
  ( &array) dup >r @ dup
  IF   1- r> cell+ faligned swap floats over +
       DO i f@ fs. -1 floats +LOOP
  ELSE r> 2drop ."  none" THEN cr normal-text ;

: FT-ERROR2  ( c-addr u -- )
  FT-ERROR1
  FT-RESULTS s"   fresults: " ft-a.
  FT-GOALS   s"     fgoals: " ft-a. ;
' FT-ERROR2 FT-ERROR-XT !

: zt-a.  ( &array c-addr u -- )
  blue-text type cr
  ( &array) dup >r @ dup
  IF   2 - r> cell+ faligned swap floats over +
       DO i z@ fswap zs. cr -2 floats +LOOP
  ELSE r> 2drop ."  none" THEN normal-text ;

: FT-ERROR3  ( c-addr u -- )
  FT-ERROR1
  FT-RESULTS s"   zresults" zt-a.
  FT-GOALS   s"   zgoals"   zt-a. cr ;

: ?.errors  ( -- )
  VERBOSE @
  IF blue-text ." XT-ERRORS " normal-text XT-#ERRORS @ .
     blue-text ." FT-ERRORS " normal-text FT-#ERRORS @ .
  THEN ;

\ ***END ERROR DISPLAY


\ Assume fp results and goals are complex numbers.
: z{  FT-ERROR-XT @  ['] FT-ERROR3 FT-ERROR-XT ! t{ ;
: }z  }t FT-ERROR-XT ! ;

PFE-HOST [IF]
  ?." Host is pfe
[THEN]
GFORTH-HOST [IF]
  ?." Host is gforth
[THEN]
IFORTH-HOST [IF]
  ?." Host is iForth
[THEN]


\ *** POLYNOMIAL UTILITIES

\ factor from fsl-util.fs
: <%>  ( "dig.dig...dig" -- f: r )
  bl word count >float 0= ABORT" NaN" ;

: z%   ( "xdig.xdig...xdig ydig.ydig...ydig" -- f: x y )
  <%> <%> STATE @ IF POSTPONE zliteral THEN ; IMMEDIATE

(
Note that FSL arrays are represented on the stack in natural
readable order as [b0 ... bn], and in memory as b0 ... bn, with
b0 at lower memory.

For polynomials, p = bn z^n + ... + b0 z^0.  Julian represents
them in memory also as b0 ... bn, with b0 at lower memory, and
we put them on the stack in natural readable order as [bn ... b0].

In the following, n+1 is the dimension of either kind of array,
sometimes call dim.  Arrays work off of the dim, while
polynomials work of the degree, one less than their dim when
stored as arrays.
)

: }za!  ( f: z0 ... zn s: 'a n+1 -- )
  dup 0 ?DO 1- 2dup 2>r } z! 2r> LOOP 2drop ;

: }za@ ( 'a n+1 -- r: z0 ... zn )
  0 ?DO ( 'a) dup i } z@ LOOP drop ; 

: }zp!  ( r: zn ... z0 s: 'a n -- )
  1+ 0 DO dup i } z! LOOP drop ;

: }zp@ ( 'a n -- r: zn ... z0 )
  dup 1+ 0 ?DO ( 'a n) 2dup i - } z@ LOOP 2drop ;

2 FLOATS DARRAY zha{

40 2 FLOATS ARRAY zpac{  \ polynomial accumulator

: }0za  ( 'a dim -- )  1+ 0 DO 0e 0e dup i } z! LOOP drop ;

: 0zpac  ( dim -- )  zpac{ swap }0za ;

: }zpac!  ( 'a deg -- )
(
Copy the polynomial into the zpac.
)
  1+ 0 DO ( a{) dup i } z@ zpac{ i } z! LOOP drop ;

: zp+  ( 'a 'b n -- )
(
Add the two polynomials of degree n into the polynomial
accumulator.  Either or both of 'a and 'b may be the
accumulator.
)
  1+ 0 DO 2dup i } z@ i } z@ z+ zpac{ i } z! LOOP 2drop ;

: zp-  ( 'a 'b n -- )
(
Subtract the two polynomials of degree n into the polynomial
accumulator.  Either or both of 'a and 'b may be the
accumulator.
)
  >r swap r> 1+ 0
  DO 2dup i } z@ i } z@ z- zpac{ i } z! LOOP 2drop ;

: zp*  ( 'a m 'b n -- )
(
Compute the product of the two polynomials of degree m and n
into the polynomial accumulator.  Neither 'a nor 'b is allowed
to be the accumulator.

If m >= n:

  P = jSum_0^[m+n] z^j iSum_[max 0,j-m]^[min n,j] a_[j-i] b_i
)
  dup 3 pick ( n m) > IF 2swap THEN
  LOCALS| n b{ m a{ |  \ m >= n
  m n + dup 0zpac
  ( m+n) 1+ 0
  DO
    i n min 1+  0 i m - max
    DO    \ coeff of z^j
      a{ j i - } z@ b{ i } z@ z*
      zpac{ j } dup z@ z+ z!
    LOOP
  LOOP ;

\ display complex polynomial
: }pz.   ( 'a deg -- )  1+ }z. ;
: }pzs.  ( 'a deg -- )  1+ }zs. ;


\ *** TESTS

TESTING  z% }za! }za@ }zp! }zp@ 0zpac }zpac! zp+ zp- zp*

\ scratch
20 2 floats ARRAY z1{
20 2 floats ARRAY z2{

z{ z% 1.0 0.0 -> 1e 0e }z

z{ z% 1.0 0.0 z% 2.0 0.5 z% 3.0 0.25 z1{ 3 }za! -> }z
z{ z1{ 0 } z@ -> 1e 0e }z
z{ z1{ 1 } z@ -> 2e 0.5e }z
z{ z1{ 2 } z@ -> 3e 0.25e }z
z{ z1{ 3 }za@ -> 1e 0e 2e 0.5e 3e 0.25e }z

z{ 4e 40e 5e 50e 6e 60e z2{ 2 }zp! -> }z
z{ z2{ 0 } z@ -> 6e 60e }z
z{ z2{ 1 } z@ -> 5e 50e }z
z{ z2{ 2 } z@ -> 4e 40e }z
z{ z2{ 2 }zp@ -> 4e 40e 5e 50e 6e 60e }z

z{ 1e 2e 1e 2e 1e 2e zpac{ 3 }za! -> }z
z{ 1 0zpac -> }z
z{ zpac{ 2 }za@ -> 0e 0e 0e 0e }z
z{ zpac{ 2 } z@ -> 1e 2e }z

z{ 2e 0e 3e 0e z1{ 1 }zp!
   4e 0e 5e 0e z2{ 1 }zp! -> }z
z{ z1{ 1 z2{ 1 zp* -> }z
z{ zpac{ 2 }zp@ -> 8e 0e 22e 0e 15e 0e }z

\  p1 = 0+i1 z^2 + 0+i2 z + 0+i3
\  p2 =            0+i4 z + 0+i5
z{ 0e 1e 0e 2e 0e 3e z1{ 2 }zp!
   0e 4e 0e 5e z2{ 1 }zp! -> }z

z{ z1{ 0 }zpac! -> }z
z{ zpac{ 0 }zp@ -> z1{ 0 }zp@ }z

z{ z1{ 2 }zpac! -> }z
z{ zpac{ 2 }zp@ -> z1{ 2 }zp@ }z

z{ z1{ z2{ 0 zp+ -> }z
z{ zpac{ 0 }zp@ -> 0e 8e }z
z{ z1{ z2{ 1 zp+ -> }z
z{ zpac{ 1 }zp@ -> 0e 6e 0e 8e }z

z{ z1{ z2{ 0 zp- -> }z
z{ zpac{ 0 }zp@ -> 0e -2e }z
z{ z1{ z2{ 1 zp- -> }z
z{ zpac{ 1 }zp@ -> 0e -2e 0e -2e }z

z{ z1{ 2  z2{ 1 zp* -> }z
\ pac = -4+i0 z^3 - 13+i0 z^2 - 22+i0 z - 15+i0
z{ zpac{ 3 }zp@ -> -4e 0e -13e 0e -22e 0e -15e 0e }z
z{ 4 0zpac -> }z
z{ z2{ 1 z1{ 2 zp* -> }z
z{ zpac{ 3 }zp@ -> -4e 0e -13e 0e -22e 0e -15e 0e }z

\ DON'T CHANGE THE ARRAY Z1{, IT'S USED JUST BELOW!

TESTING  }zpeval }zmov }>zp }0zp
(
These tests use the polynomial in Z1{  just above.
)

0e 0e zconstant 0+i0

z{ 0+i0  z1{ 0 }zpeval -> 0e 3e }z
z{ 0e 1e z1{ 0 }zpeval -> 0e 3e }z

z{ 0+i0  z1{ 1 }zpeval ->  0e  3e }z
z{ 0e 1e z1{ 1 }zpeval -> -2e  3e }z
z{ 2e 0e z1{ 1 }zpeval ->  0e  7e }z

z{ 0+i0  z1{ 2 }zpeval ->  0e  3e }z
z{ 0e 1e z1{ 2 }zpeval -> -2e  2e }z
z{ 2e 0e z1{ 2 }zpeval ->  0e 11e }z

z{ 1e 0e 4e 0e -6e 0e -4e 0e -7e 0e -48e 0e 60e 0e
   z1{ 6 }zp! -> }z

z{ z2{ 7 }0za z1{ z2{ 0 }zmov -> }z
z{ z1{ 0 }zp@ -> z2{ 0 }zp@ }z
z{ z2{ 1 } z@ -> 0e 0e }z

z{ z1{ z2{ 1 }zmov -> }z
z{ z1{ 1 }zp@ -> z2{ 1 }zp@ }z
z{ z2{ 2 } z@ -> 0e 0e }z

z{ z1{ z2{ 5 }zmov -> }z
z{ z1{ 5 }zp@ -> z2{ 5 }zp@ }z
z{ z2{ 6 } z@ -> 0e 0e }z

2 2 floats ARRAY [z-a]{  \ dedicated to }>zp
: }>zp  ( 'p n r: z1 ... zn -- )
(
Form the polynomial [z-z1]...[z-zn] and store in the polynomial
p of degree n.
)
  1e 0e zswap znegate over ( p{ ) 1 }zp!
  1 ?DO   \ i is the degree so far
    1e 0e zswap znegate [z-a]{ 1 }zp!
    ( p{ ) dup i [z-a]{ 1 zp* 
    zpac{ over ( p{ ) i 1+ }zmov
  LOOP drop ;

z{ 1e 0e z1{ 1 }>zp -> }z
z{ z1{ 1 }zp@ -> 1e 0e -1e -0e }z

z{ 1e 0e 0e 2e z1{ 2 }>zp -> }z
z{ z1{ 2 }zp@ -> 1e 0e -1e -2e 0e 2e }z

z{ 1e 0e 1e 1e 2e 3e z1{ 3 }>zp -> }z
z{ 1e 0e z1{ 3 }zpeval -> 0e 0e }z
z{ 1e 1e z1{ 3 }zpeval -> 0e 0e }z
z{ 2e 3e z1{ 3 }zpeval -> 0e 0e }z

\ clear coeffs
: }0zp  ( 'p deg -- )
  1+ 0 DO  0e 0e ( p{) dup i } z! LOOP drop ;

z{ 1e 1e z1{ 0 } z! -> }z
z{ 1e 1e z1{ 1 } z! -> }z
z{ 1e 1e z1{ 2 } z! -> }z
z{ 1e 1e z1{ 3 } z! -> }z

z{ z1{ 2 }0zp -> }z
z{ z1{ 3 } z@ -> 1e 1e }z
z{ z1{ 2 } z@ -> 0e 0e }z
z{ z1{ 1 } z@ -> 0e 0e }z
z{ z1{ 0 } z@ -> 0e 0e }z

TESTING  }zsynth

0 [IF]  \ from lagroots.fs
Test case for synthetic division:

        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 zs.

        CR b{ 5 }zs.

Results should be

 7.63200E3  + i0.00000E-1  ok

0  1.52500E3 + i0.00000E0
1  3.06000E2 + i0.00000E0
2  6.10000E1 + i0.00000E0
3  1.50000E1 + i0.00000E0
4  3.00000E0 + i0.00000E0
[THEN]

\ Julian's synthetic division polynomial
  6 2 floats ARRAY sdp{
  z%   7.0 0.0 sdp{ 0 } z!
  z%  -5.0 0.0 sdp{ 1 } z!
  z%   1.0 0.0 sdp{ 2 } z!
  z% -14.0 0.0 sdp{ 3 } z!
  z%   0.0 0.0 sdp{ 4 } z!
  z%   3.0 0.0 sdp{ 5 } z!

\ set up his synthetic division test
z{ sdp{ a{ 5 }zmov -> }z

SET-FT-MODE-EXACT
t{ a{ b{ 5 5e0 0e }zsynth -> 7.63200E3  0.00000E-1 }t
t{    5e0 0e a{ 5 }zpeval -> 7.63200E3  0.00000E-1 }t

VERBOSE @ [IF]
?." Expected coffficients::
.( 0  1.52500E3 + i 0.00000E-1) cr
.( 1  3.06000E2 + i 0.00000E-1) cr
.( 2  6.10000E1 + i 0.00000E-1) cr
.( 3  1.50000E1 + i 0.00000E-1) cr
.( 4  3.00000E0 + i 0.00000E-1) cr
?." Test results:
b{ 4 }pzs.
[THEN]

\ expected results
\ no change in dividend
z{ a{ 5 }zp@ -> 3e 0e 0e 0e -14e 0e 1e 0e -5e 0e 7e 0e }z
\ quotient
z{ b{ 4 }zp@ -> 3e 0e 15e 0e 61e 0e 306e 0e 1525e 0e }z
\ divisor
z{ 1e 0e -5e 0e z1{ 1 }zp! -> }z
\ remainder
z{ z2{ 6 }0za -> }z
z{ 5e 0e a{ 5 }zpeval z2{ 0 } z! -> }z
\ recovery of original polynomial
z{ b{ 4 z1{ 1 zp* zpac{ 5 }zp@ z1{ 5 }zp!
   z1{ z2{ 5 zp+ zpac{ 5 }zp@ -> a{ 5 }zp@ }z

\ a{} at degree 4
t{ a{ b{ 4 5e0 0e }zsynth -> 5e 0e a{ 4 }zpeval }t
\ no change in dividend
z{ a{ 4 }zp@ -> 0e 0e -14e 0e 1e 0e -5e 0e 7e 0e }z
\ divisor
z{ 1e 0e -5e 0e z1{ 1 }zp! -> }z
\ remainder
z{ z2{ 5 }0za -> }z
z{ 5e 0e a{ 4 }zpeval z2{ 0 } z! -> }z
\ recovery of original polynomial
z{ b{ 3 z1{ 1 zp* zpac{ 4 }zp@ z1{ 4 }zp!
   z1{ z2{ 4 zp+ zpac{ 4 }zp@ -> a{ 4 }zp@ }z

\ a{} at degree 1
t{ a{ b{ 1 5e0 0e }zsynth -> 5e 0e a{ 1 }zpeval }t
\ no change in dividend
z{ a{ 1 }zp@ -> -5e 0e 7e 0e }z
\ divisor
z{ 1e 0e -5e 0e z1{ 1 }zp! -> }z
\ remainder
z{ z2{ 1 }0za -> }z
z{ 5e 0e a{ 1 }zpeval z2{ 0 } z! -> }z
\ recovery of original polynomial
z{ b{ 0 z1{ 1 zp* zpac{ 1 }zp@ z1{ 1 }zp!
   z1{ z2{ 1 zp+ zpac{ 1 }zp@ -> a{ 1 }zp@ }z

TESTING  roots quadroots

0 [IF]  \ from lagroots.fs
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!

       d{ 6 10e0 0e0 1e-9 roots
       CR d{ 6 }zs.

Should give

0 -5.00000E0 - i4.55636E-13
1  7.99361E-14 - i1.73205E0
2 -2.00000E0 + i1.89263E-12
3 -1.72719E-16 + i1.73205E0
4  1.00000E0 + i0.00000E0
5  2.00000E0 + i0.00000E0
[THEN]

\ Julian's roots poly
7 2 floats ARRAY rp{
z{ 1e 0e 4e 0e -6e 0e -4e 0e -7e 0e -48e 0e 60e 0e
   rp{ 6 }zp! -> }z

\ set up test
z{ rp{ a{ 6 }zmov -> }z
\ find roots, epsilon = 1e-9
t{ d{ 6 10e 0e 1e-9 roots -> }t
VERBOSE @ [IF]
?." Expected roots:
.( 0 -5.00000E0 - i4.55636E-13) cr
.( 1  7.99361E-14 - i1.73205E0) cr
.( 2 -2.00000E0 + i1.89263E-12) cr
.( 3 -1.72719E-16 + i1.73205E0) cr
.( 4  1.00000E0 + i0.00000E0) cr
.( 5  2.00000E0 + i0.00000E0) cr
?." Test results:
z{ d{ 6 }zs. -> }z
[THEN]

\ set FT-ABS-ERROR to 0e to see errors
\ max error 1.1834E-08
1.185e-8 FT-ABS-ERROR f!
SET-FT-MODE-ABS
z{ d{ 0 } z@ rp{ 6 }zpeval -> 0e 0e }z
z{ d{ 1 } z@ rp{ 6 }zpeval -> 0e 0e }z
z{ d{ 2 } z@ rp{ 6 }zpeval -> 0e 0e }z
z{ d{ 3 } z@ rp{ 6 }zpeval -> 0e 0e }z
z{ d{ 4 } z@ rp{ 6 }zpeval -> 0e 0e }z
z{ d{ 5 } z@ rp{ 6 }zpeval -> 0e 0e }z

z{ 2e 0e 0e 3e z1{ 2 }>zp -> }z
z{ z1{ quadroots -> 2 0e 3e 2e 0e }z \ order happens to reverse
z{ z1{ a{ 2 }zmov -> }z
z{ d{ 2 10e 0e 1e-9 roots -> }z
z{ d{ 0 } z@ z1{ 2 }zpeval -> 0e 0e }z
z{ d{ 1 } z@ z1{ 2 }zpeval -> 0e 0e }z
z{ d{ 1 } z@ -> 2e 0e }z
z{ d{ 0 } z@ -> 0e 3e }z

z{ 0e 0e 1e 0e 2e 3e z1{ 2 }zp! -> }z
z{ z1{ quadroots -> 1 -2e -3e }z

z{ z1{ 3 }0za -> }z
z{ z1{ quadroots -> 0 }z

z{ 1e 2e  z1{ 0 } z! -> }z
z{ z1{ quadroots -> 0 }z

8 2 floats ARRAY 8thr{  \ 8th roots of unity
9 2 floats ARRAY cp{    \ circle polynomial
1e 2e fsqrt f/  FCONSTANT  1/sqrt2
1/sqrt2 fnegate FCONSTANT -1/sqrt2

\ load 8thr{ and cp{
z{ 1e 0e  1e       1/sqrt2 0e  1e -1/sqrt2  1/sqrt2
  -1e 0e -1/sqrt2 -1/sqrt2 0e -1e  1/sqrt2 -1/sqrt2
      8thr{ 8 }za! -> }z
z{ 8thr{ 8 }za@ cp{ 8 }>zp -> }z

\ find roots, epsilon = 1e-16
z{ cp{ a{ 8 }zmov -> }z
t{ d{ 8 2e 0e 1e-16 roots -> }t

\ set FT-ABS-ERROR to 0e to see the errors
\ max error 1.1546E-14
1.16e-14 FT-ABS-ERROR f!
SET-FT-MODE-ABS
z{ d{ 0 } z@ cp{ 8 }zpeval -> 0e 0e }z
z{ d{ 1 } z@ cp{ 8 }zpeval -> 0e 0e }z
z{ d{ 2 } z@ cp{ 8 }zpeval -> 0e 0e }z
z{ d{ 3 } z@ cp{ 8 }zpeval -> 0e 0e }z
z{ d{ 4 } z@ cp{ 8 }zpeval -> 0e 0e }z
z{ d{ 5 } z@ cp{ 8 }zpeval -> 0e 0e }z
z{ d{ 7 } z@ cp{ 8 }zpeval -> 0e 0e }z

\ set to true to see the order
0 [IF]
.( 8th roots supplied:) cr
8thr{ 8 }zs.
.( 8th roots computed:) cr
d{ 8 }zs.
[THEN]

0 [IF] \ here's the printout of the above
8th roots supplied:
0  1.00000E+00 + i0.00000E+00 
1  1.00000E+00 + i7.07107E-01 
2  0.00000E+00 + i1.00000E+00 
3 -7.07107E-01 + i7.07107E-01 
4 -1.00000E+00 + i0.00000E+00 
5 -7.07107E-01 - i7.07107E-01 
6  0.00000E+00 - i1.00000E+00 
7  7.07107E-01 - i7.07107E-01 
8th roots computed:
0 -5.55112E-16 + i1.00000E+00 
1  1.00000E+00 + i7.07107E-01 
2 -7.07107E-01 + i7.07107E-01 
3 -1.00000E+00 - i9.50277E-16 
4 -7.07107E-01 - i7.07107E-01 
5 -2.45088E-16 - i1.00000E+00 
6  7.07107E-01 - i7.07107E-01 
7  1.00000E+00 - i1.02782E-16 
[THEN]

\ set FT-ABS-ERROR to 0e to see results and goals
\ 14 SET-PRECISION
1.34e-15 FT-ABS-ERROR f!
z{ d{ 0 } z@ -> 8thr{ 2 } z@ }z
z{ d{ 1 } z@ -> 8thr{ 1 } z@ }z
z{ d{ 2 } z@ -> 8thr{ 3 } z@ }z
z{ d{ 3 } z@ -> 8thr{ 4 } z@ }z
z{ d{ 4 } z@ -> 8thr{ 5 } z@ }z
z{ d{ 5 } z@ -> 8thr{ 6 } z@ }z
z{ d{ 6 } z@ -> 8thr{ 7 } z@ }z
z{ d{ 7 } z@ -> 8thr{ 0 } z@ }z

TESTING  roots exception

[UNDEFINED] F2DROP [IF] : F2DROP  FDROP FDROP ; [THEN]
\ shouldn't matter what the polynomial is, but nevertheless...
z{ rp{ a{ 6 }zmov -> }z
\ t{ d{ 1 10e 0e 1e-9 roots -> }t
t{ d{ 1 10e 0e 1e-9 ' roots CATCH fdrop f2drop nip nip -> -2 }t
\ t{ d{ 0 10e 0e 1e-9 roots -> }t
t{ d{ 0 10e 0e 1e-9 ' roots CATCH fdrop f2drop nip nip -> -2 }t
IFORTH-HOST [IF] ?.cr [THEN]  \ for return stack dump

?." Done testing lagroots
?.errors GFORTH-HOST [IF] ?.cr [THEN]