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