(       Title:  Complex Word Set Signed Zero,
                Signed Infinity, and NaN Tests
         File:  complex-ieee-test.fs
       Author:  David N. Williams
      Version:  1.0.2
      License:  LGPL
Last revision:  November 29, 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.2
17Nov10 * Removed PRINCIPAL-ARG warning for complex.fs.
29Nov10 * Renamed ISNAN and COPYSIGN as FNAN? and FCOPYSIGN.

Version 1.0.1
18Aug09 * Added report for separate/integrated fp stack.
	* Made all reports conditional on VERBOSE, except
	  ABORT's.
	* Revised library selection logic.
	* Removed defs of [UNDEFINED] and [DEFINED], now in
	  Forth 200x.

Version 1.0
18May09 * Retired complex-ieee-test.fs.  Renamed
	  complex-ieee-ttest.fs as complex-ieee-test.fs.
	* Revised conditional compilation logic.
	* Added conditionally compiled ISNAN, revised NaN tests
	  to use it, removed +NaN, added ISNAN tests.
	* Renamed Inf as +Inf.

Version 0.9.8
 6Mar09 * Temporarily renamed as "complex-ieee-ttest.fs". 
	  Adapted to use Anton Ertl's ttester.fs instead of
	  tester.fs and ftester.fs.
	* Removed "(f:" in favor of "( f:".
	* Added tests of Z^ with NaN ouputs. 
	* Added disclaimer about testing for accuracy.
14Mar09 * Added conditional definition of NaN.  Revised
	  conditional definition of Inf and -Inf.
15Mar09 * Forced output of Z^ to positive for 0+i0^0+i0 and
	  0+i0^1+i0, to fix failure on intel found by Ed Falat.
	  This leaves unconfronted a "feature" whereby gcc
	  returns opposite NaN signs for the same arithmetic
	  operation on ppc and intel systems with pfe and
	  gforth.
	* Renamed NaN as +NaN.

Version 0.9.7
21Sep08 * Renamed F-ROT as -FROT.

Version 0.9.6
 6Sep06 * Really changed names of mixed argument words, not
	  actually done in 0.9.5.

Version 0.9.5
24Apr05 * Added conditional definition of constants Inf and
          -Inf.
        * Changed names of mixed argument, minimal computation
          words to the more telegraphic, jvn style.

Version 0.9.4
 7Feb05 * Renamed complex-szero-test.fs as complex-ieee-test.fs,
          and started moving all tests with +Inf, -Inf, and NaN
          from complex-test.fs to here.
 8Feb05 * Finished moving tests.
12Feb05 * Reorganized and added signed zero tests on the real
          axis for ZEXP, ZCOSH, ZCOS, ZTANH, ZTAN.
        * Removed all references to Z=0, etc., in favor of
          0E 0E, etc.
 4Mar05 * Added switch to skip a few tests with complex.fs.

Version 0.9.3
12Jan05 * Added ZF*, ZF/, FZ*, FZ/, ZIF*, ZIF/, IFZ*, IFZ/,
          Z2*, Z2/.
        * Added #ERRORS output, requiring ftester.fs 1.1.3 or
          above.

Version 0.9.2
27Feb03 * Started with signed zero and branch cut definitions
          from complex-test.fs.
 5Mar03 * 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 single or double floating point systems.

It is intended to test for formal correctness, not high
accuracy.

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 here, and then it's only a few.  [Search on
   occurrences of SKIP-TEST below.]  We're not sure it's
   worthwhile to be more aggressive with Kahan algorithms there.
   
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 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.

Except for DEFER and IS, this code is compatible with ANS Forth,
with an environmental dependence on lower case.
)

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

false value SKIP-TEST  \ for a few failures in complex.fs

\ USER-CONFIG: change to select a different library implementation
1 constant COMPLEX-LIB#

COMPLEX-LIB# 1 = [IF]
  s" complex.fs" included  \ defines PRINCIPAL-ARG  
  true to SKIP-TEST
  true to PRINCIPAL-ARG
[THEN]

COMPLEX-LIB# 2 = [IF]
  s" complex-kahan.fs" included
[THEN]

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

s" ttester.fs" included
true verbose !
decimal

: ?s.       ( c-addr len -- )  verbose @ IF type ELSE 2drop THEN ;
: ?emit-cr  ( -- )             verbose @ IF cr THEN ;

:NONAME  ( lib.case -- )
  COMPLEX-LIB# CASE
    1 OF s" Testing ieee 754 special values for complex.fs." ENDOF
    2 OF s" Testing ieee 754 special values for complex-kahan.fs." ENDOF
    3 OF s" Testing ieee 754 special values for pfe complex words." ENDOF
    ABORT" ***No complex words loaded"
  ENDCASE ?emit-cr ?s. ; EXECUTE

