From: HDK@psuvm.psu.edu (H. D. Knoble)
Newsgroups: comp.lang.fortran
Subject: Re: Ignoring Least Significant Bit in a Real Compare
Date: 21 Aug 92 22:13:00 GMT
 <1992Aug21.212225.2961@news.eng.convex.com>
Organization: Penn State University


                Floating-Point System Environment Parameters
                       and Tolerant (Fuzzy) Comparisons
                       H. D. Knoble, HDK@PSUVM.PSU.EDU
                  Penn State Center for Academic Computing

Fortran algorithms are given that compute tolerant (fuzzy) comparisons of
real numbers and tolerant (fuzzy) Floor/Ceiling Functions.  The driver
FUZZY FORTRAN illustrates use of these tolerant functions; it also shows
to compute 13 floating-point parameters associated with the
floating-point system being utilized by the Fortran host computer.
(Definitions are given in comments of Cody & Waite's Subroutine MACHAR).

As a programming practice, one should not compare real numbers for
equality, as quantities computed with different but algebraically
equivalent formulae will hardly ever be computed equal; yet they may be
equal within a floating-point representation (i.e., plus or minus one bit
in the last place of the mantissa).  Tolerant comparisons enable Fortran
algorithms to compare real numbers within this "tolerance" which is
sometimes called the "machine precision" (distance between floating-point
numbers on the unit interval).  Similarly when computing the greatest
integer less than or equal to an argument X, one would not want for
example, 1 for FLOOR(X=1.999999999999999).  Tolerant Floor/Ceiling
Functions resolve this problem in a most sophisticated way (equivalent to
the algorithms used in APL).

Algorithms included in FUZZY FORTRAN include:

Driver Program    Illustrates use of following subprograms as well as
                  cases where quantities are not equal with respect to
                  .EQ. but are with respect to TEQ; similarly for
                  Floor and TFLOOR. This may be IBM Fortran dependent
                  in that Z-Format (Hexadecimal) is used to display
                  some results.

TEQ(X,Y), TNE(X,Y)  Logical Arithmetic Statement Functions to tolerantly
                    compare Double Precision quantities X and Y.

TFLOOR(X,CT)  Returns Double Precision whole number that is tolerantly
              the greatest integer algebraically less than X; CT is the
              Comparison Tolerance. This function, small as it is too
              several years to evolve in the APL community.

TCEIL(X,CT)   Returns Double Precision whole number that is tolerantly
              the smallest integer algebraically greater than X; CT is
              the Comparison Tolerance.

ROUND(X,CT)   Returns Double Precision tolerant proper rounding of X;
              CT is the Comparison Tolerance.

MACHAR( parms )   Returns 13 Characteristics of the Machine's
                  floating-point system. Cody and Waite's work.

C**********************************************************************
C  ROUTINE:   FUZZY FORTRAN
C  PURPOSE:   Illustrate Hindmarsh's computation of EPS, and APL        t
C             tolerant comparisons, tolerant CEIL/FLOOR, and Tolerant
C             ROUND functions - implemented in Fortran.
C  CALLS:     none
C  PACKAGE:   SAMPLES
C  PLATFORM:  Test system IBM VS Fortran on IBM 370 and up.
C  AUTHOR:    H. D. Knoble <HDK@PSUVM> 09/22/78
C  ORGANIZATION: Penn State University Center for Academic Computing
C  REVISIONS:
C**********************************************************************
C
      DOUBLE PRECISION EPS,EPS3, X,Y,Z, TFLOOR,TCEIL, EPSNEG,XMIN,XMAX
      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.
      CALL MACHAR (IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
     1             MAXEXP,EPS,EPSNEG,XMIN,XMAX)
      WRITE(6,*) 'Floating-point environment parameters follow...'
      WRITE(6,*) 'IBETA,IT,IRND=',IBETA,IT,IRND
      WRITE(6,*) 'MACHEP,NEGEP,IEXP=',MACHEP,NEGEP,IEXP
      WRITE(6,*) 'MINEXP,MAXEXP=',MINEXP,MAXEXP
      WRITE(6,*) 'EPS,EPSNEG=',EPS,EPSNEG
      WRITE(6,*) 'XMIN,XMAX=',XMIN,XMAX
C---Accept an exact representation, or one bit on either side
C   (heuristic).
      EPS3=3.D0*EPS
      WRITE(6,1) EPS,EPS, EPS3,EPS3
1     FORMAT(' EPS=',D16.8,2X,Z16, ', EPS3=',D16.8,2X,Z16)
      WRITE(6,*) ' '
C---Illustrate Tolerant Comparisons using EPS3. Any other magnitudes will
C   behave similarly.
      Z=1.D0
      I=49
        X=1.D0/I
        Y=X*I
        WRITE(6,*) 'X=1.D0/',I,', Y=X*',I,', Z=1.D0'
        IF(Y.EQ.Z) WRITE(6,*) 'Fuzzy Comparisons: Y=Z'
        IF(Y.NE.Z) WRITE(6,*) 'Fuzzy Comparisons: Y<>Z'
        IF(TEQ(Y,Z)) THEN
          WRITE(6,*) 'but TEQ(Y,Z) is .TRUE.'
          WRITE(6,*) 'That is, Y is computationally equal to Z.'
        ENDIF
        IF(TNE(Y,Z)) WRITE(6,*) 'and TNE(Y,Z) is .TRUE.'
      WRITE(6,*) ' '
C---Evaluate Fuzzy FLOOR and CEILing Function values using a Comparison
C   Tolerance, CT, of EPS3=3*EPS.
      X=0.1D0
      Y=X*11.D0-X
      YFLOOR=TFLOOR(Y,EPS3)
      YCEIL=TCEIL(Y,EPS3)
55    Z=1.D0
      WRITE(6,*) 'X=0.1D0, Y=X*11.D0-X, Z=1.D0'
      IF(Y.EQ.Z) WRITE(6,*) 'Fuzzy FLOOR/CEIL: Y=Z'
      IF(Y.NE.Z) WRITE(6,*) 'Fuzzy FLOOR/CEIL: Y<>Z'
      IF(TFLOOR(Y,EPS3).EQ.TCEIL(Y,EPS3).AND.TFLOOR(Y,EPS3).EQ.Z) THEN
        WRITE(6,*) 'but TFLOOR(Y,EPS3)=TCEIL(Y,EPS3)=Z.'
        WRITE(6,*) 'That is, TFLOOR/TCEIL return exact whole numbers.'
      ENDIF
      STOP
      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 bewteen
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.
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
       SUBROUTINE MACHAR (IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
     1                    MAXEXP,EPS,EPSNEG,XMIN,XMAX)
       IMPLICIT REAL*8 (A-H,O-Z)
C===See "SOFTWARE MANUAL FOR THE ELEMENTARY FUNCTIONS" by William J.
C   Cody, Jr. and William Waite, Prentice Hall, ISBN 0-13-822064,
C   1980, pages 258-264.
C
C   The subroutine MACHAR is intended to dynamically determine
C   thirteen parameters relating to the floating-point
C   arithmetic system.  The first three--the radix B, the number
C   of base-B digits in the significand, and the question of
C   whether floating-point addition rounds or truncates--are
C   determined with an algorithm developed by M.  Malcolm
C   [1972].  Several of the remaining parameters are determined
C   with "folk algorithms," but most of them are determined
C   using algorithms we believe to be new.  The methods are
C   sometimes ad hoc, and the program is not guaranteed to work
C   on all machines.  Not all of the conditions under which
C   MACHAR will fail are known, but it is known to fail on
C   Honeywell machines where the active arithmetic registers are
C   wider than storage registers.  This particular failure can
C   be corrected by forcing the storage of intermediate results
C   at critical points (see Gentleman and Marovich [1974]).  We
C   have deliberately not followed that procedure in MACHAR
C   because the program is already very complicated and
C   difficult to understand.
C
C   Aside from the Honeywell machines the program has run
C   properly on all other machines we have tried.  These
C   include, in alphabetical order, BESM-6 (Russian), Burroughs
C   6700, CDC 6000/7000 series, CRAY-1, IBM 360/370 series, IBM
C   3033, PDP-10, PDP-11, Univac 1105, Varian V76, and Z-80
C   (with Microsoft fortran).  MACHAR has been translated into
C   Basic and run on the Datapoint 2200, IBM 5100 and PDP-10.
C   It has even been translated and run on the HP-67 and TI-59
C   programmable calculators.  Despite these successes we
C   strongly urge that parameter values returned by this program
C   be examined critically the first time it is run on any
C   machine.
C
       INTEGER I,IBETA,IEXP,IRND,IT,IZ,J,K,MACHEP,MAXEXP,MINEXP,
     1         MX,NEGEP,NGRD
CE     REAL A,B,BETA,BETAIN,BETAM1,EPS,EPSNEG,ONE,XMAX,XMIN,Y,Z,ZERO
       DOUBLE PRECISION A,B,BETA,BETAIN,BETAM1,EPS,EPSNEG,ONE,XMAX,
     1                  XMIN,Y,Z,ZERO
C
C      This subroutine is intended to determine the characteristics
C      of the floating-point arithmetic system that are specified
C      below.  The first three are determined according to an
C      algorithm due to M. Malcolm, CACM 15 (1972), pp. 949-951,
C      incorporating some, but not all, of the improvements
C      suggested by M. Gentleman and S. Marovich, CACM 17 (1974),
C      pp. 276-277.  The version given here is for double precision.
C      Lines containing  CE  in columns 1 and 2 can be used to
C      convert the subroutine to single precision by replacing
C      existing lines in the obvious manner.
C
C
C      IBETA   - the radix of the floating-point representation
C      IT      - the number of base ibeta digits in the floating-point
C                significand
C      IRND    - 0 if floating-point addition chops,
C                1  if floating-point addition rounds
C      NGRD    - the number of guard digits for multiplication.  It is
C                0 if  IRND=1, or if  IRND=0  and only  it  base IBETA
C                  digits participate in the post normalization shift
C                  of the floating-point significand in multiplication
C                1 if  IRND=0  and more than  it  base  IBETA  digits
C                  participate in the post normalization shift of the
C                  floating-point significand in multiplication.
C      MACHEP  - the largest negative integer such that
C                1.0+FLOAT(IBETA)**MACHEP .NE. 1.0, except that
C                MACHEP is bounded below by  -(IT+3)
C      NEGEPS  - the largest negative integer such that
C                1.0-FLOAT(IBETA)**NEGEPS .NE. 1.0, except that
C                NEGEPS is bounded below by  -(IT+3)
C      IEXP    - the number of bits (decimal places if IBETA = 10)
C                reserved for the representation of the exponent
C                (including the bias or sign) of a floating-point
C                number
C      MINEXP  - the largest in magnitude negative integer such that
C                FLOAT(IBETA)**MINEXP is a positive floating-point
C                number
C      MAXEXP  - the largest positive integer exponent for a finite
C                floating-point number
C      EPS     - the smallest positive floating-point number such
C                that  1.0+EPS .ne.  1.0. in particular, if either
C                IBETA = 2  or  IRND = 0, EPS = FLOAT(IBETA)**MACHEP.
C                otherwise,  EPS =  (FLOAT(IBETA)**MACHEP)/2
C      EPSNEG  - a small positive floating-point number such that
C                1.0-EPSNEG .ne. 1.0. in particular, if IBETA = 2
C                or  IRND = 0, EPSNEG = FLOAT(IBETA)**NEGEPS.
C                otherwise,  EPSNEG = (IBETA**NEGEPS)/2.  Because
C                NEGEPS is bounded below by -(IT+3), EPSNEG may not
C                be the smallest number which can alter 1.0 by
C                subtraction.
C      XMIN    - the smallest non-vanishing floating-point power of
C                the radix. in particular, XMIN = FLOAT(IBETA)**MINEXP
C      XMAX    - the largest finite floating-point number.  In
C                particular   XMAX = (1.0-EPSNEG)*FLOAT(IBETA)**MAXEXP
C                Note - on some machines  XMAX  will be only the
C                second, or perhaps third, largest number, being
C                too small by 1 or 2 units in the last digit of
C                the significand.
C
C      Latest revision - October 22, 1979
C
C      Author - W. J. Cody
C               Argonne National Laboratory
C
C-----------------------------------------------------------------
CE     ONE = FLOAT(1)
       ONE = DBLE(FLOAT(1))
CE     ZERO = 0.0E0
       ZERO = 0.0D0
C-----------------------------------------------------------------
C      Determine IBETA,BETA a la Malcolm
C-----------------------------------------------------------------
       A = ONE
   10  A = A + A
         IF (((A+ONE)-A)-ONE .EQ. ZERO) GO TO 10
       B = ONE
   20  B = B + B
         IF ((A+B)-A .EQ. ZERO) GO TO 20
CE    IBETA = INT((A+B)-A)
      IBETA = INT(SNGL((A + B) - A))
CE    BETA = FLOAT(IBETA)
      BETA = DBLE(FLOAT(IBETA))

C     Determine IT, IRND
C-----------------------------------------------------------------
      IT = 0
      B  = ONE
  100 IT = IT + 1
         B = B * BETA
         IF (((B+ONE)-B)-ONE .EQ. ZERO) GO TO 100
      IRND = 0
      BETAM1 = BETA - ONE
      IF ((A+BETAM1)-A .NE. ZERO) IRND = 1
C-----------------------------------------------------------------
C     Determine NEGEP, EPSNEG
C-----------------------------------------------------------------
      NEGEP = IT + 3
      BETAIN = ONE / BETA
      A  = ONE
C
      DO 200 I = 1, NEGEP
         A = A * BETAIN
  200 CONTINUE
C
      B = A
  210 IF ((ONE-A)-ONE .NE. ZERO) GO TO 220
         A = A * BETA
         NEGEP = NEGEP - 1
      GO TO 210
  220 NEGEP = -NEGEP
      EPSNEG = A
      IF ((IBETA .EQ. 2) .OR. (IRND .EQ. 0)) GO TO 300
      A = (A*(ONE+A)) / (ONE+ONE)
      IF ((ONE-A)-ONE .NE. ZERO) EPSNEG = A

C     Determine MACHEP, EPS
C-----------------------------------------------------------------
  300 MACHEP = -IT - 3
      A = B
  310 IF((ONE+A)-ONE .NE. ZERO) GO TO 320
          A = A * BETA
          MACHEP = MACHEP + 1
       GO TO 310
  320  EPS = A
       IF ((IBETA .EQ. 2) .OR. (IRND .EQ. 0)) GO TO 350
       A = (A*(ONE+A)) / (ONE+ONE)
       IF ((ONE+A)-ONE .NE. ZERO) EPS = A
C-----------------------------------------------------------------
C      Determine NGRD
C-----------------------------------------------------------------
  350  NGRD = 0
       IF ((IRND .EQ. 0) .AND. ((ONE+EPS)*ONE-ONE) .NE. ZERO) NGRD = 1
C-----------------------------------------------------------------
C      Determine IEXP, MINEXP, XMIN
C
C      Loop to determine largest I and K = 2**I such that
C         (1/BETA) ** (2**(1))
C      does not underflow.
C      Exit from loop is signaled by an underflow.
C-----------------------------------------------------------------
C---VS Fortran: allow two underflows with no messages.
C      CALL ERRSET(208,3,-1,+1)
       I = 0
       K = 1
       Z = BETAIN
  400  Y = Z
          Z = Y * Y
C-----------------------------------------------------------------
C         Check for Underflow here
C-----------------------------------------------------------------
          A = Z * ONE
CE        IF ((A+A .EQ. ZERO) .OR.  (ABS(Z) .GE. Y)) GO TO 410
          IF ((A+A .EQ. ZERO) .OR.  (DABS(Z) .GE. Y)) GO TO 410
          I = I + 1
          K = K + K
       GO TO 400
  410  IF (IBETA .EQ.  10) GO TO 420
       IEXP = I + 1
       MX = K + K
       GO TO 450
C-----------------------------------------------------------------
C      For decimal machines only
C-----------------------------------------------------------------
  420  IEXP = 2
       IZ = IBETA
  430  IF (K .LT. IZ) GO TO 440
          IZ = IZ * IBETA
          IEXP = IEXP + 1
       GO TO 430
  440  MX = IZ + IZ - 1
C-----------------------------------------------------------------
C      Loop to determine MINEXP, XMIN
C      Exit from loop is signaled by an underflow.
C-----------------------------------------------------------------
  450  XMIN = Y
          Y = Y * BETAIN
C-----------------------------------------------------------------
C         Check for underflow here
C-----------------------------------------------------------------
          A = Y * ONE
CE        IF (((A+A) .EQ. ZERO) .OR. (ABS(Y) .GE. XMIN)) GO TO 460
          IF (((A+A) .EQ. ZERO) .OR. (DABS(Y)  .GE. XMIN)) GO TO 460
          K = K + 1
       GO TO 450
  460 MINEXP = -K
C-----------------------------------------------------------------
C      Determine MAXEXP, XMAX
C-----------------------------------------------------------------
       IF ((MX .GT. K+K-3) .OR.  (IBETA .EQ. 10)) GO TO 500
       MX = MX + MX
       IEXP = IEXP + 1
  500  MAXEXP = MX + MINEXP
C-----------------------------------------------------------------
C      Adjust for machines with implicit leading
C      bit in binary significant and machines with
C      radix point at extreme right of significand
C-----------------------------------------------------------------
       I = MAXEXP + MINEXP
       IF ((IBETA .EQ. 2) .AND.  (I .EQ. 0)) MAXEXP = MAXEXP - 1
       IF (I .GT. 20) MAXEXP = MAXEXP - 1
       IF (A .NE. Y) MAXEXP = MAXEXP - 2
       XMAX = ONE - EPSNEG
       IF (XMAX*ONE .NE. XMAX) XMAX = ONE - BETA * EPSNEG
       XMAX = XMAX / (BETA * BETA * BETA * XMIN)
       I = MAXEXP + MINEXP + 3
       IF (I .LE. 0) GO TO 520
C
       DO 510 J = 1, I
           IF (IBETA .NE. 2) XMAX = XMAX * BETA
  510  CONTINUE
C
  520  RETURN
C     ---------- Last line of MACHAR ----------
      END
Newsgroups: comp.lang.fortran