C**********************************************************************
C This file: http://ftp.aset.psu.edu/pub/ger/fortran/hdk/eps.for
C  ROUTINE:   FUZZY FORTRAN OPERATORS
C  PURPOSE:   Illustrate Hindmarsh's computation of EPS, and APL
C             tolerant comparisons, tolerant CEIL/FLOOR, and Tolerant
C             ROUND functions - implemented in Fortran.
C  PLATFORM:  PC Windows Fortran, Compaq-Digital CVF 6.1a, AIX XLF90
C  TO RUN:    Windows: DF EPS.F90
C             AIX: XLF90 eps.f -o eps.exe -qfloat=nomaf
C  CALLS:     none
C  AUTHOR:    H. D. Knoble <hdkLESS at SPAM psu dot edu>
C  DATE:      22 September 1978
C  REVISIONS:
C**********************************************************************
C
      DOUBLE PRECISION EPS,EPS3, X,Y,Z, DUMACH,TFLOOR,TCEIL,EPSF90
      LOGICAL TEQ,TNE,TGT,TGE,TLT,TLE
C---Following are Fuzzy Comparison (arithmetic statement) Functions.
C
      TEQ(X,Y)=DABS(X-Y).LE.DMAX1(DABS(X),DABS(Y))*EPS3
      TNE(X,Y)=.NOT.TEQ(X,Y)
      TGT(X,Y)=(X-Y).GT.DMAX1(DABS(X),DABS(Y))*EPS3
      TLE(X,Y)=.NOT.TGT(X,Y)
      TLT(X,Y)=TLE(X,Y).AND.TNE(X,Y)
      TGE(X,Y)=TGT(X,Y).OR.TEQ(X,Y)
C
C---Compute EPS for this computer.  EPS is the smallest real number on
C   this architecture such that 1+EPS>1 and 1-EPS<1.
C   EPSILON(X) is a Fortran 90 built-in Intrinsic function. Dumach and
C   Epsilon are identically equal on IEEE architecture.
C
      EPS=DUMACH()
      EPSF90=EPSILON(X)
      IF(EPS.NE.EPSF90) THEN
        WRITE(*,2)'EPS=',EPS,' .NE. EPSF90=',EPSF90
2       FORMAT(A,Z16,A,Z16)
      ENDIF
C---Accept a representation if exact, or one bit on either side.
      EPS3=3.D0*EPS
      WRITE(*,1) EPS,EPS, EPS3,EPS3
1     FORMAT(' EPS=',D16.8,2X,Z16, ', EPS3=',D16.8,2X,Z16)
C---Illustrate Fuzzy Comparisons using EPS3. Any other magnitudes will
C   behave similarly.
      Z=1.D0
      I=49
        X=1.D0/DBLE(I)
        Y=X*DBLE(I)
        WRITE(*,*) 'X=1.D0/',I,', Y=X*',I,', Z=1.D0'
        WRITE(*,*) 'Y=',Y,' Z=',Z
        WRITE(*,3) X,Y,Z
3       FORMAT(' X=',Z16,' Y=',Z16,' Z=',Z16)
C---Floating-point Y is not identical (.EQ.) to floating-point Z.
        IF(Y.EQ.Z) WRITE(*,*) 'Fuzzy Comparisons: Y=Z'
        IF(Y.NE.Z) WRITE(*,*) 'Fuzzy Comparisons: Y<>Z'
C---But Y is tolerantly (and algebraically) equal to Z.
        IF(TEQ(Y,Z)) THEN
          WRITE(*,*) 'but TEQ(Y,Z) is .TRUE.'
          WRITE(*,*) 'That is, Y is computationally equal to Z.'
        ENDIF
        IF(TNE(Y,Z)) WRITE(*,*) 'and TNE(Y,Z) is .TRUE.'
      WRITE(*,*) ' '
C---Evaluate Fuzzy FLOOR and CEILing Function values using a Comparison
C   Tolerance, CT, of EPS3.
      X=0.11D0
      Y=((X*11.D0)-X)-0.1D0+2*EPSILON(X)
      YFLOOR=TFLOOR(Y,EPS3)
      YCEIL=TCEIL(Y,EPS3)
55    Z=1.D0
      WRITE(*,*) 'X=0.11D0, Y=X*11.D0-X-0.1D0, Z=1.D0'
      WRITE(*,*) 'X=',X,' Y=',Y,' Z=',Z
      WRITE(*,3) X,Y,Z
C---Floating-point Y is not identical (.EQ.) to floating-point Z.
      IF(Y.EQ.Z) WRITE(*,*) 'Fuzzy FLOOR/CEIL: Y=Z'
      IF(Y.NE.Z) WRITE(*,*) 'Fuzzy FLOOR/CEIL: Y<>Z'
      IF(TFLOOR(Y,EPS3).EQ.TCEIL(Y,EPS3).AND.TFLOOR(Y,EPS3).EQ.Z) THEN
C---But Tolerant Floor/Ceil of Y is identical (and algebraically equal)
C   to Z.
        WRITE(*,*) 'but TFLOOR(Y,EPS3)=TCEIL(Y,EPS3)=Z.'
        WRITE(*,*) 'That is, TFLOOR/TCEIL return exact whole numbers.'
      ENDIF
      STOP
      END
      DOUBLE PRECISION FUNCTION DUMACH ()
C***BEGIN PROLOGUE  DUMACH
C***PURPOSE  Compute the unit roundoff of the machine.
C***CATEGORY  R1
C***TYPE      DOUBLE PRECISION (RUMACH-S, DUMACH-D)
C***KEYWORDS  MACHINE CONSTANTS
C***AUTHOR  Hindmarsh, Alan C., (LLNL)
C***DESCRIPTION
C *Usage:
C        DOUBLE PRECISION  A, DUMACH
C        A = DUMACH()
C
C *Function Return Values:
C     A : the unit roundoff of the machine.
C
C *Description:
C     The unit roundoff is defined as the smallest positive machine
C     number u such that  1.0 + u .ne. 1.0.  This is computed by DUMACH
C     in a machine-independent manner.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  DUMSUM
C***REVISION HISTORY  (YYYYMMDD)
C   19930216  DATE WRITTEN
C   19930818  Added SLATEC-format prologue.  (FNF)
C   20030707  Added DUMSUM to force normal storage of COMP.  (ACH)
C***END PROLOGUE  DUMACH
C
      DOUBLE PRECISION U, COMP
      U = 1.0D0
 10   U = U*0.5D0
      CALL DUMSUM(1.0D0, U, COMP)
      IF (COMP .NE. 1.0D0) GO TO 10
      DUMACH = U*2.0D0
      RETURN
      END
      SUBROUTINE DUMSUM(A,B,C)
C     Routine to force normal storing of A + B, for DUMACH.
      DOUBLE PRECISION A, B, C
      C = A + B
      RETURN
      END
      DOUBLE PRECISION FUNCTION TFLOOR(X,CT)
C===========Tolerant FLOOR Function.
C
C    C  -  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 took five
C        years of refereed evolution (publication).
C
C    {2} L. M. Breed, "Definitions for Fuzzy Floor and Ceiling",  APL
C        QUOTE QUAD 8(3):16-23, March 1978.
C
C   H. D. KNOBLE, Penn State University.
C=====================================================================
      DOUBLE PRECISION X,Q,RMAX,EPS5,CT,FLOOR,DINT
C---------FLOOR(X) is the largest integer algegraically 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+DMAX1(CT,DMIN1(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
      DOUBLE PRECISION FUNCTION ROUND(X,CT)
C=========Tolerant Round 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