(       Title:  Complex Word Set Signed Zero,
                Signed Infinity, and NaN Tests
         File:  complex-ieee-test-xf.fs
       Author:  David N. Williams
      Version:  1.0.3
      License:  LGPL
Last revision:  December 12, 2010
)
\ Copyright (C) 2003, 2005, 2006, 2008-2010 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.

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.

Unattributed changes are by David N. Williams.  The last
revision date above may reflect cosmetic changes not logged
here.

Version 1.0.3
12Dec10 * Improved BITS/FLOAT.
	* Tuned to work with iForth double extended.

Version 1.0.2
 1Dec10 * Started from complex-ieee-test 1.0.2 29Nov10.
 8Dec10 * Released.

This code tests our revisions of Julian V. Noble's complex
arithmetic lexicon in complex.fs [FSL Algorithm #60] and
complex-kahan.fs, and our C version of the latter, which
implements pfe's COMPLEX-EXT module.

It tests the treatment of signed zero, signed infinity, and NaN,
assuming IEEE 754 double or double extended floats.  The code is
written to run on single float systems as well, but has not been
tested there.

It is intended to test for formal correctness, not high
accuracy.  At least for doubles, the accuracy is actually rather
good; but to support conclusions about accuracy would require
far more exhaustive tests.

Most of the tests are devoted to signed zero, where there are
three issues:

1. The preservation of zero for analytic functions obeying f[0]
   = 0.  There is an expectation in the literature for real
   functions of a real variable that f[+0] = +0 and f[-0] = -0
   when f'[0] > 0 [Kahan, p. 167], which we apply to complex
   functions when that's easy and seems natural.

2. Having the right sign of zero for the zero imaginary part of
   a function that is real and analytic on the real axis.  This
   property is not always easy to implement, because the
   addition and subtraction of such functions does not preserve
   the property.  For example, f = 2*z - z evaluated at z =
   x-i0, has Im f = +0, while Im z = -0.  The wrong sign can
   cause instability when a real analytic function is used in
   the argument of a function with a branch cut on the real
   axis.

   The attitude here is that we are pleasantly surprised when a
   function passes the test.  Actually only complex.fs fails any
   of the tests in this section, and then it's only a few. 
   We're not sure it's worthwhile to be more aggressive with
   Kahan algorithms in complex.fs.
   
3. Values on branch cuts on the real or imaginary axis.  Kahan's
   treatment of the complex elementary and inverse functions
   aims to handle that.  The tests here compare to explicit
   formulas worked out for the OpenMath principal expressions on
   their branch cuts.

"Gauge functions" are functions that we test against.  They are
defined independently here, sometimes in terms of already tested
functions.

This version uses a new extension of the Hayes tests to complex
numbers, called ztester, which allows individual comparison
modes and tolerances for the real and imaginary parts.  This is
especially helpful for testing signed zero output.  In all cases
where those results fail exact tests, they pass absolute
difference tests with a tolerance near the smallest nonzero
normal number, at least for doubles.

This code is compatible with Forth 200x, with an environmental
dependence on lower case.
)

[UNDEFINED] \\ [IF]
  : \\  ( -- )  -1 parse 2drop BEGIN refill 0= UNTIL ;
[THEN]
: .#  ( n -- )  cr ." #" . cr ;

1 CONSTANT NOBLE
2 CONSTANT KAHAN
3 CONSTANT PFE

\ USER-CONFIG: select complex library implementation
\  NOBLE
   KAHAN
\  PFE
  CONSTANT COMPLEX-TEST-LIB

: NOBLE? ( -- flag ) COMPLEX-TEST-LIB NOBLE = ;
: KAHAN? ( -- flag ) COMPLEX-TEST-LIB KAHAN = ;
: PFE?   ( -- flag ) COMPLEX-TEST-LIB PFE = ;

\ Ensure that fp is enabled for pfe in ttester:
  s" FLOATING-EXT" environment? [IF] ( flag) drop [ELSE]
  cr .( FLOATING POINT NOT AVAILABLE) cr ABORT [THEN]

s" ttester-xf.fs" INCLUDED
true VERBOSE !

s" tester-display.fs" INCLUDED
\ before first displayed line
PFE-HOST IFORTH-HOST or [IF] ?.cr [THEN]

PFE-HOST [IF]
  ?." Host is pfe
\ 0 REDEFINED-MSG !
[THEN]
GFORTH-HOST [IF]
  ?." Host is gforth
\ 0 WARNINGS !
[THEN]
IFORTH-HOST [IF]
  ?." Host is iForth
\ 0 WARNING !
[THEN]

NOBLE? [IF]
  s" complex.fs" INCLUDED  \ defines PRINCIPAL-ARG  
  true to PRINCIPAL-ARG
[THEN]

KAHAN? [IF]
  s" complex-kahan.fs" INCLUDED
[THEN]

PFE? [IF] \ only for pfe
  s" COMPLEX-EXT" environment?
  [IF]
    ( version) drop
    s" IEEEFP-EXT" environment? [IF] ( version) drop [THEN]
  [ELSE]
    cr .( COMPLEX-EXT not available.) cr ABORT
  [THEN]
[THEN]

s" ztester.fs" INCLUDED
decimal

[UNDEFINED] BITS/FLOAT [IF]
:NONAME  ( -- +n )
  s" MAX-FLOAT" environment? 0= ABORT" ***can't determine MAX-FLOAT"
  pad 3 represent 0= ABORT" ***MAX-FLOAT invalid" \ shouldn't happen
  ( exp plus?) drop
  CASE
      39 OF 32 ENDOF  \ .340282347e+39
     309 OF 64 ENDOF  \ .17976931348623157e+309
    4933 OF 80 ENDOF  \ .118973149535723176505e+4933
    ABORT" ***BITS/FLOAT must be 32, 64, or 80"
  ENDCASE ; EXECUTE CONSTANT BITS/FLOAT
[THEN]

?." Floats are " BITS/FLOAT ?. ?." bits 

