( Title: Borda's Mouthpiece
File: borda.fs
Author: David N. Williams
License: Unrestrictable Public Domain, see UPDL.txt
Version: 0.5.2
Last revision: January 11, 2009
Version 0.5.2
11Jan09 * Updated to work with complex.fs 0.8.5:
old new
z=1 1e 0e
z=i 0e 1e
Version 0.5.1
13Apr04 * Revised stack comments and .CONTOUR.
Version 0.5.0
27Mar04 * Started.
30Mar04 * Last revision.
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 ( f: 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
1e 0e .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
0e 1e .plotheader .contour .redplottrailer
0e -1e .plotheader .contour .redplottrailer
cr ; \ cr to avoid possible trailing "ok"