c SRTEST.FTN: A Fortran program to test any SRT divider's logic c for quotient-digit selection, usually in a PLA, c regardless of most internal details of the divider. c This version runs under MS-DOS on Intel-compatible floating-point. c My compiler provides no access to the DOS command tail, so input c parameters come from four integers in four-line file SRTEST.DTA : c 0 < n < 17 , 0 < Lx < 17 , ix < 17 , and odd dvd > 0 . c c ... Parameter n says how many leading bits of divisors 1.xxxxx... c : enter into SRT quotient-selection: -> n <- c : Over-estimating n will slow SRTEST exponentially. c : ======> Under-estimating n will Invalidate SRTEST . <======== c : ~~~~~~~~~~ c : Parameter Lx = Log2(Radix) > 0 , so Radix = 2**Lx > 1 . c : Over-estimating Lx slows SRTEST proportionately. c : Under-estimating Lx may speed SRTEST up a bit or slow c : it drastically. If unknown, try Lx = 1 first. c : c : Parameter ix increments and decrements dividends by 2**ix . c : It is still subject to experimental investigation. If on c : input ix < 0 it defaults to ix = Lx - 1 ; else try 0. c : c : For Intel's Pentium, n = 5 , Lx = 2 , ix = 1 . c : c : Parameter dvd should start at 1 but may be set to any bigger c :.. odd integer to resume SRTESTing aborted by Ctrl-Break. c c ... Every variable has a declared type: IMPLICIT NONE CHARACTER*48 aargh DATA aargh / ' TO ABORT THIS TEST, PRESS [Ctrl][Break] .' / c c ... The next few declarations are determined by the format of the c : machine's floating-point arithmetic. Intel's is five 2-byte c : words wide, with sign bit and biased exponent in the first word, c : and with an explicit leading 1 in its 2*m = 64 sig. bits. c : Later the program will set tm = 2.0**m and tmn = tm/2.0**n in c : a way invulnerable to the compiler's whimsies. The EQUIVALENCEd c :.. variables iq(...), qe, qh, dq, sq are used for bit-twiddling: INTEGER*2 m, iq(5) PARAMETER ( m = 32 ) REAL*10 qe, qh, tm, tmn REAL*8 dq REAL*4 sq EQUIVALENCE ( dq, iq, qe, qh, sq ) c --- <- qh's exp.-><--- qh's sig. bits ---> Little-Endian. c --- <-- sq --> To be set to +0.0 . c --- <------ dq ------> To be set to -0.0 . c --- <- qe's exp.-><1000 ... ... 000> To expose exponent. c --- <- iq(5) ->< iq(4),iq(3),iq(2),iq(1) > To adjust exponent. c c ... Glossary: c : count = count of divisions tested so far. c : dvd = unperturbed small odd integer dividend; 0 < dvd < dvdend. c : dvdend = 2.0**(m - 1 - ix) c : dvdn0 = normalized dvd , an m-bit integer. c : dvdn = dvdn0 + idn , an actual dividend. c : dvdx = dvdn at worst failure found so far. c : dvr0 runs through i2n1 + 1 singular divisors. c : dvr = dvr0 + s(k) is a perturbed singular divisor. c : dvrx = dvr at worst failure found so far. c : i2n1 = 2**(n-1) c : icw = integer to set Intel's Floating-Point Control Word. c : iem1 = biased exponent of +2.0**(m-1) c : k = index that runs from k0 to kl ; c : k0 = -Lx typically; kl = Lx typically. c : kount = count of failed divisions so far. c : kountx = kount at worst failure found so far. c : L counts down to print every L0-th value dvd . c : quo = quotient under test = qh + qt when split in halves. c : quox = quo at worst failure found so far. c : r = remainder dvdn - quo*dvr exactly. c : rx = r at worst failure found so far. c : s(...) = divisor perturbations {..., -2, -1, 0, 1, 2, 4, ...} . c : sd = 2.0**ix c : sdn = -sd, 0, sd to perturb dividends dvdn = dvdn0 + sdn . c : tm = 2.0**m ; tmn = 2.0**(n-m) . c : zer0m = -0.0 D0 . c :.. c INTEGER*2 icw, iem1, i2n1, ix, j, k, k0, kl INTEGER*2 L, L0, Lx, m2, n INTEGER*4 kount, kountx, i2 PARAMETER ( i2 = 2, L0 = 500, m2 = m*2 ) REAL*10 count REAL*10 dvd, dvdend, dvdn, dvdn0, dvdx, dvr, dvr0, dvrx REAL*10 quo, quox, qt, r, rx REAL*8 zer0m REAL*4 s(-16 : 16), sd, sdn c c < < < < < < < < < < < for DEBUGGING < < < < < < < < < < < < < < < < LOGICAL*2 L1, L2, L3 DATA L1, L2, L3 / .FALSE., .FALSE., .FALSE. / c > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > c c ... Initializations: DATA count, kount, kountx, rx / 0.0, 0, 0, 0.0 / DATA dvdx, L, sq / 1.0, 1, 1.0 / c OPEN( *, CARRIAGE='FORTRAN' ) OPEN( UNIT=3, FILE='SRTFAILU.RES' ) OPEN( UNIT=4, FILE='SRTEST.DTA' ) READ (4,*) n IF ( (n .LE. 0) .OR. (n .GE. 17) ) THEN PRINT *, ' SRTEST requires 0 < n < 17 , not', n STOP END IF READ (4,*) Lx IF ( (Lx .LE. 0) .OR. (Lx .GE. 17) ) THEN PRINT *, ' SRTEST needs 0 < Lx < 17 , not', Lx STOP END IF READ (4,*) ix IF ( ix .LT. 0 ) ix = Lx-1 IF ( ix .GE. 17 ) THEN PRINT *, ' SRTEST requires ix < 17 , not', ix STOP END IF sd = i2**ix READ (4,*) dvd dvd = RINT( ABS(dvd) ) IF ( RINT(0.5*dvd)*2.0 .EQ. dvd ) THEN PRINT *, ' SRTEST needs an odd integer dvd , not', dvd STOP END IF CLOSE( UNIT=4 ) WRITE(3,*) ' SRTEST: with n =', n, ', Lx =', Lx WRITE(3,*) ' and with ix =', ix, ' , has these failures:' WRITE(3,*) ' (#, Dividend, Divisor, Quotient, Remainder, Count)' c ... This CLOSE and subsequent OPEN etc. flush the disk buffer: CLOSE( UNIT=3 ) PRINT *, ' Program running now is SRTEST, a test of Intel''s' PRINT *, ' Rounded Floating-Point Division upon integers chosen' PRINT *, ' to explore ANY SRT-divider''s quotient-selection' PRINT 1, n 1 FORMAT( ' logic assuming that it samples only the first ', I5 ) PRINT *, ' sig. bits of a divisor, and assuming that all other' PRINT *, ' floating-point operations and conversions are O. K.' PRINT *, ' A guess at the divider''s Radix =', i2**Lx PRINT *, ' The dividend increment sd =', sd PRINT *, aargh PRINT * PRINT * c ... Set Floating-Point Control-Word to default round-to-nearest 64 : icw = #033Dh c < < < < < DEBUGGING TEST for FAILURE by rounding to 53 < < < < < < IF (L1) icw = #023Dh c > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > CALL LDCW87( icw ) c c ... Initialize array s(...) of divisor perturbations: c s(...) = {-Radix/2, ..., -2, -1, 0, 1, 2, 4, ..., Radix/2 }. s(0) = 0.0 DO 100, k = 1, 16 s(k) = sq s(-k) = -sq sq = sq+sq 100 END DO i2n1 = s(n) c c ... Create tm = 2.0**m to help produce divisors with at most m c nonzero leading bits and to split quotients with 2*m explicit c sig. bits into two m-bit numbers; also tmn = 2.0**(m-n) : qh = dvdx iq(5) = iq(5) + m tm = qh iem1 = iq(5) - 1 iq(5) = iq(5) - n tmn = qh iq(5) = iem1 - ix dvdend = qh c ... Create zer0m = -0.0 D0 despite the compiler : zer0m = dq c < < < < < < < < < < < for DEBUGGING < < < < < < < < < < < < < < < < IF (L2) THEN PRINT *, ' i2n1 = ', i2n1 PRINT *, ' tm = ', tm PRINT *, ' tmn = ', tmn PRINT *, ' zer0m = ', zer0m PRINT *, ' dvd = ', dvd PRINT *, ' dvdend = ', dvdend dvdend = 2 PRINT *, ' dvdend = ', dvdend PRINT *, s PRINT * END IF c > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > OPEN( UNIT=3, FILE='SRTFAILU.RES' ) REWIND(3) DO 200, k = 1, 3 READ(3,*) 200 END DO c ... Ends DO-loop with UNIT 3 positioned just before file's end. c DO 600 WHILE ( dvd .LT. dvdend ) c ... Dividend dvd = 1.0, 3.0, 5.0, 7.0, 9.0, ..., dvdend - 1 : L = L-1 IF ( L .EQ. 0 ) THEN PRINT 2, dvd c ... This line overwrites itself continually: 2 FORMAT( '+ Currently testing at dvd >', F20.1 ) L = L0 END IF c ... Normalize dividend dvd to an m-bit integer dvdn0 : qh = dvd iq(5) = iem1 dvdn0 = qh k0 = 0 kl = Lx DO 500, j = i2n1, 0, -1 c ... Run dvr0 through singular divisors: dvr0 = tm - j*tmn IF ( j .EQ. 0 ) kl = -1 DO 400, sdn = -sd, sd, sd c For each of three perturbed dividends dvdn = ... : dvdn = dvdn0 + sdn DO 300, k = k0, kl c ... Run dvr through (typically) 2*Lx+1 perturbations: dvr = dvr0 + s(k) c ... #### Here is the division operation to be tested: #### quo = dvdn/dvr count = count + 1.0 qh = quo c ... Split quo into half-precision qh + qt : sq = 0.0 qt = quo - qh c ... Remainder r = dvdn - dvr*quo EXACTLY : r = ( dvdn - dvr*qh ) - dvr*qt c < < < < < < < < < < for DEBUGGING < < < < < < < < < < < < < < < < < IF (L3) THEN PRINT *, ' count = ', count, ' dvd = ', dvd PRINT 3, 'dvdn = ', dvdn, dvdn 3 FORMAT( 1X, A7, G25.15, '= ', Z20 ) PRINT 3, 'dvr = ', dvr, dvr PRINT *, ' quo = ', quo, ' r = ', r END IF c > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > c ... Expose and decrease exponent of quo by 2*m dq = zer0m iq(5) = iq(5) - m2 c ... Test remainder against rounding error bound: IF ( ABS(r) .GE. qe*dvr ) THEN c ... Test has failed: kount = kount + 1 PRINT * PRINT *, ' Division Test has FAILED for the' PRINT *, ' #', kount, 'th time after', count PRINT *, ' divisions at dvd =', dvd, ', and ...' PRINT 4, 'Dividend ', dvdn, dvdn 4 FORMAT( 5X, A9, ' = ', F20.1, ' = ', Z20 ) PRINT 4, 'Divisor ', dvr, dvr PRINT 5, 'Remainder', r 5 FORMAT( 5X, A9, ' = ', 1PE26.17 ) PRINT 5, 'Quotient ', quo WRITE(3,6) kount, dvdn, dvr, quo, r, count 6 FORMAT( 1X, I8, 2F13.0, F22.17, 1PE22.12, 0PF15.0 ) IF ( ABS(r*dvdx) .GT. ABS(rx*dvdn) ) THEN rx = r dvdx = dvdn dvrx = dvr quox = quo kountx = kount END IF PRINT *, ' Relatively worst remainder so far is #', kountx PRINT 5, 'Remainder', rx PRINT 4, 'Dividend ', dvdx, dvdx PRINT 4, 'Divisor ', dvrx, dvrx PRINT 5, 'Quotient ', quox PRINT *, aargh PRINT * PRINT * END IF 300 END DO c ... Ends DO-loop on k from k0 to kl . 400 END DO c ... Ends DO-loop on sdn = -sd to +sd . k0 = -Lx 500 END DO c ... Ends DO-loop on j = i2n1 to 0 . dvd = dvd + 2.0 600 END DO c ... Ends WHILE-loop on dvd < dvdend . PRINT *, count, ' divisions have been tested;' IF ( kount .EQ. 0 ) THEN PRINT *, ' all PASSED the test.' ELSE PRINT *, kount, ' divisions FAILED the test;' END IF PRINT *, ' Copyright (C) W. Kahan, 1 March 1995' CLOSE( UNIT=3 ) END