(         Title:  ANS Forth Complex Arithmetic Lexicon
                  with Kahan Algorithms
           File:  complex-kahan.fs
        Version:  0.8.9
Original Author:  Julian V. Noble
    Modified by:  David N. Williams
        License:  LGPL
  Last revision:  April 24, 2005

The original code is
)
\ Copyright (C) 1998  Julian V. Noble
(
Modifications not derived from the original, including
expressions of the Kahan algorithms, are
)
\ Copyright (C) 2002-2005  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 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.

Please see the file POLITENESS included with this distribution.

ANS Forth environmental dependences:

  1. Requires the Floating-Point and Floating-Point extension
     word sets.
  2. Assumes a separate floating point stack.
  3. Does not construct a separate complex number stack.
  4. Case insensitivity for standard words.
  5. The ANS FORTH default floats must be IEEE 754/854 to work
     correctly with signed zero/infinity.
  6. Some of the Kahan algorithms are optimized for radix 2.
  7. Machine parameters like EPSILON and a few others have to be
     uncommented, or supplied if they are not 32-bit or 64-bit
     IEEE 754.  The parameter section assumes that raw floats
     are represented in memory as contiguous sets of bytes that
     are totally big-endian or little-endian, with the raw
     representation of -1E nonzero in the logical first byte and
     zero in the logical last byte.
  8. If not present in the host FORTH, SIGNBIT assumes default
     floats to be IEEE 754 32 or 64 bits, and also assumes the
     cell size to be a multiple 32 bits.

Complex numbers z=x+iy are stored on the fp stack as [f: x y],
also written as [f: z].

Angles are in radians.

This code uses higher-accuracy algorithms by William Kahan [1].
It uses the principal argument, with -pi < arg <= pi, for ARG,
ZSQRT, ZLN, and Z^.  The Kahan alorithms implement an OpenMath
compliant treatment of principal expressions, branch cuts, and
branches for the elementary functions [2], which is based on
Abramowitz and Stegun [3].

Currently the only high-level way we have to save and restore
IEEE underflow and overflow flags, and to use IEEE scaling
functions to avoid underflow/overflow, is through a pfe
environment, FPIEEE-EXT.  The parts of the Kahan algorithms that
do that are conditionally compiled here, when the environment is
present.  Kahan says that the algorithms degrade gracefully
under simple omission of these procedures, which is what happens
here when FPIEEE-EXT is absent.

Kahan pays attention to signed zero, where available in IEEE
754/854 implementations.  We address that in this file as follows.

1. ZSINH, ZASINH, ZTANH, ZATANH, ZSIN, ZASIN, ZTAN, and ZATAN
   conserve the sign of zero, and 1/Z and the minimal
   computation words Z*F, Z/F, F*Z, F/Z, Z*I*F, -I*Z/F, I*F*Z,
   I*F/Z do the right thing.

2. We would like the analytic functions that are real and
   analytic on the real axis to do the right thing for the sign
   of the zero imaginary part.  This works for all of our tests
   so far.

3. For functions having branch cuts, 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
s" FPIEEE-EXT" environment? [IF] ( flag) drop
  true [ELSE] false [THEN] constant FPIEEE immediate

-1E pad faligned dup f! c@
[IF]   true  \ memory floats are big endian
[ELSE] false \ memory floats are little endian
[THEN] constant BIG-ENDIAN immediate

\ Wil Baden's cheap, nestable locals written for floats:
: framef  ( -- addr )  here falign 8 floats allot ;
: unframe ( addr -- )  here - allot ; 

: fpframe{  ( -- )  postpone framef postpone >r ; immediate
: }frame    ( -- )  postpone r> postpone unframe ; immediate

: &fa    here [ 8 floats ] literal - ;
: &fb    here [ 7 floats ] literal - ;
: fa     here [ 8 floats ] literal - f@ ;
: fb     here [ 7 floats ] literal - f@ ;
: fc     here [ 6 floats ] literal - f@ ;
: fd     here [ 5 floats ] literal - f@ ;
: fe     here [ 4 floats ] literal - f@ ;
: ff     here [ 3 floats ] literal - f@ ;
: fg     here [ 2 floats ] literal - f@ ;
: fh     here [ 1 floats ] literal - f@ ;
: to-fa  here [ 8 floats ] literal - f! ;
: to-fb  here [ 7 floats ] literal - f! ;
: to-fc  here [ 6 floats ] literal - f! ;
: to-fd  here [ 5 floats ] literal - f! ;
: to-fe  here [ 4 floats ] literal - f! ;
: to-ff  here [ 3 floats ] literal - f! ;
: to-fg  here [ 2 floats ] literal - f! ;
: to-fh  here [ 1 floats ] literal - f! ;