BITS/FLOAT 80 = [IF]
  1.3e-18  \ 20.27 significant digits
  1e-4931  \ 3.36210314311209350626e-4932 min normal
  GFORTH-HOST [IF] 22 
  [ELSE] IFORTH-HOST
    [IF] 18 [ELSE] 21 [THEN]
  [THEN] set-precision
[THEN]

BITS/FLOAT 64 = [IF]
  1e-15   \ 16.95 significant digits
  1e-307  \ 2.2250738585072014e-308 min normal 
  GFORTH-HOST [IF] 18 [ELSE] 17 [THEN] set-precision
[THEN]

BITS/FLOAT 32 = [IF]  \ untested
  1.5e-7  \ 8.225 significant digits
  1.0e-37 \ 1.17549435e-38 min normal
  GFORTH-HOST [IF] 10 [ELSE] 9 [THEN] set-precision
[THEN]
FCONSTANT abs-err-default   FCONSTANT rel-err-default 

: SET-ABS-ERRS  ( f: err - )
  fdup ZT-XABS-ERROR f!   ZT-YABS-ERROR f! ;
: SET-REL-ERRS  ( f: err - )
  fdup ZT-XREL-ERROR f!   ZT-YREL-ERROR f! ;
rel-err-default SET-REL-ERRS
abs-err-default SET-ABS-ERRS

: SET-ZT-MODE-XYABS ( -- ) SET-ZT-MODE-XABS SET-ZT-MODE-YABS ;

:NONAME  ( lib.case -- )
  blue-text
  s" IEEE-754 special value tests for " ?type
  COMPLEX-TEST-LIB CASE
    NOBLE OF s" complex.fs" ENDOF
    KAHAN OF s" complex-kahan.fs" ENDOF
    PFE   OF s" pfe COMPLEX-EXT" ENDOF
    normal-text ABORT" ***No complex words loaded"
  ENDCASE ?type normal-text ?.cr ; EXECUTE 

[DEFINED] IEEEFP [IF] IEEEFP [IF]
  ?." Using IEEEFP-EXT environment
[ELSE]
  ?." Not using IEEEFP-EXT environment
[THEN] [THEN]


\ ***ERROR DISPLAY

VARIABLE XT-#ERRORS    0 XT-#ERRORS !

: XT-ERROR1  ( c-addr u -- )
\ Display an error message followed by the line that had the
\ error.
  red-text 1 XT-#ERRORS +! XT-ERROR-DEFAULT normal-text ;
' XT-ERROR1 XT-ERROR-XT !

: xt-a.  ( &array c-addr u -- )
  blue-text type
  ( &array) dup >r @ dup
  IF   1- r> cell+ swap cells over +
       DO i @ . -1 cells +LOOP
  ELSE r> 2drop ." none" THEN cr normal-text ;

: XT-ERROR2  ( c-addr u -- )
  XT-ERROR1
  XT-RESULTS s"   xresults: " xt-a.
  XT-GOALS   s"     xgoals: " xt-a. ;

VARIABLE FT-#ERRORS    0 FT-#ERRORS !

: FT-ERROR1  ( c-addr u -- )
\ Display an error message followed by the line that had the
\ error.
  red-text 1 FT-#ERRORS +! FT-ERROR-DEFAULT normal-text ;
\ ' FT-ERROR1 FT-ERROR-XT !

