( Title: Complex arithmetic lexicon with mpfr reliable floating point Ported from: Forth Scientific Library Algorithm #60 File: complex-rf.fs Version: 1.0.6 License: LGPL 3 Test file: complex-rf-test.fs Revised: January 13, 2012 ) \ Original author: \ Copyright (C) 1998 Julian V. Noble \ Modifications not derived from the original: \ Copyright (C) 2002, 2003, 2005, 2008-2012 David N. Williams ( 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 3 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. If you take advantage of the option in the LGPL to put a par- ticular version of this library under the GPL, the author[s] would regard it as polite if you would put any direct modifi- cations under the LGPL as well, and include a copy of this request near the beginning of the modified library source. A "direct modification" is one that enhances or extends the library in line with its original concept, as opposed to devel- oping a distinct application or library which might use it. This library ports the FSL complex arithmetic lexicon to use reliable floating point via the multiprecision forth-gmpfr library, which binds to GNU gmp and mpfr. Unattributed changes are by dnw. The revision date above may reflect cosmetic changes not logged here. Version 1.0.6 9Jan12 * Fixed ZCONSTANT to refer to LAST-ZCONSTANTP instead of LAST-RFCONSTANTP, and added missing DUP to FREE-ZCONSTANTS. 10Jan12 * Fixed ZCONSTANT to prevent multiple freeing of mfpr objects. 12Jan12 * Changed LGPL 2.1 to LGPL 3. 13Jan12 * Replaced FLOAT+ by CELL+ in Z*, complex-rf.fs. Thanks to Marcel Hendrix. Version 1.0.5 9May11 * To make porting to reliable floats easier for clients of FSL complex arithmetic, reverted to "f" nomen- clature, instead of "rf", for mixed complex and real operations. Also reverted to "znegate" instead of "-z". * Renamed this file from rfcomplex.fs to complex-rf.fs. 11May11 * Added RELIABLE-FLOATS selection flag. Version 1.0.4 6May11 * Started from FSL Algorithm #60, complex.fs 1.0.3 and the xf version of its tests, complex-test-xf.fs. * Passed rfcomplex-test.fs in principal argument mode. 7May11 * Passed in nonprincipal argument mode. * Released. ) decimal true constant RELIABLE-FLOATS immediate s" gmpfr.fs" included \ *** TABLE OF CONTENTS --- 1 INTRODUCTION 2 EXTRA RF WORDS rf2* rf2/ rf^2 -0rf= 3 VALUES AND FCONSTANTS PRINCIPAL-ARG pi pi/2 2pi 4 LOAD AND STORE COMPLEX COMPLEXES COMPLEX+ z@ z! ZVARIABLE ZCONSTANT 5 OUTPUT .f[bl|-] .f[+|-]i zs. 6 RFSTACK z= zdrop zdup zswap zover znip ztuck zrot -zrot 7 ARITHMETIC conjg real imag i* -i* znegate z* z*f z/ z/f z+ z- x+ x- y+ y- |z| 1/z z2* z2/ z^2 z^3 z^4 z^n |z|^2 8 FUNCTIONS parg nparg arg >polar polar> pln psqrt zsqrt zexp zln z^ zsinh zcosh ztanh zcoth zsin zcos ztan zcot 9 INVERSE FUNCTIONS zasinh zacosh zatanh zacoth zasin zacos zatan zacot 10 REFERENCES --- \ *** 1 INTRODUCTION --- The reliable floating point stack is called the "rfstack". Complex numbers x+iy are stored on the rfstack as (rf: -- x y), also written as (rf: -- z ). Angles are in radians. Polar representation measures angles from the positive x-axis. 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 library. 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 with versions of ZSINH and ZTANH which respect the sign of zero, and of ZCOSH which respects it on the real axis. We also write the functions having branch cuts so that signed zero in the appropriate x or y input produces correct values on the cuts. The word COMPLEX may be used with Forth Scientific Library arrays, after loading the system-dependent fsl-util file. Here is an example of an array sufficient for 20 complex numbers: 20 COMPLEX ARRAY a{ rf" 1.0" rf" 0.0" a{ 0 } z! rf" 1.0" rf" 1.0" a{ 1 } z! --- \ *** 2 EXTRA RF WORDS : rf2* (rf: x -- x*2 ) 2 rfu* ; : rf2/ (rf: x -- x/2 ) 2 rfu/ ; : rf^2 (rf: x -- x*x ) 2 rf^u ; : 1/rf (rf: x -- 1/x ) -1 rf^n ; : -0rf= (rf: x -- s: [x=-0] ) is-rf-signbit rf0= and ; \ *** 3 VALUES AND RFCONSTANTS --- The PRINCIPAL-ARG flag controls calculation with/without the principal argument for 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]. --- \ This value may be changed after this file is loaded. false VALUE PRINCIPAL-ARG \ output 0 <= arg < 2pi \ true VALUE PRINCIPAL-ARG \ output -pi < arg <= pi pi>rf rfdup rfconstant pi rfdup rf2/ rfconstant pi/2 rf2* rfconstant 2pi \ *** 4 LOAD AND STORE 2 cells constant COMPLEX : COMPLEXES ( n -- n*complex ) complex * ; : COMPLEX+ ( 'z -- 'z+complex ) complex + ; \ Hidden, *nonnestable* scratch storage for stuff from rfstack. falign here 3 complexes allot value rfnoname rfnoname 3 complexes 0 fill \ initialize for gc : z@ ( addr -- rf: -- z) dup rf@ cell+ rf@ ; : z! ( addr -- rf: z --) dup cell+ rf! rf! ; : ZVARIABLE ( "name" -- ) \ compile ( -- addr ) \ run create 0 , 0 , ; 0 value LAST-ZCONSTANTP : ZCONSTANT ( "name" rf: z -- ) rfs@ rfspace is-grec 0= IF rfspace gs>extern >rfs-copy THEN \ new ext copy rfswap (rf: y x) rfs@ rfspace is-grec 0= IF rfspace gs>extern >rfs-copy THEN \ new ext copy rfspace gs>extern ( 'x) rfspace gs>extern ( 'y) create here >r swap , , last-zconstantp , r> to last-zconstantp DOES> (rf: -- z.ext ) dup @ (>rfs) cell+ @ (>rfs) ; : free-zconstants ( -- ) last-zconstantp BEGIN dup WHILE dup @ mpfr_clear cell+ dup @ mpfr_clear cell+ @ REPEAT drop ; \ *** 5 OUTPUT : .rf[bl|-] (rf: x -- ) \ rfa = |x| rf>rfa rfa-signbit |rfa| dup invert bl and swap [char] - and + emit ; : .rf[+|-]i (rf: y -- ) \ rfa = |y| rf>rfa rfa-signbit |rfa| dup invert [char] + and swap [char] - and + emit ." i" ; \ emit complex # : zs. (rf: z -- ) rfswap .rf[bl|-] rfa. .rf[+|-]i rfa. ; \ *** 6 RFSTACK : zdrop (rf: z --) rfdrop rfdrop ; : zdup (rf: z -- z z ) rfover rfover ; : zswap (rf: z1 z2 -- z2 z1 ) 2 0 rfexchange 3 1 rfexchange ; : zover (rf: z1 z2 -- z1 z2 z1 ) 3 rfpick 3 rfpick ; : znip (rf: z1 z2 -- z2 ) zswap zdrop ; : ztuck (rf: z1 z2 -- z2 z1 z2 ) zswap zover ; : zrot (rf: z1 z2 z3 -- z2 z3 z1 ) [ rfnoname ] LITERAL z! [ rfnoname COMPLEX+ ] LITERAL z! [ rfnoname 2 COMPLEXES + ] LITERAL z! [ rfnoname COMPLEX+ ] LITERAL z@ [ rfnoname ] LITERAL z@ [ rfnoname 2 COMPLEXES + ] LITERAL z@ ; : -zrot ( f: z1 z2 z3 -- z3 z1 z2 ) [ rfnoname ] LITERAL z! [ rfnoname COMPLEX+ ] LITERAL z! [ rfnoname 2 COMPLEXES + ] LITERAL z! [ rfnoname ] LITERAL z@ [ rfnoname 2 COMPLEXES + ] LITERAL z@ [ rfnoname COMPLEX+ ] LITERAL z@ ; \ *** 7 ARITHMETIC : real (rf: x y -- x ) rfdrop ; : imag (rf: x y -- y ) rfnip ; : conjg (rf: x y -- x -y ) -rf ; : znegate (rf: x y -- -x -y ) rfswap -rf rfswap -rf ; : z= (rf: z1 z2 -- s: flag ) rfrot rf= rf= and ; : z*f (rf: x y a -- x*a y*a ) rfrot rfover rf* -rfrot rf* ; : z/f (rf: x y a -- x/a y/a ) 1/rf z*f ; : z+ (rf: z1 z2 -- z1+z2 ) rfrot rf+ -rfrot rf+ rfswap ; : z- (rf: a b x y -- a-x b-y ) ( Kahan says, to conserve signed zero, write -y+b for b-y instead of -[y-b]. ) -rf rfrot rf+ -rfrot rf- rfswap ; : z* (rf: 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 rf+ (rf: x y u v u+v) [ rfnoname ] literal rf! (rf: x y u v) rfover rf- (rf: x y u v-u) [ rfnoname cell+ ] literal rf! (rf: x y u) rfrot rfdup (rf: y u x x) [ rfnoname cell+ ] literal rf@ (rf: y u x x v-u) rf* [ rfnoname cell+ ] literal rf! (rf: y u x) rfrot rfdup (rf: u x y y) [ rfnoname ] literal rf@ (rf: u x y y u+v) rf* [ rfnoname ] literal rf! (rf: u x y) rf+ rf* rfdup (rf: u*[x+y] u*[x+y]) [ rfnoname ] literal rf@ rf- (rf: u*[x+y] x*u-y*v) rfswap [ rfnoname cell+ ] literal rf@ (rf: x*u-y*v u*[x+y] x*[v-u]) rf+ (rf: x*u-y*v x*v+y*u) ; --- These words avoid unneeded calculations on the other part that could raise gratuitous overflow or underflow signals and changes in the sign of zero (Kahan) --- : x+ (rf: x y a -- x+a y ) rfrot rf+ rfswap ; : x- (rf: x y a -- x-a y ) -rf rfrot rf+ rfswap ; : y+ (rf: x y a -- x y+a ) rf+ ; : y- (rf: x y a -- x y-a ) rf- ; : |z|^2 (rf: z -- |z|^2 ) rf^2 rfswap rf^2 rf+ ; : |z| (rf: x y -- |z|) rf-hypot ; : 1/z (rf: z -- 1/z ) ( Writing 1/z this way reduces overflow probability ) -rf zdup |z| 1/rf rfdup [ rfnoname ] LITERAL rf! z*f [ rfnoname ] LITERAL rf@ z*f ; : z/ (rf: z1 z2 -- z1/z2 ) 1/z z* ; : z2/ (rf: z -- z/2 ) rf2/ rfswap rf2/ rfswap ; : z2* (rf: z -- z*2 ) rf2* rfswap rf2* rfswap ; : i* (rf: x+iy -- -y+ix) -rf rfswap ; : -i* (rf: x+iy -- y-ix) rfswap -rf ; \ Raise z to an integer power. : z^2 (rf: z -- z^2 ) zdup z* ; : z^3 (rf: z -- z^3 ) zdup z^2 z* ; : z^4 (rf: z -- z^4 ) z^2 z^2 ; \ Use Z^ instead for n > 50 or so. : z^n ( n -- f: z -- z^n ) 1 u>rf 0>rf zswap BEGIN dup 0> WHILE dup 1 and IF ztuck z* zswap THEN z^2 2/ REPEAT zdrop drop ; \ *** 8 FUNCTIONS : parg (rf: x y -- princ.arg ) rfswap rf-atan2 ; : nparg (rf: x y -- npinc.arg ) rfdup rf0< rfswap rf-atan2 ( y<0) IF 2pi rf+ THEN ; : arg (rf: x y -- arg[x+iy] ) PRINCIPAL-ARG IF parg ELSE nparg THEN ; : >polar (rf: x+iy -- r theta ) zdup |z| -rfrot arg ; : polar> (rf: r theta -- x+iy ) rf-sincos rfrot z*f rfswap ; : pln (rf: z -- ln[z].prin ) zdup parg -rfrot |z| rf-ln rfswap ; : zln (rf: z -- ln[|z|]+iarg[z] ) >polar rfswap rf-ln rfswap ; : zexp (rf: z -- exp[z] ) rf-sincos rfswap rfrot rf-exp z*f ; : z^ (rf: z1 z2 -- [z1]^[z2] ) zswap zln z* zexp ; : psqrt (rf: z -- sqrt[z].prin ) ( Kahan's algorithm without overflow/underflow avoidance and without treating infinity. But it should handle signed zero properly. ) zdup [ rfnoname cell+ ] LITERAL (rf: y) rf! [ rfnoname ] LITERAL (rf: x) rf! |z| [ rfnoname ] LITERAL rf@ |rf| rf+ rf2/ rf-sqrt (rf: rho=sqrt[[|z|+|x|]/2]) rfdup rf0= [ rfnoname cell+ ] LITERAL rf@ (rf: rho y) 0= IF ( rho<>0) rfover rf/ rf2/ (rf: rho eta=[y/rho]/2) [ rfnoname ] LITERAL rf@ rf0< IF ( x<0) |rf| rfswap (rf: |eta| rho) [ rfnoname cell+ ] LITERAL rf@ rf-signbit IF ( y<0) -rf THEN (rf: |eta| |rho|*sgn[y]) THEN THEN ; : zsqrt (rf: z -- sqrt[z] | a b ) PRINCIPAL-ARG IF (rf: z -- sqrt[z] ) psqrt \ imag cut ELSE (rf: x y -- a b ) \ (a+ib)^2 = x+iy, real cut zdup (rf: -- z z) |z|^2 (rf: -- z |z|^2) rfdup rf0= IF rfdrop EXIT THEN (rf: -- z=0) rf-sqrt rfrot rfrot rf0< (rf: -- |z| x) ( -- sgn[y]) rftuck (rf: -- x |z| x) rf- rf2/ (rf: -- x [|z|-x]/2) rftuck rf+ (rf: -- [|z|-x]/2 [|z|+x]/2) rf-sqrt IF -rf THEN (rf: -- [|z|-x]/2 a) rfswap rf-sqrt (rf: -- a b) THEN ; \ Trigonometric and hyperbolic functions \ All stack patterns are (rf: z -- func[z] ). : zsinh ( This version is reasonably accurate and preserves signed zero. ) rf-sincos [ rfnoname ] LITERAL (rf: cos[y]) rf! [ rfnoname cell+ ] LITERAL (rf: sin[y]) rf! rfdup rf-sinh [ rfnoname ] LITERAL rf@ rf* (rf: x sh[x]cos[y]) rfswap rf-cosh [ rfnoname cell+ ] LITERAL rf@ rf* (rf: sh[x]cos[y] ch[x]sin[y]) ; : zcosh rf-sincos [ rfnoname ] LITERAL (rf: cos[y]) rf! [ rfnoname cell+ ] LITERAL (rf: sin[y]) rf! rfdup rf-cosh [ rfnoname ] LITERAL rf@ rf* (rf: x ch[x]cos[y]) rfswap rf-sinh [ rfnoname cell+ ] LITERAL rf@ rf* (rf: ch[x]cos[y] sh[x]sin[y]) ; : ztanh ( [1 + tan^2[y]] cosh[x] sinh[x] + i tan[y] tanh[z] = ----------------------------------------- 1 + [1 + tan^2[y]] sinh^2[x] ) rf-tan rfdup rf^2 1 u>rf rf+ (rf: x t=tan[y] b=1+t^2) [ rfnoname ] LITERAL rf! rfswap rfdup rf-sinh (rf: t x sh[x]) [ rfnoname cell+ ] LITERAL rf! rf-cosh [ rfnoname cell+ ] LITERAL rf@ (rf: t ch[x] sh[x]) rf* [ rfnoname ] LITERAL rf@ rf* rfswap (rf: c=ch[x]sh[x]b t) [ rfnoname ] LITERAL rf@ [ rfnoname cell+ ] LITERAL rf@ rf^2 rf* 1 u>rf rf+ (rf: c t 1+b*sh^2[x]) z/f ; : zcoth ztanh 1/z ; : zsin i* zsinh -i* ; : zcos i* zcosh ; : ztan i* ztanh -i* ; : zcot i* zcoth i* ; \ *** 9 INVERSE FUNCTIONS --- In the following, we use phrases equivalent to rf" 1.0" x+ instead of rf" 1.0" rf" 0.0" 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 (rf: z -- ln[z+sqrt[[z+i][z-i]] ) ( This is more stable for signed zero than the version with z^2+1. ) zdup 1 u>rf rf+ zover 1 u>rf rf- z* psqrt z+ pln ; : zacosh (rf: z -- 2ln[sqrt[[z+1]/2]+sqrt[[z-1]/2] ) zdup 1 u>rf x- z2/ psqrt zswap 1 u>rf x+ z2/ psqrt z+ pln z2* ; : zatanh (rf: z -- [ln[1+z]-ln[1-z]]/2 ) zdup 1 u>rf x+ pln zswap 1 u>rf x- znegate pln z- z2/ ; : zacoth (rf: z -- [ln[-1-z]-ln[1-z]]/2 ) znegate zdup 1 u>rf x- pln zswap 1 u>rf x+ pln z- z2/ ; : zasin (rf: z -- -iln[iz+sqrt[1-z^2]] ) i* zasinh -i* ; : zacos (rf: z -- pi/2-asin[z] ) zasin znegate pi/2 x+ ; : zatan (rf: z -- [ln[1+iz]-ln[1-iz]]/2i ) i* zatanh -i* ; : zacot (rf: z -- [ln[[z+i]/[z-i]]/2i ) -i* zacoth -i* ; \ *** 10 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. http://www.convertit.com/Go/ConvertIt/Reference/AMS55.ASP 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. http://www.cims.nyu.edu/~dbindel/class/cs279/slits02.pdf http://www.cims.nyu.edu/~dbindel/class/cs279/slits34.pdf http://www.cims.nyu.edu/~dbindel/class/cs279/slits56.pdf 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. http://portal.acm.org/citation.cfm?doid=362001.362023 http://www.apmaths.uwo.ca/~djeffrey/Offprints/couth.pdf ---