( Zelefunt: Test complex sine and cosine File: zsincos.fs Version: 0.9.3 Original Author: W. J. Cody [October 31, 1990] Ported by: David N. Williams License: ACM noncommercial use Last revision: January 19, 2005 This is a port of W. J. Cody's complex sine and cosine 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 CSIN/CCOS Accuracy tests are based on the identities sin[z] = sin[z-w]*cos[w] + cos[z-w]*sin[w] and cos[z] = cos[z-w]*cos[w] - sin[z-w]*sin[w] 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 five 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; XMAX - the largest finite floating-point no. REN[K] - a function subprogram returning random real numbers uniformly distributed over [0,1] RESET, TABLAT, REPORT - programs to report results. Latest revision - October 31, 1990 Author - W. J. Cody Argonne National Laboratory Porting notes: Version 0.9.3 19Jan05 * Added summary output lines using .SHORT-ZSTATS from zfunstat.fs version 0.9.3. * Replaced -FROT by FROT FROT for portability. 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 * Start. 6Apr03 * 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 s" complex-kahan.fs" included s" zfunstat.fs" included \ includes machar.fs cr .( RANDOM ARGUMENT ACCURACY TESTS) cr cr .#fdigits cr 1e 16e f/ fdup zconstant del 6.2581348413276935585e-2 6.2418588008436587236e-2 zconstant sindel -2.5431314180235545803e-6 -3.9062493377261771826e-3 zconstant cosdel-1 3 set-precision cr .( Testing the function sin[z] against the gauge sin[w]*cos[del] +) cr .( cos[w]*sin[del], z = w + del, del = 1/16 + i*1/16.) cr \ Celefunt doesn't report the comparison of the imaginary part \ in this section, although it computes the function. The code \ below reports both parts as usual. :noname ( -- ) zcurrent z@ zsin fcurrent z! ; is function :noname ( -- ) zcurrent z@ del z- (f: w) zdup zcos zswap zsin (f: c=cos[w] s=sin[w]) zdup cosdel-1 z* z+ (f: c s+s*[cosdel-1]) zswap sindel z* z+ (f: s*cosdel+c*sindel) gcurrent z! ; is gauge ' noop is purified 1e 16e f/ 0e corner1 z! 1e 1e corner2 z! get-zstats .zstats cr .( Summary:) .short-zstats cr cr .( NEW BOX) cr 16e 16e corner1 z! 17e 17e corner2 z! get-zstats .zstats cr .( Summary:) .short-zstats cr cr .( Testing the function cos[z] against the gauge cos[w]*cos[del] -) cr .( sin[w]*sin[del], z = w + del, del = 1/16 + i*1/16.) cr :noname ( -- ) zcurrent z@ zcos fcurrent z! ; is function :noname ( -- ) zcurrent z@ del z- (f: w) zdup zsin zswap zcos (f: s=sin[w] c=cos[w]) zdup cosdel-1 z* z+ (f: s c+c*[cosdel-1]) zswap sindel z* z- (f: s*cosdel-s*sindel) gcurrent z! ; is gauge get-zstats .zstats cr .( Summary:) .short-zstats cr cr .( SPECIAL TESTS) cr 1e beta f/ #digits 2/ f^n fconstant k1 6 set-precision cr .( If ) 22e k1 f- (f: x2=x-k1) 22e k1 f+ (f: x2 x1=x+k1) 12.5e zsin (f: x2 z1=sin[x1+i*12.5]) frot 12.5e zsin (f: z1 z2=sin[x2+i*12.5]) z- k1 f2* z/f (f: [z1-z2]/[2*k1]) zs. .( is not almost -134163 + i 1188,) cr .( ZSIN has the wrong period.) cr cr .( Test of the identity SIN[-Z] = -SIN[Z]:) cr .( Z) 34 spaces .( SIN[Z] + SIN[-Z]) :noname ( -- ) 5 0 DO cr ren 10e f* ren 10e f* zdup zs. 5 spaces zdup zsin zswap znegate zsin z+ zs. LOOP cr ; execute cr .( Test the significance in the real or imaginary components) cr .( near the zeros of the real sin or cos:) cr -134167e 593e zconstant sin11a -3.2928849557976510994e-1 7.8989452897413630729e-1 zconstant sin11b -1187e -134163e zconstant sin22a -5.6815858848429427210e-1 -3.8738909081835398666e-1 zconstant sin22b : .part-digit-loss ( addr len f: err -- ) ( part.s) 2>r precision 2 set-precision (f: err) fdup f0<> IF fabs fln albeta f/ ait f+ (f: w=ln[|err|]/albeta+ait) ELSE fdrop 0e (f: w=0) THEN 0e fmax (f: w=max[w,0]) cr ." Estimated loss of base " ibeta . ." digits in the " 2r> type ." part: " fdup f. cr ." This is " 1e f< IF ." acceptable." ELSE ." too large." THEN set-precision ; cr .( SIN[ ) 11e 12.5e zdup zs. .( ] = ) zsin (f: s1) zdup zs. zdup sin11a z- sin11b z- imag (f: s1 ims2=Im[[s1-sin11a]-sin11b]) frot frot imag f/ (f: c=ims2/Im[s1]) s" imaginary" .part-digit-loss cr cr .( SIN[ ) 22e 12.5e zdup zs. .( ] = ) zsin (f: s1) zdup zs. zdup sin22a z- sin22b z- real (f: s1 res2=Re[[s1-sin22a]-sin22b]) frot frot real f/ (f: c=res2/Re[s1]) s" real" .part-digit-loss cr 3 set-precision cr .( Test of error returns:) cr cr .( This should not trigger an error message:) cr .( SIN[ ) 1e xmax fln 0.5e f- zdup zs. .( ] = ) zsin zs. cr \ Celefunt says this *should* trigger an error, but it shouldn't. cr .( This should not trigger an error message:) cr .( SIN[ ) beta #digits f^n 1e zdup zs. .( ] = ) zsin zs. cr cr .( This should trigger an error message:) cr .( SIN[ ) 1e beta #digits f^n zdup zs. .( ] = ) zsin zs. cr cr .( END OF TESTS) cr [ELSE] .( ***Separate floating point stack not available.) [THEN] [ELSE] .( ***Floating point not available.) [THEN]