\ quadratic.fs
\
\ Solve for the two complex solutions of the quadratic equation:
\
\   a*x^2 + b*x + c = 0
\
\ Assume a, b, and c are real. The two complex roots are 
\ z1 = x1r + i*x1i and z2 = x2r + i*x2i
\
\  2002, K. Myneni, [email protected]
\  Permission is granted to modify and use this code
\  for any application, provided this notice is preserved.
\
\ Notes:
\
\   1. Ensure that argument "a" is not zero, or an infinity will result;
\      the correct solution of the simple linear equation will not be
\      given.
\  
\   2. The returned roots are unordered. 
\   
\ Revisions:
\
\   2002-11-02  km; first version.
\   2003-10-25  Christopher Brannon; fixed problem with calculation
\                 of complex roots.
\   2007-11-04  km; revised comments; added test code; save and restore base. 
\   2009-08-05  km; revised to preserve accuracy when the product a*c is
\                   much less than b^2; see [1]. Added new test case.
\
\ References:
\
\ 1. W.H. Press, et. al., Numerical Recipes in C, 2nd ed., pp. 183--184,
\    eqns. 5.6.4 and 5.6.5.
include lib/ansfloat.4th
include lib/ansfpio.4th

:MACRO -> => ;
:MACRO TEST-CODE? true ;

CR .( QUADRATIC         V1.1          05 August    2009   KM )
BASE @ DECIMAL

fvariable qa
fvariable 2qa
fvariable qb
fvariable qc

: solve_quadratic ( F: a b c -- x1r x1i x2r x2i )
                  ( F: a b c -- z1  z2 )
    qc f! qb f! fdup qa f! F% 2e f* 2qa f!
    qb f@ fdup f* F% 4e qa f@ f* qc f@ f* f-   \ square root term
    fdup f0<
    IF
	\ complex conjugate roots
	fabs  fsqrt   2qa f@  f/           \ imaginary part 
        qb f@ fnegate 2qa f@  f/           \ real part
	fswap 
        fover fover fnegate                \ complex conjugate
    ELSE
	\ two real roots
	fsqrt qb f@ fdup f0< IF f- ELSE f+ fnegate THEN F% 2e f/ 
	fdup  qc f@ fswap f/ F% 0e   frot qa f@ f/ F% 0e
    THEN
;

BASE !

TEST-CODE? [IF]   \ test code ==============================================
[undefined] T{      [IF]  include lib/ttester.4th  [THEN]    
BASE @ DECIMAL

F% 1e-15  rel-near F!
F% 1e-256 abs-near F!
set-near

\ Examples from:
\
\    http://www.purplemath.com/modules/quadform.htm

F% -2e F% 3e F/       fconstant -2/3
F% -3e F% 2e F/       fconstant -3/2
F%  5e fsqrt          fconstant sqrt{5}
F%  2e fsqrt F% 3e F/ fconstant sqrt{2}/3
F%  3e fsqrt F% 2e F/ fconstant sqrt{3}/2
F% 10e fsqrt F% 2e F/ fconstant sqrt{10}/2

CR
\ TESTING SOLVE_QUADRATIC
t{ F% 1e F%  3e F% -4e solve_quadratic ->  F% 1e               F% 0e   F% -4e               F% 0e   rrrr}t
t{ F% 2e F% -4e F% -3e solve_quadratic ->  F% 1e sqrt{10}/2 F- F% 0e   F%  1e sqrt{10}/2 F+ F% 0e   rrrr}t
t{ F% 1e F% -2e F% -4e solve_quadratic ->  F% 1e sqrt{5} F-    F% 0e   f%  1e sqrt{5} F+    F% 0e   rrrr}t
t{ F% 9e F% 12e F%  4e solve_quadratic ->  -2/3                F% 0e   -2/3                 F% 0e   rrrr}t
t{ F% 3e F%  4e F%  2e solve_quadratic ->  -2/3      sqrt{2}/3   -2/3 sqrt{2}/3 fnegate rrrr}t
t{ F% 1e F%  3e F%  3e solve_quadratic ->  -3/2      sqrt{3}/2   -3/2 sqrt{3}/2 fnegate rrrr}t

\ Test case which loses accuracy with ordinary quadratic formula:
\ 
\        x^2 + x + c = 0
\
\ when c << 1, the approximate solution is  x1 = -c,  x2 = -1 + c
\
t{ F% 1e F% 1e F% 1e-17 solve_quadratic -> F% -1e-17           F% 0e    F% -1e F% 1e-17 f+     F% 0e   rrrr}t


BASE !
[THEN]