\ fancy error reporting
: ft-a.  ( 'array c-addr u -- )
  blue-text type
  ( 'array) dup @ swap cell+ faligned swap ?dup
  IF   ( n) 1 swap  \ display deepest first 
       DO ( 'floats) i 1- floats over + f@
          i FT-ERROR-INDEX @ =
          IF red-text f. blue-text ELSE f. THEN
       -1 +LOOP
  ELSE ." none" THEN drop cr normal-text ;

: FT-ERROR2  ( c-addr u -- )
  FT-ERROR1
  FT-RESULTS s"   fresults: " ft-a.
  FT-GOALS   s"     fgoals: " ft-a. ;
' FT-ERROR2 FT-ERROR-XT !

VARIABLE ZT-#ERRORS    0 ZT-#ERRORS !
: ZT-ERROR1  ( c-addr u -- )
  red-text 1 ZT-#ERRORS +! ZT-ERROR-DEFAULT normal-text ;
\ ' ZT-ERROR1 ZT-ERROR-XT !

\ fancy error reporting
: zt-a.  ( 'a c-addr u -- )
  blue-text type
  ( 'a) dup @ swap cell+ faligned swap ?dup
  IF   ( n) 1 swap  \ display deepest first 
       DO ( 'z1) i 1- 2* floats over + z@
          i ZT-ERROR-INDEX @ =
          IF red-text zs. blue-text ELSE zs. THEN
       -1 +LOOP
  ELSE ." none" THEN drop cr normal-text ;

: ZT-ERROR2  ( c-addr u -- )
  ZT-ERROR1
  ZT-RESULTS s"   zresults: " zt-a.
  ZT-GOALS   s"     zgoals: " zt-a. ;
' ZT-ERROR2 ZT-ERROR-XT !

: ?.errors  ( -- )
  VERBOSE @
  IF blue-text ." XT-ERRORS: "   normal-text XT-#ERRORS @ .
     blue-text ."   FT-ERRORS: " normal-text FT-#ERRORS @ .
     blue-text ."   ZT-ERRORS: " normal-text ZT-#ERRORS @ .
  THEN ;

\ ***END ERROR DISPLAY


[UNDEFINED] d> [IF]
  : d>  2over 2over d= >r d< r> or 0= ; [THEN]

[UNDEFINED] fsgn0* [IF]
  : fsgn0*  ( f: f 0 -- f*sgn[0] )
    -0e 0e f~ IF fnegate THEN ;  [THEN]

[UNDEFINED] fsgn [IF]
  : fsgn  ( f: f -- sgn[f] )
    f0< IF -1e ELSE 1e THEN ; [THEN]

[UNDEFINED] fnan? [IF]  \ iForth and pfe IEEEFP-EXT define it
  -1e pad faligned f! pad faligned c@ 0=
  CONSTANT LITTLE-ENDIAN immediate

  : sf>raw32  ( f: sf -- s: raw32 )
    pad sfaligned sf!   pad sfaligned @ ;

  : df>raw64  ( f: df -- s: raw64 )
    pad dfaligned df!   pad dfaligned 2@
    LITTLE-ENDIAN [IF] swap [THEN] ;

  : fnan?  ( f: f -- s: nan? )  \ either quiet or signaling
    fabs
    [ BITS/FLOAT 64 = ] [IF]
      df>raw64 0 $7FF00000 d>
    [THEN]
    [ BITS/FLOAT 32 = ] [IF]    
      sf>raw32 $7F800000 >
    [THEN] ;
[THEN]

[UNDEFINED] +Inf [IF]  1e 0e f/ FCONSTANT +Inf [THEN]
[UNDEFINED] -Inf [IF] -1e 0e f/ FCONSTANT -Inf [THEN]

t{ 0e 0e f/ fnan? -> true }t
t{ 0e 0e f/ fnegate fnan? -> true }t
t{ +Inf fnan? -> false }t
t{ -Inf fnan? -> false }t
t{  3e  fnan? -> false }t
t{ -3e  fnan? -> false }t
[DEFINED] +NaN [IF]
t{ +NaN fnan? -> true }t
[THEN]
[DEFINED] -NaN [IF]
t{ -NaN fnan? -> true }t
[THEN]

\ 21 significant digits
 0.785398163397448309616E0 fconstant  pi/4
-0.785398163397448309616E0 fconstant -pi/4

[UNDEFINED] pi/2 [IF] 
 0.157079632679489661923E1 fconstant  pi/2
[THEN]
-0.157079632679489661923E1 fconstant -pi/2

 0.235619449019234492885E1 fconstant  3pi/4
-0.235619449019234492885E1 fconstant -3pi/4

[UNDEFINED] pi [IF]
 0.314159265358979323846E1 fconstant  pi
[THEN]
[UNDEFINED] -pi [IF]
-0.314159265358979323846E1 fconstant -pi
[THEN]

 0.141421356237309504880E1 fconstant  rt2
-0.141421356237309504880E1 fconstant -rt2
 0.707106781186547524401E0 fconstant  1/rt2
-0.707106781186547524401E0 fconstant -1/rt2

defer gauge
defer func

zvariable zatemp   zvariable zbtemp

: }GAUGE  ( f: z -- )
(
Compare the functions whose xt's are in FUNC and GAUGE.
) 
  fdepth 2 <> ABORT" EXACTLY ONE COMPLEX NUMBER NEEDED"
  zdup zatemp z! func z->
  zatemp z@ gauge }z ;

: }2GAUGE  ( f: z1 z2 -- )
(
Same as above with 2 complex arguments.  Note that the FUNC and
GAUGE results are still a single complex number.
) 
  fdepth 4 <> ABORT" EXACTLY TWO COMPLEX NUMBERS NEEDED"
  zbtemp z! zatemp z!
  zatemp z@ zbtemp z@ func z->
  zatemp z@ zbtemp z@ gauge }z ;

testing  F~  ZNEGATE  CONJG  I*  -I*

SET-FT-MODE-EXACT

t{ -0e  0e  0e f~ -> false }t
t{ -0e -0e  0e f~ -> true  }t

t{ 0e 0e znegate -> -0e -0e }t
t{ 0e 0e z2*     ->  0e  0e }t
t{ 0e 0e z2/     ->  0e  0e }t
t{ 0e 0e i*      -> -0e  0e }t
t{ 0e 0e -i*     ->  0e -0e }t

SET-ZT-MODE-ZEXACT

z{  0e  0e conjg z->  0e -0e }z
z{ -0e  0e conjg z-> -0e -0e }z
z{  0e -0e conjg z->  0e  0e }z
z{ -0e -0e conjg z-> -0e  0e }z

z{  0e  0e i* z-> -0e  0e }z
z{ -0e  0e i* z-> -0e -0e }z
z{  0e -0e i* z->  0e  0e }z
z{ -0e -0e i* z->  0e -0e }z

z{  0e  0e -i* z->  0e -0e }z
z{ -0e  0e -i* z->  0e  0e }z
z{  0e -0e -i* z-> -0e -0e }z
z{ -0e -0e -i* z-> -0e  0e }z

[DEFINED] i*f*z [IF]
testing  Z*F Z/F F*Z F/Z Z*I*F -I*Z/F I*F*Z I*F/Z Z2* Z2/

SET-ZT-MODE-ZEXACT

t{  2e  0e  0e frot z*f ->  0e  0e }t
t{  2e -0e  0e frot z*f -> -0e  0e }t
t{  3e  0e -0e frot z*f ->  0e -0e }t
t{  3e -0e -0e frot z*f -> -0e -0e }t
t{ -2e  0e  0e frot z*f -> -0e -0e }t
t{ -2e -0e  0e frot z*f ->  0e -0e }t
t{ -3e  0e -0e frot z*f -> -0e  0e }t
t{ -3e -0e -0e frot z*f ->  0e  0e }t

t{  2e  0e  0e frot z/f ->  0e  0e }t
t{  2e -0e  0e frot z/f -> -0e  0e }t
t{  3e  0e -0e frot z/f ->  0e -0e }t
t{  3e -0e -0e frot z/f -> -0e -0e }t
t{ -2e  0e  0e frot z/f -> -0e -0e }t
t{ -2e -0e  0e frot z/f ->  0e -0e }t
t{ -3e  0e -0e frot z/f -> -0e  0e }t
t{ -3e -0e -0e frot z/f ->  0e  0e }t

t{  2e  0e  0e f*z ->  0e  0e }t
t{  2e -0e  0e f*z -> -0e  0e }t
t{  3e  0e -0e f*z ->  0e -0e }t
t{  3e -0e -0e f*z -> -0e -0e }t
t{ -2e  0e  0e f*z -> -0e -0e }t
t{ -2e -0e  0e f*z ->  0e -0e }t
t{ -3e  0e -0e f*z -> -0e  0e }t
t{ -3e -0e -0e f*z ->  0e  0e }t

t{  2e  3e  0e f/z ->  2e  3e f/ -0e }t
t{  2e  3e -0e f/z ->  2e  3e f/  0e }t
t{  2e -3e  0e f/z ->  2e -3e f/ -0e }t
t{  2e -3e -0e f/z ->  2e -3e f/  0e }t
t{ -2e  3e  0e f/z -> -2e  3e f/  0e }t
t{ -2e  3e -0e f/z -> -2e  3e f/ -0e }t
t{ -2e -3e  0e f/z -> -2e -3e f/  0e }t
t{ -2e -3e -0e f/z -> -2e -3e f/ -0e }t

t{  2e  0e  0e frot z*i*f ->  0e  0e i* }t
t{  2e -0e  0e frot z*i*f -> -0e  0e i* }t
t{  3e  0e -0e frot z*i*f ->  0e -0e i* }t
t{  3e -0e -0e frot z*i*f -> -0e -0e i* }t
t{ -2e  0e  0e frot z*i*f -> -0e -0e i* }t
t{ -2e -0e  0e frot z*i*f ->  0e -0e i* }t
t{ -3e  0e -0e frot z*i*f -> -0e  0e i* }t
t{ -3e -0e -0e frot z*i*f ->  0e  0e i* }t

t{  2e  0e  0e frot -i*z/f ->  0e  0e -i* }t
t{  2e -0e  0e frot -i*z/f -> -0e  0e -i* }t
t{  3e  0e -0e frot -i*z/f ->  0e -0e -i* }t
t{  3e -0e -0e frot -i*z/f -> -0e -0e -i* }t
t{ -2e  0e  0e frot -i*z/f -> -0e -0e -i* }t
t{ -2e -0e  0e frot -i*z/f ->  0e -0e -i* }t
t{ -3e  0e -0e frot -i*z/f -> -0e  0e -i* }t
t{ -3e -0e -0e frot -i*z/f ->  0e  0e -i* }t

t{  2e  0e  0e i*f*z ->  0e  0e i* }t
t{  2e -0e  0e i*f*z -> -0e  0e i* }t
t{  3e  0e -0e i*f*z ->  0e -0e i* }t
t{  3e -0e -0e i*f*z -> -0e -0e i* }t
t{ -2e  0e  0e i*f*z -> -0e -0e i* }t
t{ -2e -0e  0e i*f*z ->  0e -0e i* }t
t{ -3e  0e -0e i*f*z -> -0e  0e i* }t
t{ -3e -0e -0e i*f*z ->  0e  0e i* }t

t{  2e  3e  0e i*f/z ->  2e  3e f/ -0e i* }t
t{  2e  3e -0e i*f/z ->  2e  3e f/  0e i* }t
t{  2e -3e  0e i*f/z ->  2e -3e f/ -0e i* }t
t{  2e -3e -0e i*f/z ->  2e -3e f/  0e i* }t
t{ -2e  3e  0e i*f/z -> -2e  3e f/  0e i* }t
t{ -2e  3e -0e i*f/z -> -2e  3e f/ -0e i* }t
t{ -2e -3e  0e i*f/z -> -2e -3e f/  0e i* }t
t{ -2e -3e -0e i*f/z -> -2e -3e f/ -0e i* }t

t{  0e  0e z2* ->  0e  0e }t
t{ -0e  0e z2* -> -0e  0e }t
t{  0e -0e z2* ->  0e -0e }t
t{ -0e -0e z2* -> -0e -0e }t

t{  0e  0e z2/ ->  0e  0e }t
t{ -0e  0e z2/ -> -0e  0e }t
t{  0e -0e z2/ ->  0e -0e }t
t{ -0e -0e z2/ -> -0e -0e }t
[THEN]

testing  ARG >POLAR POLAR> Z^N (n = 0, 1)

SET-ZT-MODE-ZEXACT
t{   1e   0e  arg ->  0e }t
t{   1e  -0e  arg -> -0e }t
t{  -1e   0e  arg ->  pi }t
t{  -1e  -0e  arg -> -pi }t
t{   0e   0e  arg ->  0e }t
t{  -0e   0e  arg ->  pi }t
t{   0e  -0e  arg -> -0e }t
t{  -0e  -0e  arg -> -pi }t
t{ +Inf   0e  arg ->  0e }t
t{ -Inf   0e  arg ->  pi }t
t{ +Inf  -0e  arg -> -0e }t
t{ -Inf  -0e  arg -> -pi }t
t{ +Inf   3e  arg ->  0e }t
t{ -Inf   3e  arg ->  pi }t
t{ +Inf  -5e  arg -> -0e }t
t{ -Inf  -5e  arg -> -pi }t
t{   0e +Inf  arg ->  pi/2 }t
t{  -0e +Inf  arg ->  pi/2 }t
t{   0e -Inf  arg -> -pi/2 }t
t{  -0e -Inf  arg -> -pi/2 }t
t{   3e +Inf  arg ->  pi/2 }t
t{  -3e +Inf  arg ->  pi/2 }t
t{   3e -Inf  arg -> -pi/2 }t
t{  -3e -Inf  arg -> -pi/2 }t
t{ +Inf +Inf  arg ->  pi/4 }t
t{ -Inf +Inf  arg ->  3pi/4 }t
t{ +Inf -Inf  arg -> -pi/4 }t
t{ -Inf -Inf  arg -> -3pi/4 }t

SET-ZT-MODE-ZEXACT
t{  0e  0e  >polar -> 0e  0e }t
t{  0e -0e  >polar -> 0e -0e }t
t{ -0e  0e  >polar -> 0e  pi }t
t{ -0e -0e  >polar -> 0e -pi }t
t{  1e  0e  >polar -> 1e  0e }t
t{  1e -0e  >polar -> 1e -0e }t
t{ -1e  0e  >polar -> 1e  pi }t
t{ -1e -0e  >polar -> 1e -pi }t

SET-ZT-MODE-ZEXACT
t{ 0e  0e   polar> ->  0e  0e }t
t{ 0e -0e   polar> ->  0e -0e }t
t{ 1e  0e   polar> ->  1e  0e }t
t{ 1e -0e   polar> ->  1e -0e }t

SET-ZT-MODE-ZEXACT
t{  0e  0e  0 z^n ->  1e  0e }t
t{  1e  0e  0 z^n ->  1e  0e }t
t{ -1e  0e  0 z^n ->  1e  0e }t
t{  0e  1e  0 z^n ->  1e  0e }t
t{  0e -1e  0 z^n ->  1e  0e }t

t{  0e  0e  1 z^n ->  0e  0e }t
t{  1e  0e  1 z^n ->  1e  0e }t
t{ -1e  0e  1 z^n -> -1e  0e }t
t{  0e  1e  1 z^n ->  0e  1e }t
t{  0e -1e  1 z^n ->  0e -1e }t

testing at f(0) = 0:  ZSINH ZSIN ZTANH ZTAN
?.(   ZASINH  ZASIN ZATANH ZATAN) ?.cr

\ functions with f(0) = 0 and f'(0) > 0

:noname  ; is gauge

' zsinh is func
SET-ZT-MODE-ZEXACT
z{  0e  0e  }gauge
z{  0e -0e  }gauge
z{ -0e  0e  }gauge
z{ -0e -0e  }gauge

' zsin is func
SET-ZT-MODE-ZEXACT
z{  0e  0e  }gauge
z{  0e -0e  }gauge
z{ -0e  0e  }gauge
z{ -0e -0e  }gauge

' ztanh is func
SET-ZT-MODE-ZEXACT
z{  0e  0e  }gauge
z{  0e -0e  }gauge
z{ -0e  0e  }gauge
z{ -0e -0e  }gauge

' ztan is func
SET-ZT-MODE-ZEXACT
z{  0e  0e  }gauge
z{  0e -0e  }gauge
z{ -0e  0e  }gauge
z{ -0e -0e  }gauge

' zasinh is func
SET-ZT-MODE-ZEXACT
z{  0e  0e  }gauge
NOBLE? [IF] SET-ZT-MODE-XYABS [THEN]
z{  0e -0e  }gauge
z{ -0e  0e  }gauge
z{ -0e -0e  }gauge

' zasin is func
SET-ZT-MODE-ZEXACT
NOBLE? [IF] SET-ZT-MODE-XYABS [THEN]
z{  0e  0e  }gauge
z{  0e -0e  }gauge
z{ -0e  0e  }gauge
z{ -0e -0e  }gauge

' zatanh is func
SET-ZT-MODE-ZEXACT
z{  0e  0e  }gauge
z{  0e -0e  }gauge
NOBLE? [IF] SET-ZT-MODE-XYABS [THEN]
z{ -0e  0e  }gauge
z{ -0e -0e  }gauge

' zatan is func
SET-ZT-MODE-ZEXACT
z{  0e -0e  }gauge
z{ -0e -0e  }gauge
NOBLE? [IF] SET-ZT-MODE-XYABS [THEN]
z{  0e  0e  }gauge
z{ -0e  0e  }gauge

testing real analyticity:  1/Z ZSQRT ZLN ZEXP
?.(   ZSINH ZSIN ZCOSH ZCOS ZTANH ZTAN ZASINH ZASIN) ?.cr

SET-ZT-MODE-ZEXACT
t{ 1e 0e 1/z -> 1e -0e }t

\ 1/(x+ieps) = 1/x - ieps/x^2 + ...
' 1/z is func
:noname  fnegate fswap 1/f fswap ; is gauge
SET-ZT-MODE-ZEXACT
z{  1e  0e }gauge
z{  1e -0e }gauge
z{ -2e  0e }gauge
z{ -2e -0e }gauge
z{ -1e -0e }gauge

:noname i* 1/z ; is func
:noname i* 1/f fnegate ; is gauge
SET-ZT-MODE-ZEXACT
z{  1e  0e }gauge
z{  1e -0e }gauge
z{ -2e  0e }gauge
z{ -2e -0e }gauge

\ sqrt(x+ieps) = sqrt(x) + ieps/sqrt(x) + ...
' zsqrt is func
:noname  fswap fsqrt fswap ; is gauge
SET-ZT-MODE-ZEXACT
z{ 0.5e  0e }gauge
z{ 0.5e -0e }gauge
z{  15e  0e }gauge
z{  15e -0e }gauge

\ ln(x+ieps) = ln(x) + ieps/x + ... 
' zln is func
:noname  fover f* fswap fln fswap ; is gauge
SET-ZT-MODE-ZEXACT
z{ 0.5e  0e }gauge
z{ 0.5e -0e }gauge
z{  15e  0e }gauge
z{  15e -0e }gauge

\ exp(x+ieps) = exp(x) + ieps*exp(x) + ...
' zexp is func
:noname fswap fexp fdup -frot f* ; is gauge 
SET-ZT-MODE-ZEXACT
z{  1e  0e }gauge
z{  1e -0e }gauge
z{ -2e  0e }gauge
z{ -2e -0e }gauge

\ sinh(x+ieps) = sinh(x) + ieps*cosh(x) + ...
' zsinh is func
\ :noname  fswap fsinh fswap ; is gauge
:noname fover fcosh f* fswap fsinh fswap ; is gauge 
SET-ZT-MODE-ZEXACT
z{  1e  0e }gauge
z{  1e -0e }gauge
z{ -2e  0e }gauge
z{ -2e -0e }gauge

\ sin(x+ieps) = sin(x) + ieps*cos(x) + ...
' zsin is func
:noname  fover fcos f* fswap fsin fswap ; is gauge
SET-ZT-MODE-ZEXACT
z{  1e   0e }gauge
z{  1e  -0e }gauge
z{ -0.5e 0e }gauge
z{ -2e  -0e }gauge

\ cosh(x+ieps) = cosh(x) + ieps*sinh(x) + ...
' zcosh is func
:noname fover fsinh f* fswap fcosh fswap ; is gauge 
SET-ZT-MODE-ZEXACT
z{  1e  0e }gauge
z{  1e -0e }gauge
z{ -2e  0e }gauge
z{ -2e -0e }gauge

\ cos(x+ieps) = cos(x) - ieps*sin(x) + ...
' zcos is func
:noname  fover fsin fnegate f* fswap fcos fswap ; is gauge
SET-ZT-MODE-ZEXACT
z{  1e   0e }gauge
z{  1e  -0e }gauge
z{ -0.5e 0e }gauge
z{ -2e  -0e }gauge

\ tanh(x+ieps) = tanh(x) + ieps/[cosh(x)]^2 + ...
' ztanh is func
:noname fover fcosh f^2 f/ fswap ftanh fswap ; is gauge 
SET-ZT-MODE-XREL SET-ZT-MODE-YEXACT
z{  1e  0e }gauge
z{  1e -0e }gauge
z{ -2e  0e }gauge
z{ -2e -0e }gauge

\ tan(x+ieps) = tan(x) + ieps/[cos(x)]^2 + ...
' ztan is func
:noname fover fcos f^2 f/ fswap ftan fswap ; is gauge 
SET-ZT-MODE-ZEXACT
z{  1e  0e }gauge
z{  1e -0e }gauge
z{ -2e  0e }gauge
z{ -2e -0e }gauge

\ asinh(x+ieps) = asinh(x) + ieps/sqrt(1+x^2) + ...
\ cut |y| >= 1 automatically avoided on real axis
' zasinh is func
:noname
  fover f^2 1e f+ fsqrt 1/f f* fswap fasinh fswap ; is gauge
SET-ZT-MODE-ZEXACT
z{  0.7e  0e }gauge
z{ -1.2e  0e }gauge
NOBLE? [IF] SET-ZT-MODE-YABS [THEN]
z{  0.7e -0e }gauge
z{ -1.2e -0e }gauge

\ asin(x+ieps) = asin(x) + ieps/sqrt(1-x^2) + ...
\ cut |x| >= 1 has to be avoided
' zasin is func
:noname
  fover f^2 1e f- fnegate fsqrt 1/f f* fswap fasin fswap ; is gauge
SET-ZT-MODE-ZEXACT
NOBLE? [IF] SET-ZT-MODE-XREL SET-ZT-MODE-YABS [THEN]
z{  0.7e  0e }gauge
z{  0.7e -0e }gauge
SET-ZT-MODE-XREL
z{ -0.2e  0e }gauge
z{ -0.2e -0e }gauge

testing principal branch cuts:  ZSQRT ZLN Z^ ZASINH ZACOSH
?.(   ZATANH ZACOTH ZASIN ZACOS ZATAN ZACOT) ?.cr
(
Branch cut gauges.  In the formulas that follow, eps will be
replaced by +0 or -0.  The sign of the floating-point zero
inputs in the corresponding words is significant.
)

: sqrt-cut  ( f: x 0 -- u+iv )
(
z=x+ieps,  x <= 0:
zsqrt[z] = i*sgn[eps]*sqrt[|x|]
)
  fswap fabs fsqrt fswap                ( f: sqrt[|x|] 0)
  fsgn0* 0e fswap ;                     ( f: 0 sgn[0]*sqrt[|x|])

' zsqrt    is func
' sqrt-cut is gauge
SET-ZT-MODE-ZEXACT
z{ -2e  0e }gauge
z{ -2e -0e }gauge
z{  0e  0e }gauge
z{  0e -0e }gauge
z{ -0e  0e }gauge
z{ -0e -0e }gauge

t{  0e  0e zsqrt -> 0e  0e }t
t{  0e -0e zsqrt -> 0e -0e }t
t{ -0e  0e zsqrt -> 0e  0e }t
t{ -0e -0e zsqrt -> 0e -0e }t

: ln-cut  ( f: x 0 -- u+iv )
(
z=x+ieps,  x <= 0
zln[z] = ln[|x|] + i*sgn[eps]*pi
)
  fswap fabs fln pi frot fsgn0* ;       ( f: ln|x| sgn[0]*pi)
' zln    is func
' ln-cut is gauge
SET-ZT-MODE-ZEXACT
z{ -2e  0e }gauge
z{ -2e -0e }gauge
z{ -0e  0e }gauge
z{ -0e -0e }gauge

t{  0e  0e zln -> -Inf  0e }t
t{  0e -0e zln -> -Inf -0e }t

SET-ZT-MODE-ZEXACT
t{ 0e 0e 0e 0e  z^ fnan? fnan? ->  true true }t
t{ 0e 0e 1e 0e  z^ fnan? fnan? ->  true true }t

\ zvariable a+ib
FVARIABLE ln|x|  FVARIABLE pi*sgn0  \ absent local fvariables...
FVARIABLE a      FVARIABLE b

: z^-cut  ( f: x+i0 a+ib -- u+iv )
(
z=x+ieps,  x <= 0
z^[a+ib] = exp[a*ln|x| - b*sgn[eps]*pi]
             * exp[i*a*[sgn[eps]*pi + b*ln|x|]]
)
  b f!   a f!
  pi fswap fsgn0* pi*sgn0 f!
  fabs fln ln|x| f!
  a f@   ln|x| f@   f*                  ( f: c=a*ln|x| )
  b f@   pi*sgn0 f@   f*                ( f: c d=sgn[0]*b*pi)
  f- fexp                               ( f: mod=exp[c-d])
  a f@   pi*sgn0 f@   f*                ( f: mod e=sgn[0]*a*pi)
  b f@   ln|x| f@   f* f+               ( f: mod ang=e+b*ln|x|)
  polar> ;

' z^     is func
' z^-cut is gauge
SET-ZT-MODE-ZEXACT
z{ -0.5e  0e  1e  0e     }2gauge
z{ -0.5e -0e  1e  0e     }2gauge
z{ -0.5e  0e  0e  1e     }2gauge
z{ -0.5e -0e  0e  1e     }2gauge
z{ -3.5e  0e -0.5e -0.5e }2gauge
z{ -3.5e -0e -0.5e -0.5e }2gauge
PFE? 0= [IF]
  1.5e-15 SET-REL-ERRS   \ 1.45e-15 fails
  SET-ZT-MODE-XREL SET-ZT-MODE-YREL [THEN]
z{ -1.5e  0e  0.5e  0.5e }2gauge
z{ -1.5e -0e  0.5e  0.5e }2gauge
rel-err-default SET-REL-ERRS

: asinh-cut  ( f: 0 y -- u+iv )
(
z=eps+iy,  |y| >= 1:
asinh[z] = ln[| y + sgn[eps*y]*sqrt[y^2-1] |] + i*sgn[y]*pi/2
)
  fdup frot fsgn0*                      ( f: y a=y*sgn[0])
  fsgn fover f^2 1e f- fsqrt f*         ( f: y b=sgn[a]*sqrt[y^2-1])
  fover f+ fabs fln fswap               ( f: u=ln[|y+b|] y)
  fsgn pi/2 f* ;                        ( f: u v=sgn[y]*pi/2)

' zasinh    is func
' asinh-cut is gauge
SET-ZT-MODE-ZEXACT
z{  0e  2e }gauge
z{  0e -2e }gauge
NOBLE? 0= [IF] SET-ZT-MODE-XREL [THEN]
z{ -0e  2e }gauge
z{ -0e -2e }gauge

: acosh-cut(-1,1)  ( f: x 0 -- u+iv )
(
z=x+ieps,  |x| <= 1:
acosh[z] = sgn[eps]*2i*atan[ sqrt[ [1-x]/[1+x] ] ]
)
  1e fswap fsgn0* fswap                 ( f: sgn[0] x)
  1e fover f- fswap 1e f+ f/            ( f: sgn[0] a=[1-x]/[1+x])
  fsqrt fatan f* f2* 0e fswap ;         ( f: 0 2*atan[sqrt[a]])

' zacosh          is func
' acosh-cut(-1,1) is gauge
SET-ZT-MODE-ZEXACT
SET-ZT-MODE-YREL  \ not necessary with OS X ppc
z{ -0.5e  0e }gauge
z{ -0.5e -0e }gauge
SET-ZT-MODE-YREL  \ not necessary with OS X intel
z{  0.5e  0e }gauge
z{  0.5e -0e }gauge

: acosh-cut(-Inf,-1)  ( f: x 0 -- u+iv )
(
z=x+ieps,  -Inf < x < -1:
acosh[z] = 2*ln[ [ sqrt[|x+1|] + sqrt[|x-1|] ] / sqrt[2] ]
             + i*sgn[eps]*pi
)
  pi fswap fsgn0* fswap                 ( f: v=i*pi*sgn[0] x)
  fdup 1e f+ fabs fsqrt
      fswap 1e f- fabs fsqrt f+         ( f: v a=sqrt[x+1]+sqrt[x-1])
  1/rt2 f* fln f2* fswap ;              ( f: u=2*ln[a/rt2] v)

' acosh-cut(-Inf,-1) is gauge
SET-ZT-MODE-XREL SET-ZT-MODE-YEXACT
IFORTH-HOST [IF] SET-ZT-MODE-YREL [THEN]
z{ -1.5e  0e }gauge
z{ -1.5e -0e }gauge

: atanh-cut  ( f: x 0 -- u+iv )
(
z=x+ieps,  |x| >= 1
atanh[z] = ln[| [1+x]/[1-x] |]/2 + i*sgn[eps]*pi/2
)
  pi/2 fswap fsgn0* fswap               ( f: v=i*pi/2*sgn[0] x)
  fdup 1e f+ fswap 1e f- f/ fabs        ( f: v a=|[1+x]/[1-x]|)
  fln f2/ fswap ;                       ( f: u=ln[a]/2 v)

' zatanh    is func
' atanh-cut is gauge
SET-ZT-MODE-ZEXACT
NOBLE? [IF] SET-ZT-MODE-XREL [THEN]
z{  1.5e  0e }gauge
z{  1.5e -0e }gauge
z{ -1.5e  0e }gauge
z{ -1.5e -0e }gauge

: acoth-cut  ( f: x 0 -- u+iv )
(
z=x+ieps,  |x| <= 1
acoth[z] = ln[| [1+x]/[1-x] |]/2 - i*sgn[eps]*pi/2
)
  ( f: 0) fnegate atanh-cut ;

NOBLE? [IF]  \ ZACOTH is defined
' zacoth    is func
' acoth-cut is gauge
SET-ZT-MODE-XREL SET-ZT-MODE-YEXACT
z{  0.5e  0e }gauge
z{  0.5e -0e }gauge
z{ -0.5e  0e }gauge
z{ -0.5e -0e }gauge
[ELSE]
?." ZACOTH not defined
[THEN]

: asin-cut  ( f: x 0 -- u+iv )
(
z=x+ieps,  |x| >= 1
asin[z] = sgn[x]*pi/2 - i*ln[| x - sgn[eps*x] * sqrt[x^2-1] |]
)
  fover fsgn pi/2 f* -frot              ( f: u=pi/2*sgn[x] x 0)
  fover fswap fsgn0* fsgn               ( f: u x a=sgn[0*x])
  fover f^2 1e f- fsqrt f*              ( f: u x b=a*sqrt[x^2-1])
  f- fabs fln fnegate ;                 ( f: u v=-ln|x+b|)

' zasin    is func
' asin-cut is gauge
SET-ZT-MODE-ZEXACT
NOBLE? 0= [IF] SET-ZT-MODE-YREL [THEN]
z{  1.5e  0e }gauge
z{  1.5e -0e }gauge
z{ -1.5e  0e }gauge
z{ -1.5e -0e }gauge

: acos-cut  ( f: x 0 -- u+iv )
(
z=x+ieps,  |x| >= 1
acos[z] = pi/2 - asin[z]
        = step[-x]*pi + i*ln[| x - sgn[eps*x] * sqrt[x^2-1] |]
)
  asin-cut znegate pi/2 x+ ;

' zacos    is func
' acos-cut is gauge
SET-ZT-MODE-ZEXACT
NOBLE? 0= [IF] SET-ZT-MODE-YREL [THEN]
z{  1.5e  0e }gauge
z{  1.5e -0e }gauge
z{ -1.5e  0e }gauge
z{ -1.5e -0e }gauge

: atan-cut  ( f: 0 y -- u+iv )
(
z=eps+iy,  |y| >= 1
atan[z] = sgn[eps]*pi/2 + i/2*ln[| [1+y]/[1-y] |]
)
  pi/2 frot fsgn0* fswap                ( f: u=sgn[0]*pi/2 y) 
  fdup 1e f+ fswap 1e f- f/
      fabs fln f2/ ;                    ( f: u v=ln|[1+y]/[1-y]|/2)

' zatan    is func
' atan-cut is gauge
\ set-znear
SET-ZT-MODE-ZEXACT
NOBLE? [IF] SET-ZT-MODE-YREL [THEN]
z{  0e  1.5e }gauge
z{ -0e  1.5e }gauge
z{  0e -1.5e }gauge
z{ -0e -1.5e }gauge

: acot-cut  ( f: 0 y -- u+iv )
(
z=eps+iy,  |y| <= 1
acot[z] = sgn[eps]*pi/2 + i/2*ln[ [1-y]/[1+y] ]
)
  atan-cut conjg ;

NOBLE? [IF] \ ZACOT is defined
' zacot    is func
' acot-cut is gauge
SET-ZT-MODE-XEXACT SET-ZT-MODE-YREL
z{  0e  0.5e }gauge
z{ -0e  0.5e }gauge
z{  0e -0.5e }gauge
z{ -0e -0.5e }gauge
[ELSE]
?." ZACOT not defined
[THEN]

[DEFINED] fcopysign [IF]
testing  FCOPYSIGN

SET-ZT-MODE-ZEXACT
t{   0e   0e fcopysign ->   0e }t
t{  -0e   0e fcopysign ->   0e }t
t{   0e  -0e fcopysign ->  -0e }t
t{  -0e  -0e fcopysign ->  -0e }t
t{   3e   0e fcopysign ->   3e }t
t{  -5e   0e fcopysign ->   5e }t
t{   7e  -0e fcopysign ->  -7e }t
t{ -11e  -0e fcopysign -> -11e }t
t{   0e   3e fcopysign ->   0e }t
t{  -0e   5e fcopysign ->   0e }t
t{   0e  -7e fcopysign ->  -0e }t
t{  -0e -11e fcopysign ->  -0e }t
t{   3e +Inf fcopysign ->   3e }t
t{  -5e +Inf fcopysign ->   5e }t
t{   7e -Inf fcopysign ->  -7e }t
t{ -11e -Inf fcopysign -> -11e }t
t{ +Inf   3e fcopysign -> +Inf }t
t{ -Inf   5e fcopysign -> +Inf }t
t{ +Inf  -7e fcopysign -> -Inf }t
t{ -Inf -11e fcopysign -> -Inf }t
t{ +Inf +Inf fcopysign -> +Inf }t
t{ -Inf +Inf fcopysign -> +Inf }t
t{ +Inf -Inf fcopysign -> -Inf }t
t{ -Inf -Inf fcopysign -> -Inf }t
[THEN]

[DEFINED] zbox [IF]
testing  ZBOX

SET-ZT-MODE-ZEXACT
t{   0e   0e zbox ->  1e   0e }t
t{  -0e   0e zbox -> -1e   0e }t
t{   0e  -0e zbox ->  1e  -0e }t
t{  -0e  -0e zbox -> -1e  -0e }t
t{ +Inf +Inf zbox ->  1e   1e }t
t{ -Inf +Inf zbox -> -1e   1e }t
t{ +Inf -Inf zbox ->  1e  -1e }t
t{ -Inf -Inf zbox -> -1e  -1e }t
t{ +Inf   3e zbox ->  1e   0e }t
t{ -Inf   3e zbox -> -1e   0e }t
t{ +Inf  -3e zbox ->  1e  -0e }t
t{ -Inf  -3e zbox -> -1e  -0e }t
t{   5e +Inf zbox ->  0e   1e }t
t{  -5e +Inf zbox -> -0e   1e }t
t{   5e -Inf zbox ->  0e  -1e }t
t{  -5e -Inf zbox -> -0e  -1e }t
t{ +Inf   0e zbox ->  1e   0e }t
t{ -Inf   0e zbox -> -1e   0e }t
t{ +Inf  -0e zbox ->  1e  -0e }t
t{ -Inf  -0e zbox -> -1e  -0e }t
t{   0e +Inf zbox ->  0e   1e }t
t{  -0e +Inf zbox -> -0e   1e }t
t{   0e -Inf zbox ->  0e  -1e }t
t{  -0e -Inf zbox -> -0e  -1e }t
t{   0e   5e zbox fnan? fnan? -> true true }t
t{  -0e   5e zbox fnan? fnan? -> true true }t
t{   0e  -5e zbox fnan? fnan? -> true true }t
t{  -0e  -5e zbox fnan? fnan? -> true true }t
t{   3e   0e zbox fnan? fnan? -> true true }t
t{  -3e   0e zbox fnan? fnan? -> true true }t
t{   3e  -0e zbox fnan? fnan? -> true true }t
t{  -3e  -0e zbox fnan? fnan? -> true true }t
t{   3e   5e zbox fnan? fnan? -> true true }t
t{  -3e   7e zbox fnan? fnan? -> true true }t
t{  11e  -3e zbox fnan? fnan? -> true true }t
t{ -17e -13e zbox fnan? fnan? -> true true }t
[THEN]

?.errors GFORTH-HOST [IF] ?.cr [THEN]