complex.f
\ ANS Forth Complex Arithmetic Lexicon
\ ---------------------------------------------------
\ (c) Copyright 1998 Julian V. Noble. \
\ Permission is granted by the author to \
\ use this software for any application pro- \
\ vided this copyright notice is preserved. \
\ ---------------------------------------------------
\ Environmental dependences:
\ 1. requires FLOAT and FLOAT EXT wordsets
\ 2. assumes separate floating point stack
\ 3. does not construct a separate complex number stack
\ Complex numbers x+iy are stored on the fp stack as ( f: -- x y).
\ Angles are in radians.
\ Polar representation measures angle from the positive x-axis.
\ All Standard words are in uppercase, non-Standard words in lowercase,
\ as a convenience to the user.
MARKER -complex
\ ---------------------------------------- LOAD, STORE
: z@ DUP F@ FLOAT+ F@ ; ( adr -- f: -- z)
: z! DUP FLOAT+ F! F! ; ( adr -- f: z --)
\ ------------------------------------ END LOAD, STORE
\ non-Standard fp words I have found useful
[undefined] s>f [IF] : s>f S>D D>F ; [THEN]
[undefined] f-rot [IF] : f-rot FROT FROT ; [THEN]
[undefined] fnip [IF] : fnip FSWAP FDROP ; [THEN]
[undefined] ftuck [IF] : ftuck FSWAP FOVER ; [THEN]
[undefined] 1/f [IF] : 1/f F1.0 FSWAP F/ ; [THEN]
[undefined] f^2 [IF] : f^2 FDUP F* ; [THEN]
[undefined] fpi [IF] 3.1415926535897932385E0 FCONSTANT fpi [THEN]
[undefined] f0.0 [IF] 0.0E0 FCONSTANT f0.0 [THEN]
[undefined] f1.0 [IF] 1.0E0 FCONSTANT f1.0 [THEN]
\ --------------------------------- MANIPULATE FPSTACK
: z. ( f: x y --) \ emit complex #
FSWAP FS. FDUP F0>
IF ." + i " FS.
ELSE ." - i " FABS FS. THEN
;
: z=0 f0.0 f0.0 ; ( f: -- 0 0)
: z=1 f1.0 f0.0 ; ( f: -- 1 0)
: z=i z=1 FSWAP ; ( f: -- 0 1)
: zdrop FDROP FDROP ; ( f: x y --)
: zdup FOVER FOVER ; ( f: x y -- x y x y)
\ hidden temporary storage for stuff from fpstack
FALIGN HERE VALUE noname 2 FLOATS ALLOT \ ALLOT z variable
: zswap ( f: x y u v -- u v x y)
[ noname ] LITERAL F! f-rot
[ noname ] LITERAL F@ f-rot ;
: zover ( f: x y u v -- x y u v x y )
FROT [ noname FLOAT+ ] LITERAL F! ( f: -- x u v)
FROT FDUP [ noname ] LITERAL F! ( f: -- u v x)
f-rot [ noname FLOAT+ ] LITERAL F@ ( f: -- x u v y)
f-rot [ noname ] LITERAL z@ ( f: -- x y u v x y)
;
: real STATE @ IF POSTPONE FDROP ELSE FDROP THEN ; IMMEDIATE
: imag STATE @ IF POSTPONE fnip ELSE fnip THEN ; IMMEDIATE
: conjg STATE @ IF POSTPONE FNEGATE ELSE FNEGATE THEN ; IMMEDIATE
: freal F@ ; \ used to get the real and imaginary
: fimag FLOAT+ F@ ; \ parts in a real expression
: znip zswap zdrop ;
: ztuck zswap zover ;
: z*f ( f: x y a -- x*a y*a)
FROT FOVER F* f-rot F* ;
: z/f ( f: x y a -- x/a y/a)
1/f z*f ;
: z* ( f: x y u v -- x*u-y*v x*v+y*u)
\ uses the algorithm
\ (x+iy)*(u+iv) = [(x+y)*u - y*(u+v)] + i[(x+y)*u + x*(v-u)]
\ requiring 3 multiplications and 5 additions
zdup F+ ( f: x y u v u+v)
[ noname ] LITERAL F! ( f: x y u v)
FOVER F- ( f: x y u v-u)
[ noname FLOAT+ ] LITERAL F! ( f: x y u)
FROT FDUP ( f: y u x x)
[ noname FLOAT+ ] LITERAL F@ ( f: y u x x v-u)
F*
[ noname FLOAT+ ] LITERAL F! ( f: y u x)
FROT FDUP ( f: u x y y)
[ noname ] LITERAL F@ ( f: u x y y u+v)
F*
[ noname ] LITERAL F! ( f: u x y)
F+ F* FDUP ( f: u*[x+y] u*[x+y])
[ noname ] LITERAL F@ F- ( f: u*[x+y] x*u-y*v)
FSWAP
[ noname FLOAT+ ] LITERAL F@ ( f: x*u-y*v u*[x+y] x*[v-u])
F+ ; ( f: x*u-y*v x*v+y*u)
: z+ FROT F+ f-rot F+ FSWAP ; ( f: a b x y -- a+x b+y)
: znegate FSWAP FNEGATE FSWAP FNEGATE ;
: z- znegate z+ ;
: |z|^2 f^2 FSWAP f^2 F+ ;
\ writing |z| and 1/z as shown reduces overflow probability
: |z| ( f: x y -- |z|)
FABS FSWAP FABS
zdup FMAX f-rot FMIN ( f: max min)
FOVER F/ f^2 1e0 F+ FSQRT F* ;
: 1/z fnegate zdup |z| 1/f FDUP [ noname ] LITERAL F!
z*f [ noname ] LITERAL F@ z*f ;
: z/ 1/z z* ;
: z2/ F2/ FSWAP F2/ FSWAP ;
: z2* F2* FSWAP F2* FSWAP ;
: arg ( f: x y -- arg[x+iy] )
FDUP F0< FSWAP ( : -- y<0 f: -- y x)
FATAN2
IF fpi F2* F+ THEN ;
\ tested September 27th, 1998 - 21:15
: >polar ( f: x+iy -- r phi ) zdup |z| f-rot arg ;
: polar> ( f: r phi -- x+iy ) FSINCOS FROT z*f FSWAP ;
: i* FNEGATE FSWAP ; ( f: x+iy -- -y+ix)
: (-i)* FSWAP FNEGATE ; ( f: x+iy -- y-ix)
: zln >polar FSWAP FDUP F0= ABORT" Can't take ZLN of 0" FLN FSWAP ;
: zexp ( f: z -- exp[z] ) FSINCOS FSWAP FROT FEXP z*f ;
: z^2 zdup z* ;
: z^3 zdup z^2 z* ;
: z^4 z^2 z^2 ;
: z^n ( n -- f: z -- z^n ) \ raise z to integer power
z=1 zswap
DUP 50 < IF
BEGIN DUP 0> WHILE
DUP 1 AND IF ztuck z* zswap THEN z^2
2/
REPEAT zdrop DROP
ELSE FLN S>F z*f zexp THEN ;
: z^ ( f: x y u v -- [x+iy]^[u+iv] ) zswap zln z* zexp ;
: zsqrt ( f: x y -- a b ) \ (a+ib)^2 = x+iy
zdup ( f: -- z z)
|z|^2 ( f: -- z |z|^2 )
FDUP F0= IF FDROP EXIT THEN ( f: -- z=0 )
FSQRT FROT FROT F0< ( f: -- |z| x ) ( -- sgn[y])
ftuck ( f: -- x |z| x )
F- F2/ ( f: -- x [|z|-x]/2 )
ftuck F+ ( f: -- [|z|-x]/2 [|z|+x]/2 )
FSQRT IF FNEGATE THEN ( f: -- [|z|-x]/2 a )
FSWAP FSQRT ; ( f: -- a b)
\ tested September 16th, 1999 - 13:49
\ Complex trigonometric functions
: zcosh ( f: z -- cosh[z] ) zexp zdup 1/z z+ z2/ ;
: zsinh ( f: z -- sinh[z] ) zexp zdup 1/z z- z2/ ;
: ztanh zexp z^2 i* zdup f1.0 F- zswap f1.0 F+ z/ ;
: zcoth ztanh 1/z ;
: zcos ( f: z -- cos[z] ) i* zcosh ;
: zsin ( f: z -- sin[z] ) i* zsinh (-i)* ;
: ztan ( f: z -- tan[z] ) i* ztanh (-i)* ;
\ Complex inverse trigonometric functions
\ -- after Abramowitz & Stegun, p. 80-81
\ the following is a primitive ANS-compatible data hiding mechanism
: compile! ( xt -- ) COMPILE, ; IMMEDIATE
:noname ( A|B) ( f: x y n -- a) FROT F+ |z| F2/ ; TO noname
: alpha.beta ( f: x y -- alpha beta)
zdup f1.0 [ noname ] COMPILE!
f-rot f1.0 FNEGATE [ noname ] COMPILE!
zdup F+ f-rot F- ;
' NOOP TO noname \ forget hidden xt
\ Note: the following functions have not yet been fully tested. Use
\ with caution! September 17th, 1999 - 18:14
: zasin alpha.beta FASIN FSWAP FDUP f^2 f1.0 F- FSQRT F+ FLN ;
: zacos alpha.beta FACOS FSWAP FDUP f^2 f1.0 F-
FSQRT F+ FLN FNEGATE ;
: zatan zdup FOVER |z|^2 FNEGATE f1.0 F+ F/ F2* FATAN F2/ f-rot
FSWAP f^2 FOVER f1.0 F+ f^2 ( f: -- re y x^2 [y+1]^2 )
F+ FDUP FROT F2* F- F/ FLN F2/ F2/ ;
: zasinh i* zasin (-i)* ;
: zacosh zacos i* ;
: zatanh i* zatan (-i)* ;
\ ------------------------------------------ for use with ftran2xx.f
: zvariable CREATE 2 FLOATS ALLOT ;
: cmplx ( f: x 0 y 0 -- x y) FDROP FNIP ; \ for literals
: cmplxfn ( f: -- ) ; \ for functions or expressions
\ usage: cmplxfn(a*b,atan2(x,y))