(       Title:  Borda's Mouthpiece
         File:  borda.fs
       Author:  David N. Williams
      Version:  0.5.1
      License:  Unrestrictable Public Domain, see UPDL.txt
Starting date:  March 27, 2004
Version 0.5.0:  March 30, 2004
Last revision:  April 13, 2004

ANSForth environmental dependences:

  1. Requires the Floating-Point, Floating-Point extension, and
     the Programming-Tools extension word sets.
  2. Assumes a separate floating point stack.
  3. Case insensitivity for standard words.
  4. The ANS FORTH default floats must be IEEE 754/854 to work
     correctly with signed zero/infinity.
  5. Other floating point technicalities if complex-kahan.fs is
     used.
)

\ *** PROBLEM
(
This example tests the correctness of signed zero handling on
branch cuts in our implementation of Kahan's algorithms for the
complex square root and natural logarithm [ZSQRT and ZLN], using
IEEE 754 floating point.  It is taken from "The Baleful Effect
of Computer Benchmarks upon Applied Mathematics, Physics and
Chemistry," by William Kahan, June 11, 1996:

  http://http.cs.berkeley.edu/~wkahan/ieee754status/baleful.ps.

To quote:

  Example:  Define complex analytic functions

    g[z] = z^2 + z sqrt[z^2 + 1],

  and

    F[z] = 1 + g[z] + log[g[z]].
  
  Plot the values taken by F[z] as z runs along eleven rays

    z =  r i,  z = r exp[4i pi/10], z = r exp[3i pi/10],
    z = r exp[2i pi/10], z = r exp[i pi/10], z = r

  and their complex conjugates, taking positive r from near zero
  to near infinity.

See the Kahan reference for plots of both incorrect and correct
results.  We expect only correct results here.

To get a printout suitable for cutting and pasting from the
screen to a LaTeX file for plotting with PSTricks, execute
.PLOTS after loading this file.
)

\ *** CONDITIONAL COMPILES

\ Only *one* of the following three [IF]'s should be true!

0 [IF]
  s" FLOATING-EXT" environment? [IF] ( flag) drop
    s" FLOATING-STACK" environment? [IF] ( maxdepth) drop
      s" COMPLEX-EXT" environment? [IF] ( version) drop
        cr .( Using pfe with COMPLEX-EXT.) cr
        : .provenance  ( -- )  cr
          ." % Produced with pfe's COMPLEX-EXT environment." ;
      [ELSE] .( ***COMPLEX-EXT not available.) cr bye [THEN]
    [ELSE] .( ***FLOATING-STACK not available.) cr bye [THEN]
  [ELSE] .( ***FLOATING-EXT not available.) cr bye [THEN]
[THEN]

1 [IF]
  s" complex.fs" included
  cr .( Using complex.fs.) cr
  PRINCIPAL-ARG 0= [IF]
    .( ***PRINCIPAL-ARG must be set to TRUE in complex.fs.) cr bye
  [THEN]
  : .provenance  ( -- )  cr
    ." % Produced with the ANSForth package complex.fs." ;
[THEN]

0 [IF]
  s" complex-kahan.fs" included
  cr .( Using complex-kahan.fs.) cr
  : .provenance  ( -- )  cr
    ." % Produced with the ANSForth package complex-kahan.fs." ;
[THEN]

s" [UN]" pad c! pad char+ pad c@ move 
pad find nip 0=
[IF]
: [UN]  ( "name" -- flag )
(
Leave true if name is in the search order, else leave false.
)
  bl word find nip 0= ; immediate
[THEN]

(
ZROT is defined in pfe's COMPLEX-EXT, but neither in complex.fs
nor in complex-kahan.fs.  But both of those have NONAME
temporary storage.  The code below knows that ZSWAP uses the
first of the three float slots at NONAME.
)
[UN] ZROT [IF]
  : zrot  ( f: z1 z2 z3 -- z2 z3 z1 )
    [ noname float+ ] literal z!
    zswap                               ( f: z2 z1)
    [ noname float+ ] literal z@        ( f: z2 z1 z3)
    zswap ;