s" [UN]" pad c! pad char+ pad c@ move 
pad find nip 0=
[IF]
: [UN]  ( "name" -- flag )
(
Leave true if name is in the search order, else leave false.
)
  bl word find nip 0= ; immediate
[THEN]

\ non-Standard fp words jvn has found useful
[UN] s>f    [IF]  : s>f    s>d  d>f  ;           [THEN]
[UN] f-rot  [IF]  : f-rot  frot  frot  ;         [THEN]
[UN] fnip   [IF]  : fnip   fswap  fdrop  ;       [THEN]
[UN] ftuck  [IF]  : ftuck  fswap  fover ;        [THEN]
[UN] 1/f    [IF]  : 1/f    1E  fswap  f/ ;       [THEN]
[UN] f^2    [IF]  : f^2    fdup  f*  ;           [THEN]
\ added by dnw
[UN] (f:    [IF]  : (f:   postpone ( ; IMMEDIATE [THEN]
[UN] f2*    [IF]
     FPIEEE [IF]  : f2*   1 ldexp ;
          [ELSE]  : f2*   fdup f+ ;       [THEN] [THEN]
[UN] f2/    [IF]
     FPIEEE [IF]  : f2/   -1 ldexp ;
          [ELSE]  : f2/   0.5E f* ;       [THEN] [THEN]
[UN] f0>    [IF]  : f0>   fnegate f0< ;          [THEN]
[UN] f=     [IF]  : f=    f- f0= ;               [THEN]    
[UN] f>     [IF]  : f>    fswap f< ;             [THEN]    

[UN] f^n [IF]
  : f^n    ( u f: x -- x^u )
    1E ( u) 0 ?DO fover f* LOOP fnip ;
[THEN]

[UN] signbit [IF]  \ emulate IEEE signbit(x)
  : signbit  (f: x -- s: minus? )
    fpframe{
      (f: x) to-fa
      BIG-ENDIAN 
      [IF]   &fa
      [ELSE] &fb 1- [THEN]
      c@ 128 and
    }frame ;
[THEN]

\ accuracies good up to the 21 significant digits required
\ for Intel extended double (Sparc requires 36)
 1.570796326794896619231E  fconstant pi/2
 6.283185307179586476925E  fconstant 2pi
 3.141592653589793238463E  fconstant pi
-3.141592653589793238463E  fconstant -pi
 1.414213562373095048802E  fconstant sqrt2
\ 0.7071067811865475244008E fconstant 1/sqrt2
 0.70710678118654752440084436210484903928E fconstant 1/sqrt2
\ 0.6931471805599453094172E fconstant ln2
 0.69314718055994530941723212145817656808E fconstant ln2

MARKER -EXACT-PARAMS

: hex-byte  ( caddr u -- caddr+2 u-2 byte )
  base @ >r hex
  over 0. rot 2 >number                 ( caddr u dbyte=0 addr #rem)
  ABORT" ***Not a hex string."
  ( dhi=0 addr) 2drop ( byte) >r 2 /string r> ( byte)
  r> base ! ;

: raw-hex>float  ( caddr u -- f: f )
  ( u) dup 2 floats <> ABORT" ***Wrong size for raw float string."
  pad faligned dup
BIG-ENDIAN [IF]
  float+ swap ( &end+1 &start)
  DO hex-byte i c! LOOP
[ELSE]
  float+ 1-  ( &start &end)
  DO hex-byte i c! -1 +LOOP
[THEN]
  ( caddr' 0) 2drop   pad faligned f@ ;

( If your floating-point system is not one of the two below,
you need to supply your own parameters.
)

true
[IF]    \ 64-bit IEEE 754 double precision

\ machine parameters
s" 3CB0000000000000" raw-hex>float (f: eps)
\ s" 3CA0000000000000" raw-hex>float (f: epsneg)
s" 7FEFFFFFFFFFFFFF" raw-hex>float (f: xmax)
s" 7FF0000000000000" raw-hex>float (f: +Inf)
s" FFF0000000000000" raw-hex>float (f: -Inf)
s" 7FF8000000000000" raw-hex>float (f: NaN)

\ used by Kahan algorithms
s" 4003504F333F9DE6" raw-hex>float (f: r2p1=1+sqrt[2])
s" 3CA21165F626CDD6" raw-hex>float (f: t2p1=[1+sqrt[2]]-r2p1)

[ELSE]  \ 32-bit IEEE 754 single precision

\ machine parameters
s" 34000000" raw-hex>float (f: eps)
s" 7F7FFFFF" raw-hex>float (f: xmax)
s" 7F800000" raw-hex>float (f: +Inf)
s" FF800000" raw-hex>float (f: -Inf)
s" 7FC00000" raw-hex>float (f: NaN)

\ used by Kahan algorithms
s" 401A8279" raw-hex>float (f: r2p1=1+sqrt[2])
s" 3419FCEF" raw-hex>float (f: t2p1=[1+sqrt[2]]-r2p1)

[THEN]

-EXACT-PARAMS

(f: eps xmax +Inf -Inf NaN r2p1 t2p1)
fconstant t2p1   fconstant r2p1
fconstant NaN    fconstant -Inf   fconstant Inf
fconstant xmax   fconstant epsilon

\ used by Kahan algorithms
2E epsilon f/     fconstant 2/epsilon
xmax fsqrt 4E f/  fconstant sqrt(xmax)/4
sqrt(xmax)/4 1/f  fconstant 4/sqrt(xmax)
xmax fasinh 4E f/ fconstant sinh(xmax)/4
1/sqrt2           fconstant T0  \ for ZLN
5E 4E f/          fconstant T1
3E                fconstant T2

FPIEEE [IF]
  4E 1E epsilon f- f* xmax f/ epsilon f/
      fconstant lambda/epsilon 
[THEN]

( In ANS Forth ATAN2 is ambiguous.  The following tries to
redefine it to use the principal branch.  It will be wrong if
ATAN2 uses a third convention, which we hope is unlikely.
)
  -1E 1E (f: y x) fatan2 f0>
  [IF] : fatan2  (f: y x -- atan2[y,x])
         ( x) fdup f0< FATAN2
         IF 2pi F- THEN ; [THEN]

( The x y order below follows the readability principle for
implementing functions in Forth based on other standards.
)
[UN] copysign [IF] \ emulate IEEE copysign(x,y)
  : copysign  ( f: x y -- |x|*sgn[y] )
    signbit fabs IF fnegate THEN ;
[THEN]

( Hidden, *nonnestable* temporary storage for stuff from the
fpstack.
)
  falign here 3 floats allot value noname

\ *** LOAD, STORE

: z@  ( zaddr --  f: -- z )  dup  f@  float+  f@ ;
: z!  ( zaddr --  f: z -- )  dup  float+  f!  f! ;
: x@  ( zaddr -- f: x )  f@ ;
: y@  ( zaddr -- f: y )  float+ f@ ;
: x!  ( zaddr f: x -- )  f! ;
: y!  ( zaddr f: y -- )  float+ f! ; 

\ cheap locals for complex values, sharing the floats space:
: za     here [ 8 floats ] literal - z@ ;
: zb     here [ 6 floats ] literal - z@ ;
: zc     here [ 4 floats ] literal - z@ ;
: zd     here [ 2 floats ] literal - z@ ;
: to-za  here [ 8 floats ] literal - z! ;
: to-zb  here [ 6 floats ] literal - z! ;
: to-zc  here [ 4 floats ] literal - z! ;
: to-zd  here [ 2 floats ] literal - z! ;

\ *** COMPILE NAMED DATA

: zvariable  ( "name" -- faddr )
  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.   (f: x y -- )    \ emit complex #
  fswap f.
  (f: y) 1E fover copysign f0<  \ also handles -0E
  dup invert [char] + and   swap [char] - and   +
  emit ."  i " fabs f. ;

: zs.  (f: x y -- )    \ emit complex #, scientific notation
  fswap fs.
  (f: y) 1E fover copysign f0<  \ also handles -0E
  dup invert [char] + and   swap [char] - and   +
  emit ."  i " fabs fs. ;

: real  (f: x y -- y ) fdrop ;
: imag  (f: x y -- x ) fnip ;
: conjg (f: x y -- x -y )  fnegate  ;
: zdrop (f: z --)      fdrop fdrop ;
: zdup  (f: z -- z z ) fover fover ;

: zswap  (f: z1 z2 -- z2 z1)
  [ noname ] literal  f!  f-rot
  [ noname ] literal  f@  f-rot  ;

: znip  (f: z1 z2 -- z2 )  zswap  zdrop ;

: zover     (f: z1 z2 -- z1 z2 z1 )
  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)

: ztuck (f: z1 z2 -- z2 z1 z2 ) zswap zover ;
: z=    (f: z1 z2 -- s: flag )  frot f= f= and ;
: z0=   (f: z -- s: flag )      f0= f0= and ;

\ *** MINIMAL (MIXED) OPERATIONS
(
These words implement Kahan's rules [pp. 194, 195] for minimal
computation, to avoid gratuitous overflow, underflow, and signed
zero distortion in real or imaginary parts that don't need to be
computed.
)

: x+  (f: x y a -- x+a y )  frot f+ fswap ;
: y+  (f: x y a -- x y+a )  f+ ;
: x-  (f: x y a -- x-a y )  fnegate frot f+ fswap ;
: y-  (f: x y a -- x y-a )  f- ;

: z*>real  (f: z1 z2 -- re[z1*z2] )  \ compute z* without doing Im
  frot [ noname ] literal f* f!         (f: x1 x2)
  f* [ noname ] literal f@ f- ;         (f: x1*x2-y1*y2)

: z*>imag  (f: z1 z2 -- im[z1*z2] )  \ compute z* without doing Re
  [ noname ] literal f! f* fswap        (f: y1*x2 x1)
  [ noname ] literal f@ f* f+ ;         (f: y1*x2+x1*y2)

: z*f  (f: z f -- f*z )
  frot fover f*                         (f: y f x*f)
  f-rot f* ;

: z/f  (f: z f -- z/f )
  frot fover f/                         (f: y f x/f)
  f-rot f/ ;

: f*z  (f: f z -- f*z )
  f-rot fover f*                        (f: y f x*f)
  f-rot f* ;

: (f/z)  (f: f a b -- f/[b+[a/b]a] [a/b]f/[b+[a/b]a] )
(
A factor in F/Z, used with |a| <= |b|.
)
fpframe{
  zdup f/ to-fa                         (f: f a b) \ fa=a/b
  fswap fa f* f+ f/                     (f: u=f/[b+[a/b]a])
  fdup  fa f*                           (f: u v=u[a/b])
}frame ;

: f/z  (f: f z  -- f/z )
(
Kahan algorithm *without* due attention to spurious
over/underflows and zeros and infinities.
)
  zdup fabs fswap fabs                  (f: f z |y| |x|)
  f< IF fswap (f/z)
  ELSE  (f/z) fswap
  THEN  fnegate ;                       (f: f*x/|z|^2 -f*y/|z|^2)

: z*i*f  (f: z f -- z*if )
  frot fover f*                         (f: y f x*f)
  f-rot f* fnegate fswap ;              (f: -y*f x*f)

: -i*z/f  (f: z f -- z/[if] )
  frot fover f/ fnegate                 (f: y f -x/f)
  f-rot f/ fswap ;                      (f: y/f -x/f)

: i*f*z  (f: f z -- if*z )
  f-rot fover f*                        (f: y f x*f)
  f-rot f* fnegate fswap ;              (f: -y*f x*f)                           

: i*f/z  (f: f z -- if/z )
(
Kahan algorithm *without* due attention to spurious
over/underflows and zeros and infinities.
)
  zdup fabs fswap fabs                  (f: f z |y| |x|)
  f< IF fswap (f/z)
  ELSE  (f/z) fswap
  THEN fswap ;                          (f: f*y/|z|^2 f*x/|z|^2)

\ *** COMPLEX ARITHMETIC

: z*    (f: x y u v -- x*u-y*v  x*v+y*u)
(
Uses the algorithm followed by jvn:
  [x+iy]*[u+iv] = [[x+y]*u - y*[u+v]] + i[[x+y]*u + x*[v-u]]
Requires 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  (f: x y -- -x -y )  fswap fnegate fswap fnegate ;

: z-  (f: a b x y -- a-x b-y )
(
Kahan says to conserve signed zero, write -y+b for b-y instead of
-[y-b].
)
  fnegate frot f+ f-rot f- fswap ;

: 1/z  (f: z -- 1/z )        1E f-rot f/z ;
: z/   (f: z1 z2 -- z1/z2 )  1/z  z* ;
: z2/  (f: z -- z/2 )  f2/   fswap  f2/  fswap  ;
: z2*  (f: z -- z*2 )  f2*   fswap  f2*  fswap  ;
: i*   (f: x+iy -- -y+ix)    fnegate fswap ;
: -i*  (f: x+iy -- y-ix)     fswap fnegate ;

: z^2  (f: z -- z^2 )
(
Kahan algorithm without removal of any spurious NaN created by
overflow.  It deliberately uses [x-y][x+y] instead of x^2-y^2
for the real part.
)
  zdup f+ [ noname ] literal f!
  zdup f- [ noname ] literal f@ f*      (f: x y [x-y][x+y])
  f-rot f* f2* ;

: |z|^2  (f: z -- x^2+y^2 )  f^2 fswap f^2 f+  ;

0 [IF]
: z^3   zdup  z^2  z* ;
: z^4   z^2  z^2  ;
[THEN]

: 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 ;

\ *** COMPLEX ELEMENTARY FUNCTIONS

: |z|   (f: x y -- |z|)
(
This is Kahan's algorithm, which treats infinity, the precision
of 1+sqrt[2], and the roundoff threshhold for large ratios of
the maximum and minimum absolute values of x and y.  When the
FPIEEE word list is present, it also treats the IEEE underflow
and invalid flags.
)
fpframe{
FPIEEE [IF]
  FE-INVALID fegetexceptflag            ( invalid?)
[THEN]
  fabs fswap fabs
  zdup f<
  IF  fswap
  THEN  \ a >= b >= 0 if not NaN        (f: a b )
  fdup Inf 0E f~
  IF  fnip Inf                          (f: a=inf b=inf)
  THEN
   
  zdup to-fb  to-fa
  f- fdup fa                            (f: t=a-b t a)
  fdup Inf 0E f~   0E f~ or             (f: t) 

  ( t=a.or.a=inf?)
  IF    \ b is negligible
    (f: t) fdrop 0E                     (f: s=0) 

  ELSE  \ b is not negligible
FPIEEE [IF]
    FE-UNDERFLOW fegetexceptflag        ( invalid? underflow?)
[THEN]
    fa fover                            (f: t b t)

    f< IF  \ 2 < a/b
      (f: t) fdrop
      fa fb f/                          (f: u=a/b)
      fdup 2/epsilon f<
      IF  \ a/b < 2/epsilon
        fdup f^2 1E f+ fsqrt f+         (f: s=u+sqrt[1+u^2])
      ELSE
        fdrop 0E                        (f: s=0)
      THEN

    ELSE  \ 1 <= a/b <= 2
      (f: t)
      fb f/                             (f: v=t/b)
      fdup fdup 2E f+ f*                (f: v w=[2+v]v)
      fdup 2E f+ fsqrt sqrt2 f+ f/ f+   (f: h=w/[sqrt[2+w]+sqrt2]+v)
      t2p1 f+ r2p1 f+                   (f: s=[t2p1+h]+r2p1)
      \ substitute v = u - 1 and w = u^2 - 1 to prove that
      \ this s expression is the same as the one above

    THEN
    (f: s)
    fb fswap f/                         (f: s=b/s)
FPIEEE [IF] ( invalid? underflow?)
    FE-UNDERFLOW fesetexceptflag  ( invalid?)
    \ any underflow is now deserved
[THEN]

  THEN
FPIEEE [IF]
  ( invalid?) FE-INVALID fesetexceptflag
[THEN]
  (f: s) fa  f+                         (f: s+a)
}frame ;

: zbox  (f: z -- box[z] )
(
Defined *only* for zero and infinite arguments.  This difffers
from Kahan's CBOX [1] by conserving signs when only one of x or
y is infinite, consistent with the other cases, and with its use
in his ARG.
)
  zdup z0=
  IF
    1E frot copysign fswap EXIT         (f: cps[1,x] y)
  THEN
  fover fabs Inf f=                     (f: z s: |x|=inf?)
  IF fdup fabs Inf f=                   (f: z s: |y|=inf?)
    IF
      1E frot copysign                  (f: y cps[1,x])
      1E frot copysign                  (f: cps[1,x] cps[1,y])
    ELSE
      fover 1E fswap copysign           (f: x y cps[1,x])
      f-rot fswap fabs f/               (f: cps[1,x] y/|x|)
    THEN
  ELSE fdup fabs Inf f=                 (f: z s: |y|=inf?)
    IF
      fswap fover fabs f/               (f: y x/|y|)
      1E frot copysign                  (f: x/|y| cps[1,y])
    ELSE
      zdrop 0E 0E 0E z/f                 (f: 0/0+i0/0)  \ invalid use
    THEN
  THEN ;

: arg  (f: x y -- arg[x+iy] )
(
Kahan's algorithm, including underflow signal suppression when
The FPIEEE word set is present.  He says better accuracy can be
obtained by further case reduction and identities like atan[y/x]
= pi/4 + atan[[y-x]/[y+x]].
)
fpframe{
  to-fb to-fa                           \ za=x+iy
  fa f0= fb f0= and
  IF
    1E fa copysign to-fa                \ fa=x=cps[1,x]
  THEN
  fa fabs fb fabs                       (f: |x| |y|)
  zdup Inf f= Inf f= or
  IF
    za zbox to-za                       \ za=z=zbox[z]
  THEN  \ leaves signs unchanged
  (f: |x| |y|) f<
  IF
    pi/2 fb copysign fa fb f/ fatan f-  (f: theta=cps[pi/2,y]-atan[x/y])
  ELSE
    fa f0<
    IF   \ x < 0
      pi fb copysign fb fa f/ fatan f+  (f: theta=cps[pi,y]+atan[y/x])
    ELSE
      fb fa f/ fatan                    (f: theta=atan[y/x])
    THEN
  THEN
FPIEEE [IF]
  fdup fabs 0.125E f> IF FE-UNDERFLOW feclearexcept THEN
[THEN]
}frame ;

FPIEEE [IF]
: zssqs  (f: z -- |z/2^k|^2 s: k )
(
Compute |[x+iy]/2^k|^2, scaled to avoid overflow or underflow,
and leave the scaling integer k.  Kahan's CSSQS, p. 200.
)
fpframe{
  (f: z) to-za                          \ fa=x, fb=y
  [ FE-OVERFLOW FE-UNDERFLOW or dup ]
      literal fegetexceptflag           ( flags)
      literal feclearexcept
  ( flags) >r
  fa f^2 fb f^2 f+ to-fc                \ fc=rho=x^2+y^2
  fc isnan fc isinf or
  fa isinf fb isinf or and
  IF Inf 0                              (f: rho=inf s: k=0)
  ELSE
    FE-OVERFLOW  fetestexcept
    FE-UNDERFLOW fetestexcept 0= invert
    fc (f: rho) lambda/epsilon f<
    and or
    IF
      fa fabs fb fabs fmax ilogb        ( k)
      dup negate dup
      ( -k) fa scalbn ( -k) fb scalbn
      |z|^2                             (f: rho.scaled s: k)
    ELSE fc 0                           (f: rho.unscaled s: 0)
    THEN
  THEN
  r> ( flags) [ FE-OVERFLOW FE-UNDERFLOW or ] literal
  fesetexceptflag
}frame ;
[THEN]

: >polar  (f: x+iy -- r theta )  zdup |z| f-rot  arg  ;
: polar>  (f: r theta -- x+iy )  fsincos frot z*f fswap  ;

FPIEEE [IF]
variable #k's 0 #k's !   \ DEBUG

: zlns  ( j f: z -- ln[z*2^j] )
(
Kahan, p. 202.  He says to use this for nonzero j only when |z|
is huge, and that it's good for functions like inverse trig and
hyperbolic which go like ln[2z] at huge |z|.  We blindly use his
choices of threshholds, T0 = 1/sqrt2, T1 = 5/4, T2 = 3, which he
says should depend on the accuracy of FLNP1.
)
fpframe{
  zdup zssqs to-fa                      (f: z s: j k)
                                        \ fa=rho=|z|^2-scaled
  ( k) dup IF 1 #k's +! THEN  \ DEBUG
  zdup arg to-fb                        \ fb=theta
  fabs fswap fabs                       (f: |y| |x|)
  zdup fmin to-fc                       \ fc=min[|x|,|y|]
  fmax to-fd                            \ fd=max[|x|,|y|]
  ( k) dup 0=
  fd T1 f<=   fa T2 f<   or and
  T0 fd f<    and
  IF
    ( j k) 2drop
    fd 1E f-   fd 1E f+ f*
        fc f^2 f+                       (f: a=[max-1][max+1]+min^2)
    flnp1 f2/                           (f: ln1p[a]/2)
  ELSE
    fa fln f2/                          (f: ln[rho]/2)
    ( j k) + s>d d>f ln2 f* f+
  THEN
  fb
}frame ;
[THEN]

: zln  (f: z -- ln|z|+iarg[z] )
(
Kahan with blind choices of threshholds T0 = 1/sqrt2, T1 = 5/4,
T2 = 3, which should depend on the accuracy of FLNP1, and with
large agument treatment when FPIEEE is true.
)
FPIEEE [IF]
  0 zlns ;
[ELSE]
fpframe{
  zdup |z|^2 to-fa                      \ fa=rho=|z|^2
  zdup arg to-fb                        \ fb=theta
  fabs fswap fabs                       (f: |y| |x|)
  zdup fmin to-fc                       \ fc=min[|x|,|y|]
  fmax to-fd                            \ fd=max[|x|,|y|]
  fd T1 f<=   fa T2 f<   or
  T0 fd f<    and
  IF
    fd 1E f-   fd 1E f+ f*
        fc f^2 f+                       (f: a=[max-1][max+1]+min^2)
    flnp1 f2/                           (f: ln1p[a]/2)
  ELSE
    fa fln f2/                          (f: ln[rho]/2)
  THEN
  fb
}frame ;
[THEN]

: zexp  (f: z -- exp[z] )  fsincos fswap frot fexp z*f ;

: z^  (f:  x y u v --  [x+iy]^[u+iv] )  zswap zln  z* zexp  ;

: zsqrt  (f: z -- sqrt[z] )
(
Kahan's algorithm, with overflow/underflow avoidance when FPIEEE
is true.  In either case, it should handle signed zero and
infinity properly.
)
fpframe{
  zdup to-za
FPIEEE [IF]
  zssqs                                 (f: |z|^2.scaled s: k)
  fa fa f=
  IF  \ x is not a NaN
    fsqrt fa fabs ( k) dup negate scalbn f+
  THEN
  ( k) dup 2 mod
  IF   \ k is odd
    1- 2/ ( [k-1]/2)
  ELSE \ k is even
    2/ 1- ( [k/2]-1)   f2*
  THEN
  fsqrt scalbn                          (f: rho=sqrt[[|z|+|x|]/2])
[ELSE]
  |z| fa fabs f+ f2/ fsqrt              (f: rho=sqrt[[|z|+|x|]/2])
[THEN]
  fdup f0= fb                           (f: rho y )
  0= IF  \ rho <> 0
    fdup fabs Inf f<>                   (f: rho eta=y)
    IF fover f/ f2/ THEN                (f: rho eta=[y/rho]/2)
    fa f0< IF  \ x < 0
      fabs fswap fb copysign            (f: |eta| cps[rho,y])
    THEN      
  THEN
}frame ;

\ *** TRIGONOMETRIC AND HYPERBOLIC FUNCTIONS

: zcosh  (f: z -- cosh[z] )
(
This version is reasonably accurate and preserves signed zero on
the real axis.
)
  fsincos [ noname ] literal (f: cos[y]) f!
          [ noname float+ ] literal (f: sin[y]) f!
  fdup fcosh [ noname ] literal f@ f*   (f: x cosh[x]cos[y])
  fswap fsinh
      [ noname float+ ] literal f@ f*   (f: cosh[x]cos[y] sinh[x]sin[y])
;

: zcos  (f: z -- cos[z] )  i* zcosh  ;

: zsinh  (f: z -- sinh[z] )
(
This version is reasonably accurate and preserves signed zero.
)
  fsincos [ noname ] literal (f: cos[y]) f!
          [ noname float+ ] literal (f: sin[y]) f!
  fdup fsinh [ noname ] literal f@ f*   (f: x sinh[x]cos[y])
  fswap fcosh
      [ noname float+ ] literal f@ f*   (f: sinh[x]cos[y] cosh[x]sin[y])
;

: zsin  (f: z -- sin[z] )  i* zsinh -i* ;

: ztanh  (f: z -- tanh[z] )
fpframe{
  fover fabs sinh(xmax)/4 f>
  IF  \ avoid overflow
    0E fswap copysign                   (f: x cps[0,y])
    fswap 1E fswap copysign fswap       (f: cps[1,x] cps[0,y])
  ELSE
    \ divide by zero signal should have been suppressed here
    ftan fdup to-fa                     (f: x t=tan[y]) \ fa=t
    f^2 1E f+ to-fb                     (f: x)   \ fb=beta=1+t^2=1/cos^2[y]
    fsinh fdup to-fc                    (f: s=sinh[x]) \ fc=s
    f^2 1E f+ fsqrt                     (f: rho=sqrt[1+s^2]=cosh[x])
    fa fabs Inf f=
    IF  \ |t|=Inf, signal when s=0 ok
      (f: rho) fc f/ fa 1/f             (f: rho/s+i/t)
    ELSE
      (f: rho) fb f* fc f* fa           (f: z=beta*rho*s+it)
      fb fc f^2 f* 1E f+ z/f            (f: z=z/[1+beta*s^2])  
    THEN
  THEN
}frame ;

: zcoth  (f: z -- 1/tanh[z] ) ztanh  1/z ;
: ztan   (f: z -- tan[z] )    i* ztanh -i* ;
: zcot   (f: z -- cot[z] )   -i* zcoth -i* ;

\ *** INVERSE TRIGONOMETRIC AND HYPERBOLIC FUNCTIONS

: zacos  (f: z -- u+iv=-2i*ln[sqrt[[1+z]/2]+i*sqrt[[1-z]/2]] )
(
The principal expression is computed as:

  zacos[z] = 2 * atan[ Re[sqrt[1-z]] / Re[sqrt[1+z]] ]
               + i * asinh[ Im[ [sqrt[1+z]]* * sqrt[1-z] ] ]
)
fpframe{
  fover 1E f+ fover zsqrt               (f: z c+id=sqrt[z+1])
  to-za                                 \ za=c+id
  znegate fswap 1E f+ fswap zsqrt       (f: a+ib=sqrt[1-z])
  fover fa f/                           (f: a+ib a/c)
  fatan f2* f-rot                       (f: u=2*atan[a/c] a+ib)
  za conjg z*>imag                      (f: u h=im[[a+ib][c-id]])
  fasinh                                (f: u v=asinh[h]) 
}frame ; 

: zacosh  (f: z -- u+iv=2*ln[sqrt[[z+1]/2]+sqrt[[z-1]/2]] )
(
The principal expression is computed as:

  zacosh[z] = asinh[ Re[ [sqrt[z-1]]* * sqrt[z+1] ] ]
               + 2i * atan[ Im[ sqrt[z-1] / Re[ sqrt[z+1] ] ]
)
fpframe{
  fover 1E f+ fover zsqrt               (f: z c+id=sqrt[z+1])
  to-za                                 \ za=c+id
  fswap 1E f- fswap zsqrt               (f: a+ib=sqrt[z-1])
  fdup fa f/                            (f: a+ib b/c)
  fatan f2* f-rot                       (f: v=2*atan[b/c] a+ib)
  conjg za z*>real                      (f: v h=re[[a-ib][c+id]])
  fasinh fswap                          (f: u=asinh[h] v)
}frame ;

: zasin  (f: z -- u+iv=-iln[iz+sqrt[1-z^2]] )
(
The principal expression is computed as:

  zasin[z] = atan[ x / Re[ sqrt[1-z] * sqrt[1+z] ] ]
               + i * asinh[ Im[ [sqrt[1-z]]* * sqrt[1+z] ] ]
)
fpframe{
  fover (f: x) fdup to-fc               (f: z x)
  1E f+ fover zsqrt                     (f: z c+id=sqrt[z+1])
  to-za                                 \ za=c+id
  znegate fswap 1E f+ fswap zsqrt       (f: a+ib=sqrt[1-z])
  zdup za z*>real                       (f: a+ib g=re[[a+ib]*[c+id]])
  fc                                    (f: a+ib g x)
  fswap f/ fatan f-rot conjg            (f: u=atan[x/g] a-ib)
  za z*>imag                            (f: u h=im[[a-ib]*[c+id]) 
  fasinh                                (f: u v=asinh[h])
}frame ;

: zasinh  (f: z -- ln[z+sqrt[1+z^2] )  i* zasin -i* ;

: zatanh  (f: z -- u+iv=[ln[1+z]-ln[1-z]]/2 )
fpframe{
  \ handle unsigned zero
  fover 1E fswap copysign               (f: x y beta=cps[1,x])
  fdup to-fa                            \ fa=beta
  z*f conjg                             (f: z=beta*conjg[z])

  zdup fabs sqrt(xmax)/4 f>
      sqrt(xmax)/4 f> or
  IF   \ avoid overflow
    pi/2 fover copysign f-rot           (f: v=cps[pi/2,y] z)
    1/z real fswap                      (f: u=re[1/z] v)
  ELSE
    fdup fabs 4/sqrt(xmax) f+           (f: z a=|y|+rho)
    to-fb                               \ fb=a
    fover (f: x) 1E f=                  (f: z)
    IF   \ x=1
      fnip                              (f: y)
      fdup fdup f* 4E f+ fsqrt fsqrt    (f: y b=sqrt[sqrt[4+y^2]])
      fb fsqrt f/ fln fswap             (f: u=ln[b/sqrt[a]] y)
      fb f2/ fatan pi/2 f+              (f: u y d=pi/2+atan[a/2])
      fswap copysign f2/                (f: u v=cps[d,y]/2)
    ELSE \ normal case, use ln1p(u) accurately even if u is tiny
      fswap fdup fnegate 1E f+          (f: y x 1-x)
      fover 1E f+                       (f: y x 1-x 1+x)
      fover f* to-fc                    (f: y x 1-x) \ fc=g=[1-x^2]
      f^2 fb f^2                        (f: y x h=[1-x]^2 a^2)
      fdup to-fb                        \ fb=a^2
      f+ f/ 4E f* flnp1 4E f/           (f: y u=ln1p[4x/[h+a^2]]/4)
      fswap f2*   fc  fb                (f: u 2y g a^2)
      f- fswap arg f2/                  (f: u v=arg[g-a^2+i2y]/2)
    THEN
  THEN

  conjg fa z*f                         (f: u+iv=beta*[u-iv])
}frame ;

: zatan  (f: z -- [ln[1+iz]-ln[1-iz]]/2i )  i* zatanh -i* ;

\ *** REFERENCES
(
1. 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.

2. 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.

3. 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.
)

[ELSE] .( ***Separate floating point stack not available.) [THEN]
[ELSE] .( ***Floating point not available.)                [THEN]