( Title: Collect and report statistics for complex function comparisons File: zfunstat.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 RESET, TABLAT, and REPORT subprograms in the celefunt package, from Fortran 77 to ANS Forth. As ACM Algorithm 714 [TOMS], the celefunt package 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 dependency on lower case, this code uses DEFER and IS. It also uses our port of Cody's MACHAR and REN subprograms. To quote from celefunt: 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] Porting notes: Version 0.9.3 17Jan05 * Did some factoring. Added .SHORT-ZSTATS. 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 3Apr03 * Renamed PURIFIED as PURIFIED1, and made PURIFIED an execution vector, defaulting to PURIFIED1. * Added NOOP. * Added execution vector FILTERED, defaulting to NOOP. Added FILTERED1 for use with complex power tests. * Defaulted FUNCTION and GAUGE to NOOP. * Rewrote SET-ZCURRENT to accomodate FILTERED. * ADDED zvariable WCURRENT for w arguments, and >WMAX in error accumulation structure. 4Apr03 * Factored .STATS and made other cosmetic changes. * Added ?SET-WCURRENT, depending on W?. 5Apr03 * Changed names GET-STATS and .STATS to GET-ZSTATS and .ZSTATS. Version 0.9.0 20Jan03 * Start. 15Feb03 * 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. We use x and y for the real and imaginary parts of z. A ".s" suffix in stack comments is short for the ANS Forth pair [addr len] string representation. ) decimal s" FLOATING-EXT" environment? [IF] ( flag) drop s" FLOATING-STACK" environment? [IF] ( maxdepth) drop s" machar.fs" included s" ren.fs" included s" [UN]" pad c! pad char+ pad c@ move pad find nip 0= [IF] : [UN] ( "name" -- flag ) bl word find nip 0= ; immediate [THEN] [UN] (f: [IF] : (f: ( -- ) postpone ( ; immediate [THEN] [UN] fnip [IF] : fnip fswap fdrop ; [THEN] [UN] f>s [IF] : f>s (f: r -- s: n ) f>d drop ; [THEN] [UN] s>f [IF] : s>f ( n -- f: r ) 0 d>f ; [THEN] [UN] f0> [IF] : f0> (f: x y -- s: flag ) 0E fswap f< ; [THEN] [UN] f0<> [IF] : f0<> (f: x y -- s: flag ) f0= 0= ; [THEN] [UN] f2dup [IF] : f2dup (f: x y -- x y x y ) fover fover ; [THEN] [UN] cell [IF] : cell ( -- /cell ) 1 cells ; [THEN] [UN] float [IF] : float ( -- /cell ) 1 floats ; [THEN] [UN] complex [IF] : complex ( -- /complex ) 2 floats ; [THEN] [UN] z0= [IF] : z0= (f: z -- s: flag ) f0= f0= and ; [THEN] [UN] z. [IF] : z. (f: x y -- ) fswap f. ." + i " f. ; [THEN] [UN] zs. [IF] : zs. (f: x y -- ) fswap fs. ." + i " fs. ; [THEN] [UN] f2/ [IF] : f2/ (f: x -- x/2 ) 0.5e f* ; [THEN] [UN] f^n [IF] : f^n ( u f: x -- x^u ) 1E ( u) 0 ?DO fover f* LOOP fnip ; [THEN] [UN] x@ [IF] : x@ ( zaddr -- f: x ) f@ ; [THEN] [UN] y@ [IF] : y@ ( zaddr -- f: y ) float+ f@ ; [THEN] [UN] x! [IF] : x! ( zaddr f: x -- ) f! ; [THEN] [UN] y! [IF] : y! ( zaddr f: y -- ) float+ f! ; [THEN] [UN] noop [IF] : noop ( -- ) ; [THEN] [UN] \\ [IF] : \\ BEGIN -1 parse 2drop refill 0= UNTIL ; [THEN] beta #digits 2/ 4 + f^n fconstant scale beta fln fconstant albeta #digits s>f fconstant ait 2000 constant #args -999e fconstant all9 : ln_beta (f: f -- ln_beta[f] ) fdup f0= IF fdrop all9 ELSE fln albeta f/ THEN ; : field: ( offset size -- offset+size ) create over , + DOES> ( struc -- struc+offset ) @ + ; \ Error accumulation structure/union 0 complex field: >zmax complex field: >wmax complex field: >fmax complex field: >gmax float field: >max-err float field: >sum-squares cell field: >#f.ne.g ( size) dup constant /vec-err cell - \ overlap onto >#F.NE.G cell field: >#f.lt.g cell field: >#f.gt.g constant /part-err : mre ( err.struc -- f: max.rel.err ) >max-err f@ ; : rms-err ( err.struc -- f: rms.err ) >sum-squares f@ #args s>f f/ fsqrt ; : ulp (f: error -- ulp ) ln_beta ait f+ 0E (f: #digs+ln_beta[err] 0) fmax ; : mre-ulp ( err.struc -- f: ulp ) mre ulp ; : rms-ulp ( err.struc -- f: ulp ) rms-err ulp ; : vec-#same ( vec.err.struc -- #same ) >#f.ne.g @ #args swap - ; : part-#same ( part.err.struc -- #same ) dup >#f.lt.g @ swap >#f.gt.g @ + #args swap - ; falign here /vec-err allot value vec-err falign here /part-err allot value real-err falign here /part-err allot value imag-err zvariable zcurrent zvariable wcurrent false value w? zvariable fcurrent zvariable gcurrent defer function ' noop is function defer gauge ' noop is gauge defer filtered ' noop is filtered defer purified \ the default is PURIFIED1, see below \ opposite corners for box of random z's zvariable corner1 \ AX + iAY zvariable corner2 \ BX + iBY fvariable deltax fvariable deltay : reset-error ( -- ) 0e vec-err >max-err f! 0e vec-err >sum-squares f! 0 vec-err >#f.ne.g ! 0e real-err >max-err f! 0e real-err >sum-squares f! 0e imag-err >max-err f! 0e imag-err >sum-squares f! 0 real-err >#f.lt.g ! 0 real-err >#f.gt.g ! 0 imag-err >#f.lt.g ! 0 imag-err >#f.gt.g ! ; : update-error ( err.struc f: |err| -- ) fdup dup mre f> IF fdup dup >max-err f! zcurrent z@ dup >zmax z! wcurrent z@ dup >wmax z! \ does no harm when W? is false fcurrent z@ dup >fmax z! gcurrent z@ dup >gmax z! THEN f^2 >sum-squares dup f@ f+ f! ; : handle-vector-error ( -- ) ( We don't understand the celefunt notion of error when the function is zero, but we do it that way for comparison. Note that an accurate version of |Z| should be used here. I really think the equality test shouldn't use |Z| at all. ) fcurrent z@ z0= IF 1e 0e (f: z=1) ELSE fcurrent z@ gcurrent z@ z- fcurrent z@ z/ THEN |z| fdup f0<> IF 1 vec-err >#f.ne.g +! THEN (f: |err|) vec-err update-error ; : handle-part-error ( &fpart &gpart part.err.struc -- ) locals| part.err gpart fpart | fpart f@ f0= IF 1e ELSE fpart f@ gpart f@ f- fpart f@ f/ THEN fdup f0< IF 1 part.err >#f.lt.g +! ELSE fdup f0> IF 1 part.err >#f.gt.g +! THEN THEN fabs part.err update-error ; 0.95e fconstant c1 1.05e fconstant c2 : filtered1 (f: z -- z' ) ( If c1 < |y/x| < c2 is true, reselect y at random until false. ) BEGIN zdup fswap f/ fabs (f: |y/x|) fdup c1 f> (f: |y/x|) c2 f< and WHILE (f: y) fdrop deltay f@ ren f* corner1 y@ f+ (f: x deltay*ren+AY) REPEAT ; : purified1 (f: f -- [f*scale+f]-f*scale ) ( Chop the last few bits to make a representable number. ) fdup scale f* fdup frot (f: f*scale f*scale f) f+ fswap f- ; ' purified1 is purified : set-zcurrent (f: XL AY -- ) ( Select a z at random and process it with the execution vectors FILTERED and PURIFIED before storing it in ZCURRENT. ) fswap \ to use the same order of randoms as celefunt ren deltax f@ f* f+ fswap (f: x=ren*deltax+XL AY) ren deltay f@ f* f+ (f: x y=ren*deltay+AY) filtered purified zcurrent y! purified zcurrent x! ; : ?set-wcurrent ( -- ) ( If the second complex argument flag W? is true, select a w at random and store it in WCURRENT. This version does not filter or purify w, in accord with Cody. ) w? IF deltay f@ ren f* corner1 x@ f+ (f: u=deltay*ren+AX) deltay f@ ren f* corner1 y@ f+ (f: u v=deltay*ren+AY) wcurrent z! THEN ; : set-deltas ( -- ) ( Set DELTAX and DELTAY from CORNER1, CORNER2, and #ARGS. DELTAX is plus or minus the width of the box divided by the number of points, and DELTAY is plus or minus the height of the box. ) corner2 x@ corner1 x@ f- #args s>f f/ deltax f! corner2 y@ corner1 y@ f- deltay f! ; : get-zstats ( -- ) ( Assume CORNER1 and CORNER2 have been set. Evaluate and compare FUNCTION and GAUGE at #ARGS random points in the box, and accumulate error information in the VEC-ERR, REAL-ERR, and IMAG-ERR structures. The real and imaginary parts of the points are computed differently: x is stepped by DELTAX and increased by a random variation in the next DELTAX interval, while y is a random value in the interval corresponding to the vertical sides of the box. I think Cody does it this way for filtering on y only. ) reset-error set-deltas corner1 x@ (f: XL=AX) #args 0 DO fdup corner1 y@ (f: XL XL AY) set-zcurrent ?set-wcurrent function gauge (f: XL) handle-vector-error fcurrent gcurrent real-err handle-part-error fcurrent float+ gcurrent float+ imag-err handle-part-error deltax f@ f+ (f: XL=XL+deltax) LOOP (f: XL) fdrop ; : .#fdigits ( -- ) ." There are " #digits . ." base " ibeta . ." significant digits in a floating point number." ; : .part-#same ( err.struc part.s -- ) rot ( err.struc) >r cr ( part.s) type ." part of function: larger " r@ >#f.gt.g @ 4 .r ." times" cr ." same " r@ part-#same 4 .r ." times" cr ." smaller " r> >#f.lt.g @ 4 .r ." times" cr ; : prec-f. ( prec f: r -- ) precision swap set-precision f. set-precision ; : prec-fs. ( prec f: r -- ) precision swap set-precision fs. set-precision ; : .error ( exp.prec err.prec f: err -- ) ( err.prec) fdup prec-fs. ." = " ibeta 2 .r ." ^" ( exp.prec) ln_beta prec-f. ; : .digit-loss ( prec f: err -- ) ." Estimated loss of base " ibeta . ." significant digits: " ulp prec-f. ; : ?.wmax ( err.struc -- ) w? IF cr ." w = " ( err.struc) >wmax z@ zs. ELSE drop THEN ; : .max-err ( err.struc -- ) ( err.struc) >r precision 5 set-precision cr ." Maximum relative error: " r@ mre (f: w) fdup 2 3 .error cr ." at z = " r@ >zmax z@ zs. r@ ?.wmax cr ." f = " r> >fmax z@ zs. cr (f: w) 2 .digit-loss cr set-precision ; : .rms-err ( err.struc -- ) cr ." Root mean square relative error: " ( err.struc) rms-err (f: w) fdup 2 3 .error cr (f: w) 2 .digit-loss cr ; : .zstats ( -- ) ( Assume the VEC-ERR, REAL-ERR, and IMAG-ERR structures, and #ARGS. Display a statistical report. ) precision 3 set-precision cr ." The function and the gauge were compared for " #args . ." random arguments" cr ." in the box with opposite corners: z = " corner1 z@ zs. cr ." z = " corner2 z@ zs. cr cr ." Comparison: the same " vec-err vec-#same 4 .r ." times" cr ." different " vec-err >#f.ne.g @ 4 .r ." times" cr vec-err .max-err vec-err .rms-err real-err s" Real" .part-#same real-err .max-err real-err .rms-err imag-err s" Imag" .part-#same imag-err .max-err imag-err .rms-err set-precision ; : .short-stats-header ( -- ) ." N= MRE ULP RMS ULP" ; : .short-errs ( err.struc -- ) >r r@ mre 2 prec-fs. space r@ mre-ulp 2 prec-f. space r@ rms-err 2 prec-fs. space r> rms-ulp 2 prec-f. ; : .short-vec-stats ( -- ) ." vector " vec-err vec-#same 4 .r 2 spaces vec-err .short-errs ; : .short-part-stats ( part.err.struc part.name.s -- ) type dup part-#same 4 .r 2 spaces .short-errs ; : .short-zstats ( -- ) ( Assume the VEC-ERR, REAL-ERR, and IMAG-ERR structures, and #ARGS. Display an abbreviated statistical report. ) cr .short-stats-header cr .short-vec-stats cr real-err s" real " .short-part-stats cr imag-err s" imaginary " .short-part-stats ; [ELSE] .( ***Separate floating point stack not available.) [THEN] [ELSE] .( ***Floating point not available.) [THEN]