[THEN]

3.141592653589793238463e fconstant pi

\ *** PSTRICKS CODE

: .plotheader  ( -- )
  cr  cr ." \savedata{\contourdata}[{%" ;

\ The spacing below is good for pfe's default precision, but not
\ for gforth's.
: .dataheader  ( f: sin[theta] -- )  \ just comments
(
The spacing here for x and y column heads is tuned for pfe's
default float output format.
)
  cr ." % angle: "
  ( sin) fasin pi f/ fs. ." pi"
  cr ." %      x              y" ;

: .blueplottrailer  ( -- )
  cr ." }]"
  cr ." \dataplot[plotstyle=curve,linecolor=blue]{\contourdata}" ;

: .redplottrailer  ( -- )
  cr ." }]"
  cr ." \dataplot[plotstyle=curve,linecolor=red]{\contourdata}" ;

\ *** CALCULATION
(
The complex lexicon uses the float stack, assumed to be separate
from the data stack, as the complex number stack, and notations
like
)
     ( f: z1 -- z2 )
(
are short for
)
    ( f: x1 y1 -- x2 y2 )

: g  ( f: z -- z^2+z*sqrt[z^2+1] )
  zdup z^2 zdup            ( f: z z^2 z^2)
  1e0 x+ zsqrt             ( f: z z^2 sqrt[z^2+1])
  zrot                     ( f: z^2 sqrt[z^2+1] z)
  z* z+ ;

: F ( f: z -- 1+g[z]+ln[g[z]] ) 
  g zdup zln z+ 1e x+ ;

(
The ray parameters R-MIN and R-MAX are tuned to make the 11
contours sufficiently overlap the plot region, defined below.
)
  8e-3 fconstant  r-min
  3e   fconstant  r-max

  50   constant   #points
  r-max r-min f- #points 1- 0 d>f f/
       fconstant  delta-r

(
Here are the plot clipping limits, based on the plot region in
Kahan's article mentioned above.  We use X-MAX instead of XMAX
to avoid name collision with complex-kahan.fs.  We clip the plot
output in this code, because we're not fluent enough to see how
to do it in PSTricks.
)
  -3.9e fconstant x-min  7.9e fconstant x-max
  -6.9e fconstant y-min  6.9e fconstant y-max

: inclip?  ( f: z -- s: flag )
(
Leave true if z is in the plot region.
)
  fdup y-min f> y-max f< and
  fdup x-min f> x-max f< and and ;

: xy.  ( f: x y -- )
(
Emit x and y in scientific notation.
)
  fswap
  space fdup f0< 0= IF space THEN fs.
  space fdup f0< 0= IF space THEN fs. ;

fvariable r  \ work around missing float locals

: .contour  ( f: cos[theta] sin[theta] -- )
(
Print out the data for the part of the F contour corresponding
to fixed theta, from R-MIN to R-MAX, but within the plot region.
)
  ( sin[theta]) fdup .dataheader
  r-min r f!
  #points 0 DO
    zdup r f@ fdup delta-r f+ r f!
    z*f F                     ( r: cos[theta] sin[theta] F[z])
    zdup inclip? IF cr xy. ELSE zdrop THEN
  LOOP ( f: cos sin) zdrop ;


: .plots  ( -- )
(
Run this to get the PSTricks code printout for insertion into
borda.tex.
)
  .provenance
  .plotheader
  z=1 .contour  \ only one middle ray 
  .blueplottrailer
  .1e            \ initial fraction of pi
  4 0 DO
    fdup pi f*                ( f: theta/pi theta)
    fsincos fswap             ( f: theta/pi cos[theta] sin[theta])
    .plotheader zdup    .contour .blueplottrailer
    .plotheader fnegate .contour .blueplottrailer
    .1e f+
  LOOP ( f: theta/pi) fdrop
  \ rays at plus or minus pi/2 are special
  z=i    .plotheader .contour .redplottrailer
  0e -1e .plotheader .contour .redplottrailer
  cr ;  \ cr to avoid possible trailing "ok"