\ COMPLEX2 An alternative (demonstration) definition of a complex number type \ and associated operations on it. This version uses a seperate \ complex number stack. \ 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 .( COMPLEX2 V1.2 21 December 1994 EFC ) Private: FALSE CONSTANT PLAUGER-CODE? \ set to TRUE or FALSE as desired \ before compiling \ a COMPLEX STACK where the numbers are stored 16 CONSTANT ZSTACK-SIZE VARIABLE Z-PTR -1 Z-PTR ! ZSTACK-SIZE 1+ FLOAT ARRAY zre{ ZSTACK-SIZE 1+ FLOAT ARRAY zim{ : Z-INC ( -- ptr ) Z-PTR @ 1+ DUP ZSTACK-SIZE > ABORT" Complex stack error, too deep " DUP Z-PTR ! ; : Z-DEC ( -- ptr ) Z-PTR @ DUP 0< ABORT" Complex stack error, empty " DUP 1- Z-PTR ! ; : ZPUSH ( -- ) ( f: re im -- ) Z-INC zim{ OVER } F! zre{ SWAP } F! ; : (zre-stack) Z-INC zre{ SWAP } F! ; : (zim-stack) zim{ Z-PTR @ } F! ; : ZPOP ( -- ) ( f: -- re im ) Z-DEC zre{ OVER } F@ zim{ SWAP } F@ ; : ZRPUSH ( -- ) ( f: im re -- ) Z-INC zre{ OVER } F! zim{ SWAP } F! ; : ZRPOP ( -- ) ( f: -- im re ) Z-DEC zim{ OVER } F@ zre{ SWAP } F@ ; \ pushes following value to the complex stack as the real part : (zre%) BL WORD COUNT >FLOAT 0= ABORT" NAN" STATE @ IF POSTPONE FLITERAL POSTPONE (zre-stack) ELSE (zre-stack) THEN ; IMMEDIATE \ pushes following value to the complex stack as the imaginary part : (zim%) BL WORD COUNT >FLOAT 0= ABORT" NAN" STATE @ IF POSTPONE FLITERAL POSTPONE (zim-stack) ELSE (zim-stack) THEN ; IMMEDIATE FVARIABLE scale-factor FVARIABLE t-real FVARIABLE t-imag PLAUGER-CODE? [IF] % 2.0 FSQRT FCONSTANT sqrt2.0 sqrt2.0 % 1.4142 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 F2DUP F< IF FSWAP THEN ; : set-scale-factor ( -- , f: x y -- x y ) FOVER % 1.0 F> IF % 4.0 F* FSWAP % 4.0 F* FSWAP % 0.25 scale-factor F! ELSE % 0.25 F* FSWAP % 0.25 F* FSWAP % 4.0 scale-factor F! THEN ; : zabs-for-small-y ( -- , f: x y -- z ) F2DUP F/ FDUP FDUP F* % 1.0 F+ FSQRT F+ F/ F+ ; : zabs-for-mid-xy ( -- , f: x y -- z ) F2DUP F- FOVER F/ FDUP FDUP % 2.0 F+ F* FDUP % 2.0 F+ FSQRT sqrt2.0 F+ F/ little-bits F+ F+ % 2.4142 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 ) F2DUP F/ FROT FOVER F* FROT F+ FDUP F0= IF zero-divide-error ELSE scale-factor F! FROT F2DUP 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 F2DUP F/ FROT FOVER F* FROT F+ FDUP F0= IF zero-divide-error ELSE scale-factor F! FROT F2DUP 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 ) % 0.0 scale-factor F! zabs-xy FOVER F0= IF FDROP ELSE set-scale-factor F2DUP FOVER FSWAP F- F= IF FDROP ELSE F2DUP F- FOVER F< IF zabs-for-small-y ELSE zabs-for-mid-xy THEN THEN THEN ; [THEN] Public: : ZDEPTH Z-PTR @ 1+ ; 2 FLOATS CONSTANT COMPLEX \ the size of the type : Z@ ( addr -- ) ( f: -- z ) DUP F@ FLOAT+ F@ ZPUSH ; : Z! ( addr -- ) ( f: z -- ) DUP FLOAT+ ZPOP F! F! ; : ZDUP ( -- ) ( f: z -- z1 z2 ) zre{ Z-PTR @ } F@ zim{ Z-PTR @ } F@ ZPUSH ; : ZSWAP ( -- ) ( f: z1 z2 -- z2 z1 ) Z-PTR @ DUP 1 < ABORT" Complex stack error, need at least two numbers " zre{ OVER } F@ zre{ OVER 1- } F@ zre{ OVER } F! zre{ OVER 1- } F! zim{ OVER } F@ zim{ OVER 1- } F@ zim{ OVER } F! zim{ SWAP 1- } F! ; : ZDROP ( -- ) ( f: z -- ) Z-DEC DROP ; : REAL ( -- ) ( f: re im -- re ) ZPOP FDROP ; : IMAG ( -- ) ( f: re im -- im ) ZRPOP FDROP ; : ZCONJUGATE ( -- ) ( f: z -- z* ) Z-PTR @ zim{ OVER } F@ FNEGATE zim{ SWAP } F! ; : REAL>Z ( -- ) ( f: x -- z ) % 0.0 ZPUSH ; \ convert float to complex : Z->R,I ( -- ) ( F: z -- re im ) ZPOP ; : Z->I,R ( -- ) ( F: z -- im re ) ZRPOP ; : R,I->Z ( -- ) ( F: re im -- z ) ZPUSH ; : I,R->Z ( -- ) ( F: im re -- z ) ZRPUSH ; \ convert float to pure imaginary complex value : IMAG>Z ( -- ) ( F: r -- z ) % 0.0 ZRPUSH ; : ZVARIABLE ( "name" -- ) CREATE FALIGN COMPLEX ALLOT DOES> FALIGNED ; : ZCONSTANT ( "name" -- ) CREATE Z->I,R F, F, \ Note: assumes that F, does an FALIGN ! DOES> FALIGNED Z@ ; : Z. ( F: z -- ) ZRPOP ." ( " F. F. ." ) " ; : }ZPRINT ( n addr -- ) \ for printing complex arrays SWAP 0 DO I print-width @ MOD 0= I AND IF CR THEN DUP I } Z@ Z. LOOP DROP ; \ get the next two real numbers as a complex number : Z% POSTPONE (zre%) POSTPONE (zim%) ; IMMEDIATE \ useful complex constants Z% 1.0 0.0 ZCONSTANT 1+0i Z% 0.0 0.0 ZCONSTANT 0+0i Z% 0.0 1.0 ZCONSTANT 0+1i : Z+ ( -- , f: z1 z2 -- z3 ) Z->R,I Z-PTR @ zim{ OVER } DUP F@ F+ F! zre{ SWAP } DUP F@ F+ F! ; : Z- ( -- , f: z1 z2 -- z3 ) Z->R,I Z-PTR @ zim{ OVER } DUP F@ FSWAP F- F! zre{ SWAP } DUP F@ FSWAP F- F! ; : Z* ( --, f: z1 z2 -- z3 ) Z-PTR @ DUP 1 < ABORT" Complex stack error, need at least two numbers " zre{ OVER } F@ zre{ OVER 1- } F@ F* zim{ OVER } F@ zim{ OVER 1- } F@ F* F- zim{ OVER } F@ zre{ OVER 1- } F@ F* zre{ OVER } F@ zim{ SWAP 1- } F@ F* F+ Z-DEC Z-DEC 2DROP R,I->Z ; PLAUGER-CODE? [IF] : ZABS ( -- , f: z -- r ) Z->R,I ZABS-SCALED normalize ; : Z/ ( -- , f: z1 z2 -- z3 ) ZSWAP Z->R,I Z->R,I FOVER FABS FOVER FABS F< IF smaller-real-div ELSE smaller-imag-div THEN R,I->Z ; [ELSE] : ZABS ( -- , 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] Reset_Search_Order