verbose @ [IF]
:noname  ( -- fp.separate? )
  depth >r 1e depth >r fdrop 2r> = ; execute
cr .( Floating-point stack is )
[IF] .( *separate*) [ELSE] .( *not separate*) [THEN]
.(  from the data stack.) cr
[THEN]

: near-defaults  ( -- )
  1E-15 rel-near f!
  1E-15 abs-near f! ;
near-defaults

variable #errors   0 #errors !

:noname  ( c-addr u -- )
(
Display an error message followed by the line that had the
error.
)
  1 #errors +! error1 ; error-xt !

decimal

[UNDEFINED] -frot [IF]
  : -frot  frot frot ; [THEN]

[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]
  s" ADDRESS-UNIT-BITS" environment? 0=
  [IF] cr .( ***Can't determine ADDRESS-UNIT-BITS.) ABORT [THEN]
  1 floats * 64 = constant 64BIT-FLOATS immediate

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

  hex
  : fnan?  ( f: f -- s: nan? )  \ either quiet or signaling
    fabs
    64BIT-FLOATS [IF]
      df>raw64 0 7FF00000 d>
    [ELSE]    
      sf>raw32 7F800000 >
    [THEN] ;
  decimal
[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]

 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

\ comparison type enumeration
  0 constant znear
  1 constant xexact
  2 constant yexact
  3 constant zexact

znear value z-comparison-type 

: set-znear   znear  to z-comparison-type ;
: set-xexact  xexact to z-comparison-type ;
: set-yexact  yexact to z-comparison-type ;
: set-zexact  zexact to z-comparison-type ;

: type-error  ( c-addr u -- )
(
Display an error message followed by the line that had the
error.
)
  1 #errors +!
  type source type cr ;

: ?near  ( f: a b -- )
  fnearly= 0=
  IF s" incorrect near float result: " type-error THEN ;

: ?exact  ( f: a b -- )
  fexactly= 0=
  IF s" incorrect exact float result: " type-error THEN ;

: actual-x@  ( f: -- x )  actual-fresults faligned float+ f@ ;
: actual-y@  ( f: -- y )  actual-fresults faligned f@ ;

: z{  ( -- )  t{ ;  \ ttester stack cursors not actually used

: }z  ( -- )
(
Issue an error message if the float stack [expected] contents
and the saved [actual] contents do not both have depths of two,
or if the expected and actual contents of the data stack do not
have depth zero.  Compare the float stack contents using
FEXACTLY= and/or FNEARLY= according to the comparison type
enumeration in Z-COMPARISON-TYPE, and issue an error message if
there is not agreement.
)
\  depth actual-depth @ or
\  IF s" data stack not empty " type-error THEN
  fdepth actual-fdepth @ =   fdepth 2 = and
  IF  \ if depths are both 2
    fswap ( f: expected-y expected-x)
    z-comparison-type
    CASE
    znear  OF actual-x@ ?near   actual-y@ ?near  ENDOF
    xexact OF actual-x@ ?exact  actual-y@ ?near  ENDOF
    yexact OF actual-x@ ?near   actual-y@ ?exact ENDOF
    zexact OF actual-x@ ?exact  actual-y@ ?exact ENDOF
    cr ." }Z:  comparison type not set" cr
    ENDCASE
  ELSE \ depths not both 2
    s" wrong number of float results: " error
  THEN ;

defer gauge
defer func

zvariable zatemp   zvariable zbtemp

: ?gauge  ( f: z -- )
(
Compare the functions whose xt's are in FUNC and GAUGE.
) 
  zdup zatemp z! gauge   ->   zatemp z@ func  ;

: ?2gauge  ( f: z1 z2 -- )
(
Same as above with 2 complex arguments.  Note that the FUNC and
GAUGE results are still a single complex number.
) 
  zbtemp z! zatemp z!
  zatemp z@ zbtemp z@ gauge -> zatemp z@ zbtemp z@ func ;

testing  F~  ZNEGATE  CONJG  I*  -I*

set-exact
t{ -0E  0E  0E f~ -> false }t
t{ -0E -0E  0E f~ -> true  }t

set-exact
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-zexact
z{  0E  0E conjg ->  0E -0E }z
z{ -0E  0E conjg -> -0E -0E }z
z{  0E -0E conjg ->  0E  0E }z
z{ -0E -0E conjg -> -0E  0E }z

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

z{  0E  0E -i* ->  0E -0E }z
z{ -0E  0E -i* ->  0E  0E }z
z{  0E -0E -i* -> -0E -0E }z
z{ -0E -0E -i* -> -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-exact
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-exact
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-exact
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-exact
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-exact
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
s"   ZASINH  ZASIN ZATANH ZATAN" ?s. ?emit-cr

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

set-zexact
:noname  ; is gauge

' zsinh is func
z{  0E  0E  ?gauge }z
z{  0E -0E  ?gauge }z
z{ -0E  0E  ?gauge }z
z{ -0E -0E  ?gauge }z

' zsin is func
set-zexact
z{  0E  0E  ?gauge }z
z{  0E -0E  ?gauge }z
z{ -0E  0E  ?gauge }z
z{ -0E -0E  ?gauge }z

' ztanh is func
set-zexact
z{  0E  0E  ?gauge }z
z{  0E -0E  ?gauge }z
z{ -0E  0E  ?gauge }z
z{ -0E -0E  ?gauge }z

' ztan is func
set-zexact
z{  0E  0E  ?gauge }z
z{  0E -0E  ?gauge }z
z{ -0E  0E  ?gauge }z
z{ -0E -0E  ?gauge }z

' zasinh is func
set-znear
z{  0E  0E  ?gauge }z
z{  0E -0E  ?gauge }z
z{ -0E  0E  ?gauge }z
z{ -0E -0E  ?gauge }z

' zasin is func
set-znear
z{  0E  0E  ?gauge }z
z{  0E -0E  ?gauge }z
z{ -0E  0E  ?gauge }z
z{ -0E -0E  ?gauge }z

' zatanh is func
set-znear
z{  0E  0E  ?gauge }z
z{  0E -0E  ?gauge }z
z{ -0E  0E  ?gauge }z
z{ -0E -0E  ?gauge }z

' zatan is func
set-znear
z{  0E  0E  ?gauge }z
z{  0E -0E  ?gauge }z
z{ -0E  0E  ?gauge }z
z{ -0E -0E  ?gauge }z

testing real analyticity:  1/Z ZSQRT ZLN ZEXP
s"   ZSINH ZSIN ZCOSH ZCOS ZTANH ZTAN ZASINH ZASIN" ?s. ?emit-cr
set-zexact
t{ 1E 0E   1/z   ->  1E -0E  }t

set-yexact  \ to test signed zero at real analytic points

\ 1/(x+ieps) = 1/x - ieps/x^2 + ...
' 1/z is func
:noname  fnegate fswap 1/f fswap ; is gauge
z{  1E  0E ?gauge }z
z{  1E -0E ?gauge }z
z{ -2E  0E ?gauge }z
z{ -2E -0E ?gauge }z

z{ -1E -0E ?gauge }z

:noname i* 1/z ; is func
:noname i* 1/f fnegate ; is gauge
z{  1E  0E ?gauge }z
z{  1E -0E ?gauge }z
z{ -2E  0E ?gauge }z
z{ -2E -0E ?gauge }z

\ sqrt(x+ieps) = sqrt(x) + ieps/sqrt(x) + ...
' zsqrt is func
:noname  fswap fsqrt fswap ; is gauge
z{ 0.5E  0E ?gauge }z
z{ 0.5E -0E ?gauge }z
z{  15E  0E ?gauge }z
z{  15E -0E ?gauge }z

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

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

\ 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 
z{  1E  0E ?gauge }z
z{  1E -0E ?gauge }z
z{ -2E  0E ?gauge }z
z{ -2E -0E ?gauge }z

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

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

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

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

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

\ 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
z{  0.7E  0E ?gauge }z
z{ -1.2E  0E ?gauge }z
SKIP-TEST [IF]
s" Skipping two real analyticity tests for ZASINH." ?s. ?emit-cr
[ELSE]
z{  0.7E -0E ?gauge }z
z{ -1.2E -0E ?gauge }z
[THEN]

\ 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
SKIP-TEST [IF]
s" Skipping two real analyticity tests for ZASIN." ?s. ?emit-cr
[ELSE]
z{  0.7E  0E ?gauge }z
z{ -0.2E  0E ?gauge }z
[THEN]
z{  0.7E -0E ?gauge }z
z{ -0.2E -0E ?gauge }z

testing principal branch cuts:  ZSQRT ZLN Z^ ZASINH ZACOSH
s"   ZATANH ZACOTH ZASIN ZACOS ZATAN ZACOT" ?s. ?emit-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+i0,  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|])

: ln-cut  ( f: x 0 -- u+iv )
(
z=x+i0,  x <= 0
zln[z] = ln[|x|] + i*sgn[eps]*pi
)
  fswap fabs fln pi frot fsgn0* ;       ( f: ln|x| sgn[0]*pi)

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

: z^-cut  ( f: x+iy a+ib -- u+iv )
(
z=x+i0,  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> ;

: asinh-cut  ( f: 0 y -- u+iv )
(
z=0+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)

: acosh-cut(-1,1)  ( f: x 0 -- u+iv )
(
z=x+i0,  |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]])

