(      Zelefunt:  Test complex power
           File:  zpow.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 18, 2005

This is a port of W. J. Cody's complex power 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 complex power [exponentiation]

   Accuracy tests are based on the identity

             z = z ** [1,0]

             z * z = Z ** [2,0]
      and
             z ** w = [z*z] ** [w/2]

   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;
               MAXEXP - the smallest integer power of FLOAT[IBETA]
                        that causes overflow;
               XMIN   - the smallest non-vanishing floating-point
                        power of the radix.

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

      RESET, TABLAT, REPORT - programs to report results.

   Latest revision - November 29, 1990

   Author - W. J. Cody
            Argonne National Laboratory

Porting notes:

Version 0.9.3
18Jan05 * 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
11Mar03  * Start.
 5Apr03  * 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  \ 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 the function z^[1+i0] against the gauge z.)
cr

:noname  ( -- )  zcurrent z@ z=1 z^ fcurrent z! ;
is function

:noname  ( -- )  zcurrent z@ gcurrent z! ;
is gauge

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

cr .( Testing the function z^[2+i0] against the gauge z*z.)
cr

:noname ( -- )  zcurrent z@ 2e 0e z^ fcurrent z! ;
is function

:noname ( -- )  zcurrent z@ zdup z* gcurrent z! ;
is gauge

' filtered1 is filtered
' purified1 is purified
1e 0e corner1 z!   10e 10e corner2 z!
get-zstats .zstats
cr .( Summary:) .short-zstats  cr

cr .( Testing the function z^w against the gauge [z*z]^[w/2].)
cr
cr .( *****************************************************)
cr .( * Very large MRE values are the norm for this test; *)
cr .( * the frequency of exact agreement is meaningful.   *)
cr .( *****************************************************)
cr

true to w?

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

:noname ( -- )
  zcurrent z@ zdup z*
  wcurrent z@ f2/ fswap f2/ fswap z^ gcurrent z!
; is gauge

' noop is filtered
4e 4e corner1 z!   10e 10e corner2 z!
get-zstats .zstats
cr .( Summary:) .short-zstats  cr

cr .( SPECIAL TESTS)
cr
cr .( Test of the identity Z^W = [1/Z]^[-W]:)
cr .(  Z) 26 spaces  .( W) 25 spaces .( [Z^W - [1/Z]^[-W]]/|Z^W|)  

:noname  ( -- )
  3 set-precision
  5 0 DO
    ren 10e f* ren 10e f*
    ren 10e f* ren 10e f* wcurrent z!
    cr zdup zs. 2 spaces wcurrent z@ zs. 2 spaces
    zdup wcurrent z@ z^                 (f: z f=z^w)
    zswap 1/z wcurrent z@ znegate z^    (f: f g=[1/z]^[-w])
    zover z- znegate zswap              (f: f-g f)             
    |z| z/f zs.     
  LOOP cr
; execute

6 set-precision
cr .( Test near an extreme argument.  This should not, but may)
cr .( trigger an error message:)
cr .( [BETA + i*BETA*XMIN]^[MAXEXP-1+i0]        = )
   beta beta xmin f*   maxexp 1- s>f 0e z^ zdup zs.
\ Next is BETA factored out times the first two terms
\ of the binomial expansion:
cr .( BETA^[MAXEXP-1] * [1 + i*[MAXEXP-1]*XMIN] = )
   beta maxexp 1- dup f^n
   1e s>f xmin f* frot z*f zdup zs. 
cr .( difference                                = )
   z- zs.
cr
cr .( Test of error returns:)
cr .( [0+i0]^[0+i0] = ) 0e 0e zdup z^ zs.
cr
cr .( END OF TESTS)
cr

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