complex.f
\ ANS Forth Complex Arithmetic Lexicon
\ --------------------------------------------------------------
\ Copyright (C) 1998 Julian V. Noble \
\ \
\ This library is free software; you can redistribute it \
\ and/or modify it under the terms of the GNU Lesser General \
\ Public License as published by the Free Software Foundation; \
\ either version 2.1 of the License, or at your option any \
\ later version. \
\ \
\ This library is distributed in the hope that it will be \
\ useful, but WITHOUT ANY WARRANTY; without even the implied \
\ warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR \
\ PURPOSE. See the GNU Lesser General Public License for more \
\ details. \
\ \
\ You should have received a copy of the GNU Lesser General \
\ Public License along with this library; if not, write to the \
\ Free Software Foundation, Inc., 59 Temple Place, Suite 330, \
\ Boston, MA 02111-1307 USA. \
\ --------------------------------------------------------------
\ 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.
\ Modifications: David N. Williams
\ Version: 0.8.2
\ Last revision: March 5, 2003
\ Version 0.8.2
\ 5Mar03 * Release with changes under 0.8.1 Passes all our
\ tests, except for a few "exotic" signed zero
\ properties.
\ Version 0.8.1
\ 21Feb03 * Rewrote PSQRT for stability on the branch cut.
\ * Added X+ and X-. Used to make ZASINH, ZACOSH,
\ ZATANH, and ZACOTH valid on their cuts, labeled by
\ signed zero.
\ * Passed complex-test.fs on MacOS X, including signed
\ zero and branch cut tests.
\ 22Feb03 * Replaced Z. and ZS. with jvn's code for Krishna
\ Myneni's suggestion to factor out the sign of the
\ imaginary part.
\ Version 0.8.0
\ 15Dec02 * Started revision of jvn's complex.f.
\ 20Feb03 * Release.
\ The basic modifications have been a few changes in
\ floating-point alignment and the completion of the definitions
\ of the inverse functions.
\ The definitions here coincide with Abramowitz and Stegun [1],
\ Kahan [2], and the OpenMath standard [3], which produce
\ principal branches given principal branches for square roots
\ and natural logs. The formulas, or "principal expressions",
\ are selected from [3], with choices among equivalent, formally
\ correct possibilities based mainly on computational
\ conciseness, with a nod to numerical stability where the
\ authors mention it. Those authors do not claim to analyze
\ numerical stability, but Kahan does, and we implement his
\ algorithms in Forth in a separate file.
\ The original Noble code uses the convention common among
\ physicists for branch cuts, with arguments between zero and
\ 2pi, especially for logs and noninteger powers. The numerical
\ analysis community is pretty unanimous about using principal
\ branches instead.
\ Everybody seems to agree on the nontriviality of the branch
\ cuts for the inverse functions and to follow Abramowitz and
\ Stegun, who define them in terms of principal branches.
\ In this code we include a PRINCIPAL-ARG switch to select
\ between the two common conventions for arg, log, and
\ nonintegral powers, but we use only principal arguments for
\ the inverse functions.
\ Kahan pays attention to signed zero, where available in IEEE
\ 754/854 implementations. We address that a couple of ways in
\ this file. One is to provide uncommentable versions of ZSINH
\ and ZTANH which respect the sign of zero. The other is to
\ write the functions having branch cuts so that signed zero in
\ the appropriate x or y input produces correct values on the
\ cuts.
s" FLOATING-EXT" environment? [IF] ( flag) drop
s" FLOATING-STACK" environment? [IF] ( maxdepth) drop
\ MARKER -complex
1.570796326794896619231E FCONSTANT pi/2
6.283185307179586476925E FCONSTANT 2pi
\ The PRINCIPAL-ARG flag controls conditional compilation
\ with/without the principal argument for FATAN2, ARG, >POLAR,
\ ZSQRT, ZLN, and Z^. The inverse functions are always defined
\ with principal arguments, in accord with Abramowitz and
\ Stegun [1], Kahan [2], and OpenMath [3].
false CONSTANT PRINCIPAL-ARG \ output 0 <= arg < 2pi
\ true CONSTANT PRINCIPAL-ARG \ output -pi < arg <= pi
IMMEDIATE
\ In ANS Forth FATAN2 is ambiguous. The following assumes it
\ either gives -pi < arg < pi (or maybe <=) or 0 <= arg < 2pi.
\ Then PARG is defined to force the principal arg for later use
\ in the forced principal arg function PLN, and FATAN2 is
\ redefined according to the PRINCIPAL-ARG flag for later use in
\ defining ARG.
-1E 0E ( f: y x) FATAN2 F0<
( arg<0) DUP
[IF] \ FATAN2 principal
: parg ( f: x y -- princ.arg ) FSWAP FATAN2 ;
[ELSE] \ FATAN2 not principal arg
: parg ( f: x y -- princ.arg )
FDUP F0< FSWAP FATAN2
( y<0) IF 2pi F- THEN ;
[THEN]
PRINCIPAL-ARG [IF]
( arg<0) 0= [IF] \ FATAN2 not principal, redefine
: fatan2 FOVER ( f: y x y) F0< FATAN2
( y<0) IF 2pi F- THEN ;
[THEN]
[ELSE]
( arg<0) [IF] \ FATAN2 principal, redefine
: fatan2 FOVER ( f: y x y) F0< FATAN2
( y<0) IF 2pi F+ THEN ;
[THEN]
[THEN]
s" [undefined]" pad c! pad char+ pad c@ move
pad find nip 0=
[IF]
\ Leave true if name is in the search order, else leave false.
: [undefined] ( "name" -- flag )
bl word find nip 0= ; immediate
[THEN]
\ non-Standard fp words jvn has found useful
[undefined] f0.0 [IF] 0.0E FCONSTANT f0.0 [THEN]
[undefined] f1.0 [IF] 1.0E FCONSTANT f1.0 [THEN]
[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]
\ added by dnw
[undefined] f2* [IF] : f2* FDUP F+ ; [THEN]
[undefined] f2/ [IF] : f2/ 0.5E F* ; [THEN]
\ --------------------------------- LOAD, STORE
: z@ DUP F@ FLOAT+ F@ ; ( adr -- f: -- z)
: z! DUP FLOAT+ F! F! ; ( adr -- f: z --)
: zvariable FALIGN HERE 2 FLOATS ALLOT CONSTANT ;
: zconstant ( "name" f: z -- )
FALIGN HERE 2 FLOATS ALLOT DUP z! CREATE ,
DOES> ( f: -- z ) @ z@ ;
\ --------------------------------- MANIPULATE FPSTACK
: z. FSWAP F. ( F: x y -- ) \ emit complex #
FDUP F0<
DUP INVERT [CHAR] + AND SWAP [CHAR] - AND +
EMIT ." i " FABS F. ;
: zs. FSWAP FS. ( F: x y -- ) \ emit complex #, scientific notation
FDUP F0<
DUP INVERT [CHAR] + AND SWAP [CHAR] - AND +
EMIT ." i " FABS FS. ;
: 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 3 FLOATS ALLOT VALUE noname
: 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
: znip zswap zdrop