: acosh-cut(-Inf,-1)  ( f: x 0 -- u+iv )
(
z=x+i0,  -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)

: atanh-cut  ( f: x 0 -- u+iv )
(
z=x+i0,  |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)

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

: asin-cut  ( f: x 0 -- u+iv )
(
z=x+i0,  |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|)

: acos-cut  ( f: x 0 -- u+iv )
(
z=x+i0,  |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+ ;

: atan-cut  ( f: 0 y -- u+iv )
(
z=0+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)

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

' zsqrt    is func
' sqrt-cut is gauge
set-znear
z{ -2E  0E ?gauge }z
z{ -2E -0E ?gauge }z
z{  0E  0E ?gauge }z
z{  0E -0E ?gauge }z
z{ -0E  0E ?gauge }z
z{ -0E -0E ?gauge }z
set-exact
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

' zln    is func
' ln-cut is gauge
set-znear
z{ -2E  0E ?gauge }z
z{ -2E -0E ?gauge }z
z{ -0E  0E ?gauge }z
z{ -0E -0E ?gauge }z
set-exact
t{  0E  0E zln -> -Inf  0E }t
t{  0E -0E zln -> -Inf -0E }t

set-zexact
t{ 0E 0E 0E 0E  z^ fnan? fnan? ->  true true }t
t{ 0E 0E 1E 0E  z^ fnan? fnan? ->  true true }t

' z^     is func
' z^-cut is gauge
set-znear
z{ -0.5E  0E  1E  0E     ?2gauge }z
z{ -0.5E -0E  1E  0E     ?2gauge }z
z{ -0.5E  0E  0E  1E     ?2gauge }z
z{ -0.5E -0E  0E  1E     ?2gauge }z
z{ -1.5E  0E  0.5E  0.5E ?2gauge }z
z{ -1.5E -0E  0.5E  0.5E ?2gauge }z
z{ -3.5E  0E -0.5E -0.5E ?2gauge }z
z{ -3.5E -0E -0.5E -0.5E ?2gauge }z

' zasinh    is func
' asinh-cut is gauge
set-znear
z{  0E  2E ?gauge }z
z{ -0E  2E ?gauge }z
z{  0E -2E ?gauge }z
z{ -0E -2E ?gauge }z

' zacosh          is func
' acosh-cut(-1,1) is gauge
set-znear
z{  0.5E  0E ?gauge }z
z{  0.5E -0E ?gauge }z
z{ -0.5E  0E ?gauge }z
z{ -0.5E -0E ?gauge }z
' acosh-cut(-Inf,-1) is gauge
z{ -1.5E  0E ?gauge }z
z{ -1.5E -0E ?gauge }z

' zatanh    is func
' atanh-cut is gauge
set-znear
z{  1.5E  0E ?gauge }z
z{  1.5E -0E ?gauge }z
z{ -1.5E  0E ?gauge }z
z{ -1.5E -0E ?gauge }z

[DEFINED] zacoth [IF]
' zacoth    is func
' acoth-cut is gauge
set-znear
z{  0.5E  0E ?gauge }z
z{  0.5E -0E ?gauge }z
z{ -0.5E  0E ?gauge }z
z{ -0.5E -0E ?gauge }z
[THEN]

' zasin    is func
' asin-cut is gauge
set-znear
z{  1.5E  0E ?gauge }z
z{  1.5E -0E ?gauge }z
z{ -1.5E  0E ?gauge }z
z{ -1.5E -0E ?gauge }z

' zacos    is func
' acos-cut is gauge
set-znear
z{  1.5E  0E ?gauge }z
z{  1.5E -0E ?gauge }z
z{ -1.5E  0E ?gauge }z
z{ -1.5E -0E ?gauge }z

' zatan    is func
' atan-cut is gauge
set-znear
z{  0E  1.5E ?gauge }z
z{ -0E  1.5E ?gauge }z
z{  0E -1.5E ?gauge }z
z{ -0E -1.5E ?gauge }z

[DEFINED] zacot [IF]
' zacot    is func
' acot-cut is gauge
set-znear
z{  0E  0.5E ?gauge }z
z{ -0E  0.5E ?gauge }z
z{  0E -0.5E ?gauge }z
z{ -0E -0.5E ?gauge }z
[THEN]

[DEFINED] fcopysign [IF]
testing  FCOPYSIGN

set-exact
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-exact
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]

verbose @ [IF] .( #ERRORS: ) #errors @ . cr [THEN]