(* * LANGUAGE : ANS Forth * PROJECT : Forth Environments * DESCRIPTION : Simulation tools; Complex numbers * CATEGORY : Tool * AUTHOR : Marcel Hendrix * AUTHOR : parts from http://galileo.phys.virginia.edu/classes/551.jvn.fall01/complex_williams.htm * LAST CHANGE : Friday, February 20, 2004 7:20 PM, Marcel Hendrix, bug in POLAR> * LAST CHANGE : Friday, February 21, 2003 7:30 PM, Marcel Hendrix, bug in FEXCHANGE on 'deep' stack? * LAST CHANGE : Thursday, February 20, 2003 8:57 PM, Marcel Hendrix added I0(z), I1(z) * LAST CHANGE : July 31, 2000, Marcel Hendrix added Z!+ Z+!+ * LAST CHANGE : July 18, 1995, Marcel Hendrix *) NEEDS -miscutil NEEDS -blas REVISION -cplx_fsl "ÄÄÄ COMPLEX numbers Version 1.06 ÄÄÄ" PRIVATES DOC COMPLEX (* A definition of a complex number type and associated operations on it. 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 specific platforms or compilers. *) ENDDOC 2 DFLOATS =: =ZFLOAT : ZFLOATS ( u1 -- u2 ) =ZFLOAT * ; : COMPLEX[] ( addr1 index -- addr2 ) XFLOAT[] ; : []COMPLEX ( addr1 index -- addr2 ) []XFLOAT ; 1 FLOATS 1 DFLOATS = [IF] : Z@ F@+ F@ ; ( addr -- ) ( F: -- re im ) : Z@+ F@+ F@+ ; ( addr1 -- addr2 ) ( F: -- re im ) : Z! DUP FLOAT+ F! F! ; ( addr -- ) ( F: re im -- ) : Z!+ FSWAP F!+ F!+ ; ( addr1 -- addr2 ) ( F: re im -- ) : Z0! DUP F0! FLOAT+ F0! ; ( addr -- ) [ELSE] : Z@ DF@+ DF@ ; ( addr -- ) ( F: -- re im ) : Z@+ DF@+ DF@+ ; ( addr1 -- addr2 ) ( F: -- re im ) : Z! DUP DFLOAT+ DF! DF! ; ( addr -- ) ( F: re im -- ) : Z!+ FSWAP DF!+ DF!+ ; ( addr1 -- addr2 ) ( F: re im -- ) : Z0! DUP DF0! DFLOAT+ DF0! ; ( addr -- ) [THEN] \ We need these for the methods array in fsl_util.frt : _Z! ( addr -- ) ( F: re im -- ) Z! ; : _Z@ ( addr -- ) ( F: -- re im ) Z@ ; 1 [IF] : ZCOPY ( a1 incx a2 incy n -- ) >S 2>R 2>R S> 2R> 2R> 5 BLAS-LIB: izcopy FOREIGN DROP ; [ELSE] : ZCOPY ( a1 incx a2 incy n -- ) 0 0 =ZFLOAT LOCALS| sz incx*sz incy*sz n incy a2 incx a1 | incx 1 = incy 1 = AND IF a1 a2 n sz * MOVE EXIT ENDIF incx sz * TO incx*sz incy sz * TO incy*sz incx 1 = IF a1 a2 n incy*sz * BOUNDS ?DO Z@+ I Z! incy*sz +LOOP DROP EXIT ENDIF incy 1 = IF a2 a1 n incx*sz * BOUNDS ?DO I Z@ Z!+ incx*sz +LOOP DROP EXIT ENDIF a1 a2 n incy*sz * BOUNDS ?DO DUP Z@ I Z! incx*sz + incy*sz +LOOP DROP ; [THEN] : ZDUP F2DUP ; ( F: re im -- re im re im ) : ZOVER F2OVER ; ( F: z1 z2 -- z1 z2 z1 ) : ZROT F2ROT ; ( F: z1 z2 z3 -- z2 z3 z1 ) : ZSWAP F2SWAP ; ( F: z1 z2 -- z2 z1 ) : ZDROP F2DROP ; ( F: re im -- ) : ZNIP ZSWAP ZDROP ; ( F: z1 z2 -- z2 ) : Z2DUP ZOVER ZOVER ; ( F: z1 z2 -- z1 z2 z1 z2 ) : ZTUCK ZSWAP ZOVER ; ( F: z1 z2 -- z2 z1 z2 ) : REAL FDROP ; ( F: re im -- re ) : IMAG FNIP ; ( F: re im -- im ) : ZCONJUGATE FNEGATE ; ( F: re im -- re -im ) : ZNEGATE FNEGATE FSWAP ( F: re im -- -re -im ) FNEGATE FSWAP ; 1 [IF] : ZCOPY' ( a1 incx a2 incy n -- ) >S 2>R 2>R S> 2R> 2R> 5 BLAS-LIB: izcopycj FOREIGN DROP ; [ELSE] : ZCOPY' ( a1 incx a2 incy n -- ) 0 0 =ZFLOAT LOCALS| sz incx*sz incy*sz n incy a2 incx a1 | incx sz * TO incx*sz incy sz * TO incy*sz incx 1 = IF a1 a2 n incy*sz * BOUNDS ?DO Z@+ ZCONJUGATE I Z! incy*sz +LOOP DROP EXIT ENDIF incy 1 = IF a2 a1 n incx*sz * BOUNDS ?DO I Z@ ZCONJUGATE Z!+ incx*sz +LOOP DROP EXIT ENDIF a1 a2 n incy*sz * BOUNDS ?DO DUP Z@ ZCONJUGATE I Z! incx*sz + incy*sz +LOOP DROP ; [THEN] : FEXP(jX) FSINCOS FSWAP ; ( F: phi -- cos[phi] +sin[phi] ) : FEXP(-jX) FNEGATE FEXP(jX) ; ( F: phi -- cos[phi] -sin[phi] ) -- convert float to complex 0e FCONSTANT REAL>Z -- convert float to pure imaginary complex value : IMAG>Z 0e FSWAP ; ( F: r -- z ) : Z->R,I ; ( F: x y -- re im ) : Z->I,R FSWAP ; ( F: x y -- im re ) : R,I->Z ; ( F: x y -- re im ) : I,R->Z FSWAP ; ( F: x y -- im re ) -- calculate (x1+j*y1) * (x2+j*y2) : Z* ( F: z1 z2 -- z3 ) 3 FPICK 2 FPICK F* ( x1 y1 x2 y2 x1*x2 ) 3 FPICK 2 FPICK F* F- ( x1 y1 x2 y2 x1*x2-y1*y2 ) 4 FEXCHANGE ( x1*x2-y1*y2 y1 x2 y2 x1 ) F* -FROT ( x1*x2-y1*y2 x1*y2 y1 x2 ) F* F+ ; ( x1*x2-y1*y2 x1*y2+y1*x2 ) : Z2/ F2/ FSWAP F2/ FSWAP ; ( F: re im -- re/2 im/2 ) : Z2* F2* FSWAP F2* FSWAP ; ( F: re im -- re*2 im*2 ) : ZSCALE FROT FOVER F* -FROT F* ; ( F: re1 im1 scale -- re2 im2 ) : ZABS^2 FSQR FSWAP FSQR F+ ; ( F: z -- r^2 ) : ZABS FSQR FSWAP FSQR F+ FSQRT ; ( F: z -- r ) : ZARG FSWAP FATAN2 ; ( F: z -- theta ) : Z= FROT F= F= AND ; ( F: z1 z2 -- ) ( -- bool ) : Z< ZSWAP ZABS -FROT ZABS F< ; ( F: z1 z2 -- ) ( -- bool ) : Z<= ZSWAP ZABS -FROT ZABS F<= ; ( F: z1 z2 -- ) ( -- bool ) : Z> ZSWAP ZABS -FROT ZABS F> ; ( F: z1 z2 -- ) ( -- bool ) : Z>= ZSWAP ZABS -FROT ZABS F>= ; ( F: z1 z2 -- ) ( -- bool ) : ZVARIABLE ( -- ) CREATE FALIGN =ZFLOAT ALLOT DOES> ; : ZCONSTANT ( F: re im -- ) CREATE FSWAP DF, DF, \ Note: assumes that DF, does an FALIGN ! DOES> Z@ ; -- useful complex constants 1e 0e ZCONSTANT 1+0i 0e 0e ZCONSTANT 0+0i 0e 1e ZCONSTANT 0+1i 0e -1e ZCONSTANT 0-1i 1e 1e ZCONSTANT 1+1i 1e -1e ZCONSTANT 1-1i : Z+ ( F: z1 z2 -- z3 ) FROT F+ -FROT F+ FSWAP ; : Z+! ( addr -- ) ( F: re im -- ) DUP Z@ Z+ Z! ; : Z+!+ ( addr1 -- addr2 ) ( F: re im -- ) DUP Z@ Z+ Z!+ ; : Z- ( F: z1 z2 -- z3 ) FROT FSWAP F- -FROT F- FSWAP ; : Z-! ( addr -- ) ( F: re im -- ) FNEGATE FSWAP FNEGATE FSWAP Z+! ; : ZVALUE ( F: x y "str" -- ) CREATE IMMEDIATE FALIGN HERE =ZFLOAT ALLOT Z! DOES> FALIGNED %VAR @ 0 %VAR ! CASE ( -to) -2 OF ALITERAL EVAL" Z-! " ENDOF ( +to) -1 OF ALITERAL EVAL" Z+! " ENDOF ( from) 0 OF ALITERAL EVAL" Z@ " ENDOF ( to) 1 OF ALITERAL EVAL" Z! " ENDOF ( 0to) 2 OF ALITERAL EVAL" Z0! " ENDOF ( 'of) 3 OF ALITERAL ENDOF ( sizeof) 4 OF DROP 2 DFLOATS ILITERAL ENDOF ( /of) 5 OF DROP 1 ILITERAL ENDOF DUP #-21 ?ERROR" ZVALUE: undefined message" ENDCASE ; : POLAR> ( F: r t -- x y ) \ convert from polar to Cartesian form FSINCOS FROT FSWAP FOVER F* -FROT F* ; : >POLAR ( F: x y -- r t ) \ convert from Cartesian to polar form F2DUP ZABS -FROT FSWAP FATAN2 ; : Z/ ( F: z1 z2 -- z3 ) ZDUP FSQR FSWAP FSQR F+ 1/F FLOCALS| 1/scale t-imag t-real | FOVER t-real F* FOVER t-imag F* F+ 1/scale F* FROT t-imag F* FNEGATE FROT t-real F* F+ 1/scale F* ; -- Extra : ZLN >POLAR FSWAP FLN FSWAP ; ( F: z -- ln[|z|]+iarg[z] ) : ZEXP FSINCOS FSWAP FROT FEXP ZSCALE ; ( F: z -- exp[z] ) : ZSQR ( re im -- r i ) FLOCALS| im(c) re(c) | re(c) FSQR im(c) FSQR F- re(c) im(c) F* F2* ; -- Complex Square Root : ZSQRT ( re im -- r i ) F2/ FLOCAL b F2/ FLOCAL a a FSQR b FSQR F+ FSQRT FLOCAL w a F0> IF a w F+ FSQRT FDUP TO w F0= IF 0+0i ELSE w b w F/ ENDIF ELSE w a F- FSQRT b F0< IF FNEGATE ENDIF FDUP TO w F0= IF 0+0i ELSE b w F/ w ENDIF ENDIF ; : 1/Z ( F: z1 -- 1/z ) 1+0i ZSWAP Z/ ; : Z^N ( F: z -- z^n ) ( n -- ) 1+0i ZSWAP BEGIN DUP 0> WHILE DUP 1 AND IF ZTUCK Z* ZSWAP ENDIF ZSQR 2/ REPEAT DROP ZDROP ; : Z** ( F: x y u v -- [x+iy]^[u+iv] ) ZSWAP ZLN Z* ZEXP ; CREATE 1/fac^2 PRIVATE #10 FLOATS ALLOT CREATE 1/facgam PRIVATE #10 FLOATS ALLOT 0+0i ZVALUE zzz PRIVATE : do-fac ( n -- ) ( F: -- r ) 1 SWAP 0 ?DO I 1+ * LOOP S>F FSQR 1/F ; : fac! #11 1 DO I do-fac 1/fac^2 I 1- FLOAT[] F! LOOP #11 1 DO I do-fac I 1+ S>F F/ 1/facgam I 1- FLOAT[] F! LOOP ; fac! FORGET do-fac -- Modified Bessel functions of 0 and 1st kind, for any complex arguments < 5. -- Not very accurate: error < 1e-7 [when |z| < 5] : I0(z) ( F: z -- z2 ) ZSQR 0.25e ZSCALE ZDUP TO zzz 1/fac^2 1+0i #10 0 DO ZOVER F@+ ZSCALE Z+ ZSWAP zzz Z* ZSWAP LOOP DROP ZNIP ; : I1(z) ( F: z -- z2 ) ZDUP ZSQR 0.25e ZSCALE ZDUP TO zzz 1/facgam 1+0i #10 0 DO ZOVER F@+ ZSCALE Z+ ZSWAP zzz Z* ZSWAP LOOP DROP ZNIP Z* 0.5e ZSCALE ; -- Bessel functions of 0 and 1st kind, for any complex arguments < 5. -- Not very accurate: error < 1e-7 [when |z| < 5] : J0(z) ( F: z -- z2 ) ZSQR -0.25e ZSCALE ZDUP TO zzz 1/fac^2 1+0i #10 0 DO ZOVER F@+ ZSCALE Z+ ZSWAP zzz Z* ZSWAP LOOP DROP ZNIP ; : J1(z) ( F: z -- z2 ) ZDUP ZSQR -0.25e ZSCALE ZDUP TO zzz 1/facgam 1+0i #10 0 DO ZOVER F@+ ZSCALE Z+ ZSWAP zzz Z* ZSWAP LOOP DROP ZNIP Z* 0.5e ZSCALE ; -- Extra trigonometric functions : I* FNEGATE FSWAP ; PRIVATE : (-I)* FSWAP FNEGATE ; PRIVATE : ZSINH ZEXP ZDUP 1/Z Z- Z2/ ; : ZCOSH ZEXP ZDUP 1/Z Z+ Z2/ ; : ZTANH ZEXP ZSQR I* ZDUP F1- ZSWAP F1+ Z/ ; : ZCOTH ZTANH 1/Z ; : ZSIN I* ZSINH (-I)* ; : ZCOS I* ZCOSH ; : ZTAN I* ZTANH (-I)* ; : ZCOT I* ZCOTH I* ; -- Complex inverse trigonometric functions \ In the following, we use phrases like "1E x+" instead of \ "z=1 z+", for stability on branch cuts involving signed zero. \ This follows a suggestion by Kahan [2], and it actually makes \ a difference in every one of the functions. \ To avoid unneeded calculations on the other part that could \ raise gratuitous overflow or underflow signals and changes in \ the sign of zero (Kahan) : x+1 ( F: x y -- x+a y ) FSWAP F1+ FSWAP ; PRIVATE : x-1 ( F: x y -- x-a y ) FSWAP F1- FSWAP ; PRIVATE : parg ( F: z -- t ) FSWAP FATAN2 ; PRIVATE : pln ( F: z -- ln[z].prin ) ZDUP parg -FROT ZABS FLN FSWAP ; PRIVATE : psqrt ( F: z1 -- z2 ) ZSQRT ; PRIVATE : ZASINH ( F: z -- ln[z+sqrt[[z+i][z-i]] ) ZDUP F1+ ZOVER F1- Z* psqrt Z+ pln ; : ZACOSH ( F: z -- 2ln[sqrt[[z+1]/2]+sqrt[[z-1]/2] ) ZDUP x-1 Z2/ psqrt ZSWAP x+1 Z2/ psqrt Z+ pln Z2* ; : ZATANH ( z -- [ln[1+z]-ln[1-z]]/2 ) ZDUP x+1 pln ZSWAP x-1 ZNEGATE pln Z- Z2/ ; : ZACOTH ( z -- [ln[-1-z]-ln[1-z]]/2 ) ZNEGATE ZDUP x-1 pln ZSWAP x+1 pln Z- Z2/ ; : ZASIN I* ZASINH (-I)* ; ( z -- -iln[iz+sqrt[1-z^2]] ) : ZACOS PI/2 0e ZSWAP ZASIN Z- ; ( z -- pi/2-asin[z] ) : ZATAN I* ZATANH (-I)* ; ( z -- [ln[1+iz]-ln[1-iz]]/2i ) : ZACOT (-I)* ZACOTH (-I)* ; ( z -- [ln[[z+i]/[z-i]]/2i ) : +ZE.R ( n -- ) ( F: re im -- ) FSWAP ." (" SPACE DUP +E.R SPACE +E.R ." )" ; : +ZCE.R ( n -- ) ( F: re im -- ) >S FROUND F>S FROUND F>S ." (" SPACE S .R SPACE S> .R ." )" ; DEPRIVE NESTING @ 1 = [IF] : Z. ( F: re im -- ) 6 +ZE.R ; 0+0i zvalue c1 ( define complex value c1 ) 1e 0e TO c1 ( store 1+i0 in c1 ) 1+0i TO c1 ( equivalent to previous line ) CR c1 z. .( -> 1 + i 0 ) CR 2e 3e TO c1 .( -> store 2+i3 in c1 ) CR c1 zconjugate z. .( -> 2 - i 3 ) CR c1 znegate z. .( -> -2 - i 3 ) CR c1 0+1i z* z. .( -> -3 + i 2 ) CR c1 0-1i z* z. .( -> 3 - i 2 ) CR c1 zsqr z. .( -> -5 + i 12 ) CR c1 zsqrt z. .( -> 1.67415 + i 0.895977 ) CR c1 zabs^2 f. .( -> 13 ) CR c1 zabs f. .( -> 3.60555 ) CR c1 7 z^n z. .( -> 6554 + i 4449 ) CR c1 7e 0e z** z. .( -> 6554 + i 4449 ) CR c1 2e zscale z. .( -> 4 + i 6 ) CR c1 2e 1/f zscale z. .( -> 1 + i 1.5 ) CR c1 1/z z. .( -> 0.153846 - i 0.230769 ) CR c1 zarg f. .( -> 0.982794 ) CR c1 >polar f. f. .( -> 0.982794 3.60555 ) CR c1 >polar polar> z. .( -> 2 + i 3 ) CR c1 zexp z. .( -> -7.31511 + i 1.04274 ) CR c1 zln z. .( -> 1.28247 + i 0.982794 ) 0+0i zvalue c2 2e 3e TO c1 1e 5e TO c2 CR c1 c2 z+ z. .( -> 3 + i 8 ) CR c1 c2 z- z. .( -> 1 - i 2 ) CR c1 c2 z* z. .( -> -13 + i 13 ) CR c1 c2 z/ z. .( -> 0.653846 - i 0.269231 ) FORGET Z. [THEN] (* End of File *)