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

This is a port of W. J. Cody's CLOG 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 CLOG

  Accuracy tests are based on the identity

            log[z] = log[z*z]/2
     and
            log[z] = log[17z/16] - log[17/16]

  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 - November 20, 1990

  Author - W. J. Cody
           Argonne National Laboratory

Porting notes:

Version 0.9.3
17Jan05 * 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  * Cosmetic changes to align with version 0.9.1 of
           zfunstat.fs.

Version 0.9.0
20Jan03  * Start.
11Mar03  * Last revision.

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  \ be sure to make PRINCIPAL-ARG true
s" complex-kahan.fs" included
s" zfunstat.fs"      included  \ includes machar.fs

6 set-precision

cr .( RANDOM ARGUMENT ACCURACY TESTS)
cr
cr .#fdigits
cr
cr .( Testing function ln[z] against gauge ln[z*z]/2.)
cr

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

:noname  ( -- )  zcurrent z@ z^2 zln z2/ gcurrent z! ;
is gauge

2e 0e corner1 z!   10e 10e corner2 z!
get-zstats .zstats
cr .( Summary:) .short-zstats  cr

cr .( NEW BOX)
cr
1000e -1000e corner1 z!   2000e -4000e corner2 z!
get-zstats .zstats
cr .( Summary:) .short-zstats  cr

cr .( NEW BOX)
cr
eps eps fnegate corner1 z!   0.25e -0.25e corner2 z!
get-zstats .zstats
cr .( Summary:) .short-zstats  cr

cr .( SPECIAL TESTS)
cr

:noname  ( -- )
  cr ." Test of the identity LN[Z] = -LN[1/Z]:"
  cr ."   Z" 34 spaces ." LN[Z] + LN[1/Z]"
  5 0 DO
    ren 10e f* ren 10e f*
    cr zdup zs. 5 spaces
      zdup zln zswap 1/z zln z+ zs.
  LOOP cr
; execute

cr .( Test of special argument:)
cr .( LN[ 1.0 + i0.0 ] = ) z=1 zln zs.
cr

cr .( Tests near extreme arguments.  These should not, but may)
cr .( trigger error messages:)
cr .( LN[ XMIN + i*XMIN ] = LN[ ) xmin xmin zdup zs. .( ])
cr 20 spaces .( = ) zln zs.
cr .( LN[ XMAX + i*XMAX ] = LN[ ) xmax xmax zdup zs. .( ])
cr 20 spaces .( = ) zln zs.
cr
cr .( Test of error returns:)
cr .( LN[ 0.0 + i0.0 ] = ) z=0 zln zs.
cr
cr .( END OF TESTS)
cr

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