\ 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 ; : ztuck zswap zover ; : z= ( f: x y u v -- s: flag ) FROT F= F= AND ; : 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+ ( f: a b x y -- a+x b+y) FROT F+ f-rot F+ FSWAP ; : znegate FSWAP FNEGATE FSWAP FNEGATE ; : z- znegate z+ ; \ 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+ ( f: x y a -- x+a y ) frot f+ fswap ; : x- ( f: x y a -- x-a y ) fnegate frot f+ fswap ; \ : y+ ( f: x y a -- x y+a ) f+ ; \ : y- ( f: x y a -- x y-a ) f- ; : |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 FDUP F0= IF FDROP zdrop 0E ELSE f-rot FMIN ( f: max min) FOVER F/ f^2 1E F+ FSQRT F* THEN ; : 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] ) FSWAP fatan2 ; : >polar ( f: x+iy -- r theta ) zdup |z| f-rot arg ; : polar> ( f: r theta -- x+iy ) FSINCOS FROT z*f FSWAP ; : i* FNEGATE FSWAP ; ( f: x+iy -- -y+ix) : (-i)* FSWAP FNEGATE ; ( f: x+iy -- y-ix) : pln ( f: z -- ln[z].prin ) zdup parg f-rot |z| FLN FSWAP ; : zln ( f: z -- ln[|z|]+iarg[z] ) >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 \ Use Z^ instead for n > 50 or so. z=1 zswap BEGIN DUP 0> WHILE DUP 1 AND IF ztuck z* zswap THEN z^2 2/ REPEAT zdrop DROP ; : z^ ( f: x y u v -- [x+iy]^[u+iv] ) zswap zln z* zexp ; : psqrt ( f: z -- sqrt[z].prin ) ( Kahan's algorithm without overflow/underflow avoidance and without treating infinity. But it should handle signed zero properly. ) zdup [ noname FLOAT+ ] LITERAL ( f: y) F! [ noname ] LITERAL ( f: x) F! |z| [ noname ] LITERAL F@ fabs f+ f2/ fsqrt ( f: rho=sqrt[[|z|+|x|]/2]) fdup f0= [ noname FLOAT+ ] LITERAL F@ ( f: rho y) 0= IF \ rho <> 0 fover f/ f2/ ( f: rho eta=[y/rho]/2) [ noname ] LITERAL F@ f0< IF \ x < 0 fabs fswap ( f: |eta| rho) [ noname FLOAT+ ] LITERAL F@ FDUP F0< -0e 0e F~ OR IF \ y < 0 FNEGATE THEN ( f: |eta| cps[rho,y]) THEN THEN ; PRINCIPAL-ARG [IF] : zsqrt psqrt ; [ELSE] : 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) [THEN] \ Complex trigonometric functions \ All stack patterns are ( f: z -- func[z] ). 0 [IF] : zsinh zexp zdup 1/z z- z2/ ; [ELSE] \ This version preserves signed zero. : zsinh FSINCOS [ noname ] LITERAL ( f: cos[y]) F! [ noname FLOAT+ ] LITERAL ( f: sin[y]) F! FDUP FSINH [ noname ] LITERAL F@ F* ( f: x sh[x]cos[y]) FSWAP FCOSH [ noname FLOAT+ ] LITERAL F@ F* ( f: sh[x]cos[y] ch[x]sin[y]) ; [THEN] : zcosh zexp zdup 1/z z+ z2/ ; 0 [IF] : ztanh zexp z^2 i* zdup f1.0 F- zswap f1.0 F+ z/ ; [ELSE] \ This version, based on Kahan, preserves signed zero. : ztanh ( [1 + tan^2[y]] cosh[x] sinh[x] + i tan[y] tanh[z] = ----------------------------------------- 1 + [1 + tan^2[y]] sinh^2[x] ) FTAN FDUP f^2 1E F+ ( f: x t=tan[y] b=1+t^2) [ noname ] LITERAL F! FSWAP FDUP FSINH ( f: t x sh[x]) [ noname FLOAT+ ] LITERAL F! FCOSH [ noname FLOAT+ ] LITERAL F@ ( f: t ch[x] sh[x]) F* [ noname ] LITERAL F@ F* FSWAP ( f: c=ch[x]sh[x]b t) [ noname ] LITERAL F@ [ noname FLOAT+ ] LITERAL F@ f^2 F* 1E F+ ( f: c t 1+b*sh^2[x]) z/f ; [THEN] : 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. : zasinh ( f: z -- ln[z+sqrt[[z+i][z-i]] ) \ This is more stable than the version with z^2+1. zdup 1E F+ zover 1E F- z* psqrt z+ pln ; : zacosh ( f: z -- 2ln[sqrt[[z+1]/2]+sqrt[[z-1]/2] ) zdup 1E x- z2/ psqrt zswap 1E x+ z2/ psqrt z+ pln z2* ; : zatanh ( f: z -- [ln[1+z]-ln[1-z]]/2 ) zdup 1E x+ pln zswap 1E x- znegate pln z- z2/ ; : zacoth ( f: z -- [ln[-1-z]-ln[1-z]]/2 ) znegate zdup 1E x- pln zswap 1E x+ pln z- z2/ ; : zasin ( f: z -- -iln[iz+sqrt[1-z^2]] ) i* zasinh (-i)* ; : zacos ( f: z -- pi/2-asin[z] ) pi/2 0E zswap zasin z- ; : zatan ( f: z -- [ln[1+iz]-ln[1-iz]]/2i ) i* zatanh (-i)* ; : zacot ( f: z -- [ln[[z+i]/[z-i]]/2i ) (-i)* zacoth (-i)* ; \ --------------------------------- for use with ftran2xx.f : cmplx ( f: x 0 y 0 -- x y) FDROP fnip ; \ --------------------------------- REFERENCES ( 1. M. Abramowitz and I. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, US Government Printing Office, 10th Printing December 1972, Secs. 4.4, 4.6. 2. William Kahan, "Branch cuts for complex elementary functions", The State of the Art in Numerical Analysis, A. Iserles and M.J.D. Powell, eds., Clarendon Press, Oxford, 1987, pp. 165-211. 3. Robert M. Corless, James H. Davenport, David J. Jeffrey, Stephen M. Watt, "'According to Abramowitz and Stegun' or arcoth needn't be uncouth", ACM SIGSAM Bulletin, June, 2000, pp. 58-65. ) [ELSE] .( ***Separate floating point stack not available.) [THEN] [ELSE] .( ***Floating point not available.) [THEN]