( Zelefunt: Test complex absolute value File: zabs.fs Version: 0.9.3 Original Author: W. J. Cody [November 7, 1990] Ported by: David N. Williams License: ACM noncommercial use Last revision: January 19, 2005 This is a port of W. J. Cody's CABS program in the celefunt package from Fortran 77 to ANS Forth. As ACM Algorithm 714 [TOMS], that package is presumed to be under the ACM license for noncommercial use: http://www.acm.org/pubs/copyright_policy/softwareCRnotice.html There is an ANS Forth environmental dependence on lower case. Quote from the original package: Program to test CABS Accuracy tests are based on the identities cabs[[3x,4x]] = 5x and cabs[[5x,12x]] = 13x 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 six 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] Latest revision - November 7, 1990 Author - W. J. Cody Argonne National Laboratory Porting notes: Version 0.9.3 19Jan05 * Changed large argument in the last calculation from xmax/4 + i xmax to [xmax/4]*3 + i xmax to agree with celefunt. 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, including alignment of the version number with the rest of zelefunt. Version 0.9.0 10Jan03 * Start. 10Feb03 * 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 t for Cody's x, and we use x and y for the real and imaginary parts of z. ) decimal s" FLOATING-EXT" environment? [IF] ( flag) drop s" FLOATING-STACK" environment? [IF] ( maxdepth) drop loadm complex s" machar.fs" included xmin xmax beta ibeta #digits -MACHAR constant #digits constant ibeta fconstant beta fconstant xmax fconstant xmin s" ren.fs" included \ loadm complex \ pfe \ s" complex.fs" included s" complex-kahan.fs" included \ s" zabs.seq" included \ FSL s" [UNDEFINED]" pad c! pad char+ pad c@ move pad find nip 0= [IF] : [UNDEFINED] ( "name" -- flag ) bl word find nip 0= ; immediate [THEN] [UNDEFINED] (f: [IF] : (f: ( -- ) postpone ( ; immediate [THEN] [UNDEFINED] f>s [IF] : f>s (f: r -- s: n ) f>d drop ; [THEN] [UNDEFINED] s>f [IF] : s>f ( n -- f: r ) 0 d>f ; [THEN] [UNDEFINED] f0<> [IF] : f0<> (f: x y -- s: flag ) f0= 0= ; [THEN] [UNDEFINED] f2dup [IF] : f2dup (f: x y -- x y x y ) fover fover ; [THEN] [UNDEFINED] z. [IF] : z. (f: x y -- ) fswap f. ." + i " f. ; [THEN] [UNDEFINED] zs. [IF] : zs. (f: x y -- ) fswap fs. ." + i " fs. ; [THEN] 0 [IF] \ naive calcution : |z| ( x y -- sqrt|x^2+y^2| ) fdup f* fswap (f: y^2 x) fdup f* (f: y^2 x^2) f+ fsqrt ; cr ." Using naive |z|." cr [THEN] beta fln fconstant albeta #digits s>f fconstant ait 1e fconstant a 20e fconstant b 2000 constant #args #args s>f fconstant f#args b a f- f#args f/ fconstant del -999e fconstant all9 fvariable tl \ origin for random t variation fvariable xfac \ x=xfac*t, e.g., 3e*t fvariable yfac \ y=yfac*t, e.g., 4e*t fvariable absfac \ |z|=absfac*t, e.g., 5e*t fvariable max-err \ maximum magnitude of relative error fvariable tmax \ t value at maximum relative error fvariable sum-sqrs \ sum of squares of relative errors : x (f: t -- xfac*t ) xfac f@ f* ; : y (f: t -- yfac*t ) yfac f@ f* ; : targ (f: t -- absfac*t) absfac f@ f* fabs ; : calc (f: t -- |xfac*t+iyfac*t| ) fdup x fswap y |z| ; : rel-err (f: t -- [calc-targ]/calc ) fdup targ fswap calc fdup (f: targ calc calc) frot f- (f: calc calc-targ) fswap f/ ; : get-stats ( -- #above #below ) ( Assume XFAC, YFAC, and ABSFAC to be set, with XFAC^2 + YFAC^2 = ABSFAC^2. Step TL over the range of t values, adding a random variation to get the t value actually used. Compute |XFAC*t + i YFAC*t| and compare to the expected value ABSFAC*t. Set MAX-ERR, TMAX, and SUM-SQS; and leave the number of results above and the number below the expected value. ) 0e tl f! 0e max-err f! 0e tmax f! 0e sum-sqrs f! 0 0 ( #above=0 #below=0) #args 0 DO ren del f* tl f@ f+ (f: t=ren*del+tl) fdup 8e f* fswap fover f+ f- (f: t=[t+8t]-8t) fdup rel-err (f: t err=[calc-targ]/calc) fdup f0< IF ( #below) 1+ ELSE fdup f0<> IF swap ( #above) 1+ swap THEN THEN fdup fdup f* (f: t err err^2) sum-sqrs dup f@ f+ f! (f: t err) fabs max-err f@ fover (f: t |err| max.err |err|) f< IF max-err f! tmax f! ELSE fdrop fdrop THEN tl f@ del f+ tl f! LOOP ; 0 [IF] \ for debugging : .vars ( -- ) cr ." tl " tl f@ f. cr ." max-err " max-err f@ fs. cr ." tmax " tmax f@ fs. cr ." sum-sqrs " sum-sqrs f@ fs. cr ." calc at tmax " tmax f@ calc fs. cr ." targ at tmax " tmax f@ targ fs. ; [THEN] : log_2 (f: f -- log_2[f] ) fdup f0= IF fdrop all9 ELSE fln albeta f/ THEN ; : .error (f: err log_2[err] -- log_2[err] ) fswap fs. ." = " ibeta 2 .r ." **" fdup f. ; : .digit-loss (f: log_2[err] -- ) ait f+ 0e (f: #digs+log_2[err] 0) fmax ." Estimated loss of base " ibeta . ." significant digits: " f. ; : .stats ( #above #below -- ) ( Assume XFAC, YFAC, ABSFAC, MAX-ERR, TMAX, and SUM-SQRS. Display a statistical report. ) cr ." Statistics for |" xfac f@ f>s . ." t + i " yfac f@ f>s . ." t| vs. " absfac f@ f>s . ." t" cr cr #args . ." random t values were drawn from the interval (" a f. ." , " b f. ." )." 2dup + #args swap - rot ( #below #same #above) cr cr ." The calculated |z| was:" cr ." larger " 4 .r ." times" cr ." same " 4 .r ." times" cr ." smaller " 4 .r ." times" cr cr ." The maximum relative error of " max-err f@ fdup log_2 (f: err=max w=log_2[err]) .error (f: w) cr ." occurred for t = " tmax f@ f. ." ." cr (f: w) .digit-loss cr cr ." The root mean square relative error was " sum-sqrs f@ f#args f/ fsqrt fdup log_2 (f: err=rms w=log_2[err]) .error ." ." (f: w) cr .digit-loss cr ; cr .( RANDOM ARGUMENT TESTS) cr cr .( There are ) #digits . .( base ) ibeta . .( significant digits in a floating-point number.) cr 3e xfac f! 4e yfac f! 5e absfac f! get-stats .stats 5e xfac f! 12e yfac f! 13e absfac f! get-stats .stats cr .( SPECIAL ARGUMENT TESTS) cr cr .( Testing |3*xmin + i 4*xmin| = 5*xmin,) cr .( which should *not* trigger an error message:) cr xmin 5e f* xmin 3e f* xmin 4e f* |z| (f: 5*xmin |3*xmin+i*4*xmin|) cr .( |3*xmin + i 4*xmin| = ) fdup fs. .( = |z|) cr .( 5*xmin = ) fover fs. cr .( |z| - 5*xmin = ) fswap f- fs. cr \ dnw: next test not in celefunt cr .( Testing |3*[xmax/5] + i 4*[xmax/5]| = xmax,) cr .( which should *not* trigger an error message:) cr xmax 5e f/ fdup 3e f* fswap 4e f* |z| cr .( |3*[xmax/5] + i 4*[xmax/5]| = ) fdup fs. .( = |z|) cr .( xmax = ) xmax fdup fs. cr .( |z| - xmax = ) f- fs. cr cr .( FLOATING POINT EXCEPTION TESTS) cr xmax 16e f/ fdup 5e f* fswap 12e f* f2dup cr .( Calling |z| with the argument ) cr 2 spaces zs. .( ,) cr .( which should give 13*[xmax/16] *without* an error message:) cr cr .( |z| = ) |z| fdup fs. cr .( 13*[xmax/16] = ) xmax 16e f/ 13e f* fdup fs. cr .( |z| - 13*[xmax/16] = ) f- fs. cr xmax 4e f/ 3e f* xmax cr .( Calling |z| with the argument ) cr 2 spaces .( [xmax/4]*3 + i xmax = ) f2dup zs. .( ,) cr .( which *should* trigger an error message or overflow result:) cr cr .( Value returned: |z| = ) |z| fs. cr cr .( END OF TESTS) cr [ELSE] .( ***Separate floating point stack not available.) [THEN] [ELSE] .( ***Floating point not available.) [THEN]