C This file: http://ftp.aset.psu.edu/pub/ger/fortran/hdk/fuzzloop.for
C
C====This is a program to illustrate how DO loops with REAL
C parameters can compromise the implied trip count. The reason is
C significance loss in the DO index introduced by summing errors
C in representation at each trip (exaserbated by finite precision
C representation of numbers in the machine base).
C
C A better loop specification (see loop 3) and the use of the well
C established "fuzzy" ROUND function are shown to illustrate
C circumvention of this numerical problem.
C MAKE CERTAIN THAT ALL COMPILER OPTIMIZATION IS OFF.
C
C H. D. KNOBLE, PENN STATE UNIVERSITY - AUGUST - 1982.
C hdkLESS at SPAM psu dot edu
C================================================================
CHARACTER*5 FLAG
CHARACTER*16 MACH
LOGICAL TEQ
REAL XSTART,XINC,X
INTEGER I,N,NTRIPS,TRIPS,DUMMY
DOUBLE PRECISION ROUND,CT,DX,DSTART,DEND,DINC,EPS,D1MACH,U,V
DATA MACH /'PC (WINDOWS)'/
TEQ(U,V)=DABS(U-V).LE.MAX(DABS(U),DABS(V))*CT
C--------EPS, the "Machine Precision" is defined as:
C the smallest single precision REAL number such that
C 1+EPS>1 AND 1-EPS<1.
EPS=D1MACH(DUMMY)
C--------CT=3*EPS essentially means: "accept a floating-point
C representation if it is either exact or one floating-point
C representation to either side (within a bit) of the number
C in question.
CT=3.*EPS
C
C--------Hypothetical problem: plot a smooth curve of a function, F(X),
C on a plotter with increment of 0.1 inches (rougher than typical
C line plotter accuracy). Define the curve for X over the closed
C interval, 0 .LE. X .LE. N.
C
N=10000
XSTART=0.
XEND=N
DEND=N
XINC=.1
DINC=.1D0
TRIPS=ROUND(REAL(N)/DINC,CT)+1
WRITE(*,100) MACH,EPS,TRIPS,XSTART,XEND,XINC
100 FORMAT('Computer: ',A16,' EPS=',E22.15,' True Trip Count=',I6//
* 'Loop is: DO X=XSTART,XEND,XINC where:'/
* ' XSTART=',F8.2,', XEND=',F8.2,', XINC=',F8.2//
* ' Number Final',
* ' Last Value'/
* ' DO Loop Algorithm of Trips Index Value',
* ' of (XEND-X)'/' ',77('-')/)
C
C--------Loop 1 is equivalent to the Fortran 77 specification of:
C DO 1 X=XSTART,XEND,XINC
C
C Significance loss in representation of the increment and starting
C value accumulates; that is the index progressively includes larger
C amounts of error for each successive trip. This algorithm requires
C essentially NTRIPS additions. But note that the probability for
C NTRIPS being calculated wrong is quite high.
C
NTRIPS=MAX(INT((XEND-XSTART+XINC)/XINC),0)
X=XSTART-XINC
DO 1 I=1,NTRIPS
X=X+XINC
C Evaluate F(X) (and hypothetically plot scaled point (X,F(X))).
C CALL PLOT(X,F(X))
1 CONTINUE
FLAG=' '
IF(NTRIPS.EQ.TRIPS .AND. TEQ(DBLE(X),DBLE(XEND))) FLAG=' <---'
WRITE(*,10) NTRIPS,X,DBLE(XEND)-DBLE(X),FLAG
10 FORMAT(' FORTRAN 77 ANSI X3.9 DO loop (SP):',
* T40,I7,T50,F9.2,1PE14.3,A5)
C--------Loop 2 shows the current implementation of the REAL
C parameter DO loop. It does a compare like:
C IF(X.GT XEND) THEN "quit this loop."
C Like Loop 1 above, this will quite often take the wrong number
C of trips. The intermediate values of the loop index X are as
C bad as loop 1, and the final value of the index is more is
C accurate, but only by chance alone in this case.
C
NTRIPS=0
DO 2 X=XSTART,XEND,XINC
NTRIPS = NTRIPS + 1
C Evaluate F(X) (and hypothetically plot scaled point (X,F(X))).
C CALL PLOT(X,F(X))
2 CONTINUE
FLAG=' '
IF(NTRIPS.EQ.TRIPS .AND. TEQ(DBLE(X),DBLE(XEND))) FLAG=' <---'
WRITE(*,20) NTRIPS,X,DBLE(XEND)-DBLE(X),FLAG
20 FORMAT(' PC WINDOWS Fortran DO loop(SP):',
* T40,I7,T50,F9.2,1PE14.3,A5)
C
C--------Loop 3 is an optimal calculation of the number of trips and
C will always effect an evaluation of the loop index X that is as
C accurate or more accurate than loops 1 or 2 above. Significance
C loss is not accumulative, but constant (error of representation)
C for each trip. The error in the final index value is optimal and
C the number of trips is optimally accurate. This algorithm requires
C NTRIPS conversions from INTEGER to REAL and one call to the
C Tolerant rounding function, ROUND. This "fuzzy" ROUND function
C makes the probability of calculating the correct number of trips
C higher. The difference in CPU time between loop 3 and loops 1
C or 2 is negligible.
C
NTRIPS=MAX(ROUND((XEND-DBLE(XSTART)+XINC)/XINC, CT),0.D0)
DO 3 I=1,NTRIPS
X=XSTART+(I-1)*XINC
C Evaluate F(X) (and hypothetically plot scaled point (X,F(X))).
C CALL PLOT(X,F(X))
3 CONTINUE
FLAG=' '
IF(NTRIPS.EQ.TRIPS .AND. TEQ(DBLE(X),DBLE(XEND))) FLAG=' <---'
WRITE(*,30) NTRIPS,X,DBLE(XEND)-DBLE(X),FLAG
30 FORMAT(/' Optimally accurate DO loop:',
* T40,I7,T50,F9.2,1PE14.3,A5/)
C
C---------Loop 4 is a Double Precision implementation of Loop 2. This
C results in a correct value for the number of trips (by chance) but
C still a less accurate final index value than loop 3 here.
C
DSTART=0.D0
DINC=.1D0
NTRIPS=0
DO 4 DX=DSTART,DEND,DINC
NTRIPS = NTRIPS + 1
C Evaluate F(X) (and hypothetically plot scaled point (X,F(X))).
C CALL PLOT(X,F(X))
4 CONTINUE
FLAG=' '
IF(NTRIPS.EQ.TRIPS .AND. TEQ(DX,DEND)) FLAG=' <---'
WRITE(*,40) NTRIPS,DX,DEND-DX,FLAG
40 FORMAT(' PC WINDOWS Fortran DO loop (DP):',
* T40,I7,T50,F9.2,1PE14.3,A5)
C
C---------Loop 5 is a Double Precision implementation of Loop 1. It
C is better, but like Loop 4 will break down "sooner" than Loop 3.
DSTART = 0.D0
DINC=.1D0
NTRIPS=MAX(INT((DEND-DSTART+DINC)/DINC),0)
DX=DSTART-DINC
DO 5 I=1,NTRIPS
DX=DX+DINC
C Evaluate F(X) (and hypothetically plot scaled point (X,F(X))).
C CALL PLOT(X,F(X))
5 CONTINUE
FLAG=' '
IF(NTRIPS.EQ.TRIPS .AND. TEQ(DX,DEND)) FLAG=' <---'
WRITE(*,50) NTRIPS,DX,DEND-DX,FLAG
50 FORMAT(' Fortran 77 ANSI X3.9 DO loop (DP):',
* T40,I7,T50,F9.2,1PE14.3,A5)
C
WRITE(*,200)
200 FORMAT(' ',77('-')/ ' The symbol "<---" indicates that ',
* 'these are the mathematically correct values.')
WRITE(*,*) '** Loop Example Ended **'
STOP
END
DOUBLE PRECISION FUNCTION D1MACH (IDUM)
INTEGER IDUM
C===This routine computes the unit round-off of the machine in Double
C Precision. This is defined as the smallest positive machine real
C real number, EPS, such that (1+EPS .GT. 1) & (1-EPS .LT. 1).
C This particular computation of EPS is the work of Alan C. Hindmarsh.
C
DOUBLE PRECISION EPS, COMP
EPS = 1.0D0
10 EPS = EPS*0.5D0
C---Note from Richard Hanson:
C Make sure the following statement isn't optimized out of the loop.
C could do this by making it a (separate) function and call here.
COMP = 1.0D0 + EPS
IF (COMP .NE. 1.0D0) GO TO 10
D1MACH = EPS*2.0D0
RETURN
END
DOUBLE PRECISION FUNCTION ROUND(X,CT)
C===Tolerant rounding function
C See Knuth, Art Of Computer Programming, Vol. 1, Problem 1.2.4-5.
DOUBLE PRECISION TFLOOR,X,CT
ROUND=TFLOOR(X+0.5D0,CT)
RETURN
END
DOUBLE PRECISION FUNCTION TFLOOR(X,CT)
C===Tolerant FLOOR function.
C
C X - is given as a Double Precision argument to be operated on.
C It is assumed that X is represented with M mantissa bits.
C CT - is given as a Comparison Tolerance such that
C 0.LT.CT.LE.3-SQRT(5)/2. If the relative difference between
C X and A whole number is less than CT, then TFLOOR is
C returned as this whole number. By treating the
C floating-point numbers as a finite ordered set note that
C the heuristic EPS=2.**(-(M-1)) and CT=3*EPS causes
C arguments of TFLOOR/TCEIL to be treated as whole numbers
C if they are exactly whole numbers or are immediately
C adjacent to whole number representations. Since EPS, the
C "distance" between floating-point numbers on the unit
C interval, and M, the number of bits in X'S mantissa, exist
C on every floating-point computer, TFLOOR/TCEIL are
C consistently definable on every floating-point computer.
C
C For more information see the following references:
C (1) P. E. Hagerty, "More On Fuzzy Floor And Ceiling," APL QUOTE
C QUAD 8(4):20-24, June 1978. Note that TFLOOR=FL5.
C (2) L. M. Breed, "Definitions For Fuzzy Floor And Ceiling", APL
C QUOTE QUAD 8(3):16-23, March 1978. This paper cites FL1 through
C FL5, the history of five years of evolutionary development of
C FL5 - the seven lines of code below - by open collaboration
C and corroboration of the mathematical-computing community.
C
C Penn State University Center for Academic Computing
C H. D. Knoble - August, 1978.
C=====================================================================
DOUBLE PRECISION X,Q,RMAX,EPS5,CT,FLOOR,DINT
C---------FLOOR(X) is the largest integer algebraically less than
C or equal to X; that is, the unfuzzy FLOOR function.
DINT(X)=X-DMOD(X,1.D0)
FLOOR(X)=DINT(X)-DMOD(2.D0+DSIGN(1.D0,X),3.D0)
C---------Hagerty's FL5 function follows...
Q=1.D0
IF(X.LT.0)Q=1.D0-CT
RMAX=Q/(2.D0-CT)
EPS5=CT/Q
TFLOOR=FLOOR(X+MAX(CT,MIN(RMAX,EPS5*DABS(1.D0+FLOOR(X)))))
IF(X.LE.0 .OR. (TFLOOR-X).LT.RMAX) RETURN
TFLOOR=TFLOOR-1.D0
RETURN
END
DOUBLE PRECISION FUNCTION TCEIL(X,CT)
C===Tolerant ceiling function.
C See TFLOOR.
DOUBLE PRECISION X,CT,TFLOOR
TCEIL= -TFLOOR(-X,CT)
RETURN
END