(
This file contains the subset of the FSL complex.seq needed to
define ZABS for testing.  It is renamed here as |Z|.  There are
also cosmetic modifications to run under ANS Forth without some
of the special FSL utilities.  dnw 12Jan03
)
\ COMPLEX   A definition of a complex number type and associated operations on it.
\           These are the words that might care about the implementation details
\           of the Complex type.  This file specifies the SYNTAX of operations
\           on the complex data type.  The internals of the implementation here
\           are typical but should not be assumed by other programs as they
\           very likely could be rewritten for efficiency on specfic platforms
\           or compilers.

\ For the sake of possible optimizations, the following words are also defined here
\ in spite of the fact that they can be defined without knowledge of the
\ details of implementation:
\ Z+      sum of two complex numbers
\ Z-      difference of two complex numbers
\ Z*      product of two complex numbers
\ ZABS    complex magnitude
\ Z/      complex division


\ The Z/ routine does the complex division (x1,y1)/(x2,y2)
\ The ZABS routine finds Magnitude of the pair(x,y)

\ If the CONSTANT PLAUGER-CODE? is defined to be TRUE then,
\     the routines ZABS and Z/ are designed to preserve as much precision
\     as possible. They are based upon the algorithms described in:
\        Plauger, P.J., 1994; Complex Made Simple, Embedded Systems Programming,
\        July, pp. 85-88
\ otherwise the obvoius straightfoward algorithm is used.


\ This is an ANS Forth program requiring:
\      1. The Floating-Point word set
\      2. The immediate word '%' which takes the next token
\         and converts it to a floating-point literal
\         : %   BL WORD  COUNT  >FLOAT  0= ABORT" NAN"
\               STATE @  IF POSTPONE FLITERAL  THEN  ; IMMEDIATE
\
\      3. The word 'F,' to store an floating point number at (FALIGNED) HERE
\         : F,  ( f: x -- )        HERE 1 FLOATS ALLOT F! ;

\      4. Uses words 'Private:', 'Public:' and 'Reset_Search_Order'
\         to control the visibility of internal code.

\  (c) Copyright 1994 Everett F. Carter.  Permission is granted by the
\  author to use this software for any application provided this
\  copyright notice is preserved.

CR .( COMPLEX           V1.4           21 December 1994   EFC )

TRUE CONSTANT PLAUGER-CODE?             \ set to TRUE or FALSE as desired
                                        \ before compiling

FVARIABLE t-real
FVARIABLE t-imag
FVARIABLE scale-factor

PLAUGER-CODE? [IF]


2e0 FSQRT         FCONSTANT sqrt2.0
sqrt2.0 1.4142e0  F- FCONSTANT little-bits

: zabs-xy ( -- , f: x y -- x y )

        FABS FSWAP FABS        \ first set absolute values ( x y -- y x )

        \ now put smallest on top of the stack
        FOVER FOVER F< IF FSWAP THEN
        
        
;

: set-scale-factor ( -- , f: x y -- x y )

           FOVER 1e0 F> IF
                            4e0 F* FSWAP
                            4e0 F* FSWAP
                            0.25e0 scale-factor F!
                          ELSE
                            0.25e0 F* FSWAP
                            0.25e0 F* FSWAP
                            4e0 scale-factor F!
                          THEN
;

: zabs-for-small-y ( -- , f: x y -- z )
                 FOVER FOVER F/
                 FDUP FDUP F*
                 1e0 F+ FSQRT
                 F+ F/ F+
;

: zabs-for-mid-xy ( -- , f: x y -- z )
                FOVER FOVER F- FOVER F/
                FDUP FDUP 2e0 F+ F*
                FDUP 2e0 F+ FSQRT

                sqrt2.0 F+
                F/
                little-bits F+
                F+
                2.4142e0 F+
                F/ F+
;

: normalize ( -- , f: z -- z )
       scale-factor F@ F*        
;

: zero-divide-error   FDROP FDROP FDROP FDROP ." zero divide error" ABORT ;

: smaller-real-div ( -- , f: x1 y1 x2 y2 -- x3 y3 )
       FOVER FOVER F/
       FROT FOVER F* FROT F+

       FDUP F0= IF zero-divide-error
                ELSE
                   scale-factor F!
                   FROT FOVER FOVER F* t-real F!
                   FROT FDUP t-real F@ F+ scale-factor F@ F/ t-real F!
                   FROT F* FSWAP F- scale-factor F@ F/
                   t-real F@ FSWAP
                THEN
;

: smaller-imag-div ( -- , f: x1 y1 x2 y2 -- x3 y3 )
       FSWAP FOVER FOVER F/
       FROT FOVER F* FROT F+

       FDUP F0= IF zero-divide-error
                ELSE
                   scale-factor F!
                   FROT FOVER FOVER F* t-imag F!
                   FROT FDUP t-imag F@ F- scale-factor F@ F/ t-imag F!
                   FROT F* F+ scale-factor F@ F/
                   t-imag F@
                THEN

;

: ZABS-SCALED ( -- , f: x y -- z )

         0e0 scale-factor F!
         zabs-xy

         FOVER F0= IF FDROP
                   ELSE
                        set-scale-factor
                        FOVER FOVER FOVER FSWAP F-
                        F= IF FDROP
                           ELSE
                               FOVER FOVER F-
                               FOVER
                               F< IF
                                      zabs-for-small-y
                                   ELSE
                                      zabs-for-mid-xy
                                   THEN
                           THEN
                   THEN    

;

[THEN]

: ZDUP ( -- ) ( f: re im -- re im re im )     FOVER FOVER ;

: Z->R,I ( -- ) ( F: x y -- re im )     ; IMMEDIATE
: Z->I,R ( -- ) ( F: x y -- im re )   FSWAP ;
: R,I->Z ( -- ) ( F: x y -- re im )     ; IMMEDIATE

PLAUGER-CODE? [IF]

: |Z| ( -- , f: z1 -- r )
	    Z->R,I

            ZABS-SCALED
            normalize
;

: Z/ ( -- , f: z1 z2 -- z3 )

         FOVER FABS FOVER FABS
         F< IF   smaller-real-div
            ELSE smaller-imag-div
            THEN

         R,I->Z
;

[ELSE]

: |Z| ( -- , f: z1 -- r )
	Z->R,I
	FDUP F* FSWAP FDUP F* F+  FSQRT
;

: Z/ ( -- , f: z1 z2 -- z3 )
        
	ZDUP   

	Z->R,I
        FDUP F* FSWAP FDUP F* F+ scale-factor F!

	Z->R,I
	t-imag F!   t-real F!

        Z->R,I
	FOVER t-real F@ F*
        FOVER t-imag F@ F* F+   scale-factor F@ F/

	FROT t-imag F@ F* FNEGATE
        FROT t-real F@ F* F+    scale-factor F@ F/

	R,I->Z

;

[THEN]