(      Zelefunt:  Test complex exponential
           File:  zexp.fs
        Version:  0.9.3
Original Author:  W. J. Cody [November 29, 1990]
      Ported by:  David N. Williams
        License:  ACM noncommercial use
  Last revision:  January 19, 2005

This is a port of W. J. Cody's complex exponential test program
in the celefunt package, from Fortran 77 to ANS Forth.  As ACM
Algorithm 714 [TOMS], celefunt is presumed to be under the ACM
license for noncommercial use:

http://www.acm.org/pubs/copyright_policy/softwareCRnotice.html

Besides an ANS Forth environmental dependence on lower case,
this code uses DEFER and IS.

The code uses a complex number lexicon mostly compatible with
that of Everett Carter in the Forth Scientific Library [FSL],
with extensions by Julian Noble and myself [dnw].  In
particular, a separate floating point stack is assumed, on which
complex numbers like z=x+iy appear as [ x y], with y nearest the
top.

Quote from the original celefunt package:

   Program to test CEXP

   Accuracy tests are based on the identity

             exp[z] = exp[z-w]*exp[w]

   Data required:

      None

   Subprograms required from this package:

      MACHAR - An environmental inquiry program providing
               information on the floating point arithmetic
               system.  Note that the call to MACHAR can
               be deleted provided the following four
               parameters are assigned the values indicated:

               IBETA  - the radix of the floating point system;
               IT     - the number of base-IBETA digits in the
                        significand of a floating-point number;
               XMIN   - the smallest non-vanishing floating-point
                        power of the radix;
               XMAX   - the largest finite floating-point no.

      REN[K] - a function subprogram returning random real
               numbers uniformly distributed over [0,1]

      RESET, TABLAT, REPORT - programs to report results.

   Latest revision - October 31, 1990

   Author - W. J. Cody
            Argonne National Laboratory

Porting notes:

Version 0.9.3
19Jan05 * Added summary output lines using .SHORT-ZSTATS from
          zfunstat.fs version 0.9.3.

Version 0.9.2
 8Apr03  * Adjusted to work with machar.fs 0.9.5 and later,
           which defines machine parameters as constants.

Version 0.9.1
 5Apr03  * Start and finish.

Our Forth port of MACHAR uses #DIGITS for IT, and also returns the
floating point conversion of IBETA, called BETA.

The port of REN does not take an argument.  The argument K in
the original Fortran has no effect on the value returned by REN.

The port of RESET, TABLAT, and REPORT is contained in zfunstat.fs,
which also loads the ports of MACHAR and REN.
)

decimal
s" FLOATING-EXT"   environment? [IF] ( flag) drop
s" FLOATING-STACK" environment? [IF] ( maxdepth) drop

\ loadm complex  \ pfe
\ s" complex.fs"       included
s" complex-kahan.fs" included
s" zfunstat.fs"      included  \ includes machar.fs

cr .( RANDOM ARGUMENT ACCURACY TESTS)
cr
cr .#fdigits
cr

1e 16e f/ fdup               zconstant del
6.2416044877018563681e-2
    6.6487597751003112768e-2 zconstant expdel-1

3 set-precision
cr .( Testing exp[z] against exp[z-del]*exp[del], del = 1/16 + i*1/16.)
cr

:noname  ( -- )  zcurrent z@ zexp fcurrent z! ;
is function

:noname  ( -- )
  zcurrent z@ del z- zexp zdup expdel-1 z* z+ gcurrent z! ;
is gauge

' noop is purified
0e 1e 16e f/ corner1 z!   1e 1e corner2 z!
get-zstats .zstats
cr .( Summary:) .short-zstats  cr

cr .( NEW BOX)
cr
1.625e 1.625e corner1 z!   3e 3e corner2 z!
get-zstats .zstats
cr .( Summary:) .short-zstats  cr

cr .( NEW BOX)
cr
16e 16e corner1 z!   17e 17e corner2 z!
get-zstats .zstats
cr .( Summary:) .short-zstats  cr

cr .( SPECIAL TESTS)
cr
cr .( Test of the identity EXP[Z]*EXP[-Z] = 1:)
cr .(  Z) 34 spaces  .( EXP[Z]*EXP[-Z] - 1) 

6 set-precision

:noname  ( -- )
  5 0 DO cr
    beta ren f* beta ren f* zdup zs. 5 spaces
    zdup zexp zswap znegate zexp z* 1e x- zs. 
  LOOP cr
; execute

cr .( Special argument:  EXP[0+i0] - 1 = ) z=0 zexp 1e x- zs. 
cr
cr .( Test of cos[y], sin[y], and exp[x]:)
cr .( [REAL[EXP[0+i10]] - COS[10]]/COS[10] = )
    0e 10e zexp fover 10e fcos f- 10e fcos f/ fs.
cr .( [IMAG[EXP[0+i10]] - SIN[10]]/SIN[10] = )
    10e fsin f- 10e fsin f/ fs. fdrop
cr .(  [EXP[EXP[10+i0]] - EXP[10]]/EXP[10] = )
    10e 0e zexp fdrop 10e fexp f- 10e fexp f/ fs.
cr

5 set-precision
cr .( Tests near extreme arguments:)
cr
cr .( The real part of the first result below may underflow:)
cr .( EXP[ ) xmin fln
    fround    \ celefunt uses AINT() here, which truncates
    1e zdup zs. .( ] = ) zexp zs.
cr .( EXP[ ) xmax fln
    fround 1e f- \ celefunt uses AINT() here
    fdup 1e zdup zs. .( ] = ) zexp zs.
cr

cr .( There is an argument reduction error if)
cr .( the following are not about the same:)
cr .( EXP[ )
    (f: aint[ln[xmax]]) 0.5e f* 10e zdup zs. .( ]   = )
   zdup zexp zs.
cr .( EXP[ ) 0.5e z*f zdup zs. .( ]^2 = )
   zexp z^2 zs.
cr

3 set-precision
cr .( Test of error returns:)
cr

cr .( The real part of the argument below is so negative that the)
cr .( real exponential should underflow to zero, giving 0+i0:)
cr .( EXP[ ) -1e xmin fsqrt f/ 1e zdup zs. .( ] = ) zexp zs.
cr
cr .( The imaginary part of the argument below is too large for)
cr .( the sine and cosine computations to make sense:)
cr .( EXP[ ) 10e xmax fsqrt zdup zs. .( ] = ) zexp zs.
cr
cr .( The real part of the argument below is so large that the real)
cr .( exponential should overflow, giving Inf-i*Inf and/or an error:)
cr .( EXP[ ) 1e xmin fsqrt f/ -1e zdup zs. .( ] = ) zexp zs.
cr

cr .( END OF TESTS)
cr

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