\ 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 angles 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.3 \ Last revision: April 24, 2005 \ Version 0.8.3 \ 3Mar05 * Revised ZCOSH for fairly good accuracy with \ conservation of signed zero on the real axis. \ 4Mar05 * Rewrote ZACOS to use minimal computation. \ 24Apr05 * Renamed (-I)* as -I*. Removed simple floating point \ and complex constants, leaving them to the user. \ * Added Y+ and Y-, syntactic sugar for F+ and F-. \ 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 optional versions of ZSINH and \ ZTANH which respect the sign of zero, and of ZCOSH which \ respects it on the real axis. 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. \ In fact, with the optional versions mentioned above, the \ functions here pass all tests in complex-ieee-test.fs, except \ for signed zero in the imaginary part at points of real \ analyticity in ZASIN and ZASINH. The more aggressive Kahan \ algorithms in complex-kahan.fs and pfe pass those as well. s" FLOATING-EXT" environment? [IF] ( flag) drop s" FLOATING-STACK" environment? [IF] ( maxdepth) drop DECIMAL \ important for fp input 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] 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. ; : 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. 1E 0E 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 is reasonably accurate and 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] 0 [IF] : zcosh zexp zdup 1/z z+ z2/ ; [ELSE] \ This version is reasonably accurate and preserves \ signed zero on the real axis. : zcosh FSINCOS [ noname ] LITERAL ( f: cos[y]) F! [ noname float+ ] LITERAL ( f: sin[y]) F! FDUP FCOSH [ noname ] LITERAL F@ F* ( f: x ch[x]cos[y]) FSWAP FSINH [ noname float+ ] LITERAL F@ F* ( f: ch[x]cos[y] sh[x]sin[y]) ; [THEN] 0 [IF] : ztanh zexp z^2 i* zdup 1E F- zswap 1E 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 \ "1E 0E 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 for signed zero 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] ) zasin znegate pi/2 x+ ; : 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]