1!********************************************************************** 2! ROUTINE: FUZZY FORTRAN OPERATORS 3! PURPOSE: Illustrate Hindmarsh's computation of EPS, and APL 4! tolerant comparisons, tolerant CEIL/FLOOR, and Tolerant 5! ROUND functions - implemented in Fortran. 6! PLATFORM: PC Windows Fortran, Compaq-Digital CVF 6.1a, AIX XLF90 7! TO RUN: Windows: DF EPS.F90 8! AIX: XLF90 eps.f -o eps.exe -qfloat=nomaf 9! CALLS: none 10! AUTHOR: H. D. Knoble <hdk@psu.edu> 22 September 1978 11! REVISIONS: 12!********************************************************************** 13! 14 DOUBLE PRECISION EPS,EPS3, X,Y,Z, D1MACH,TFLOOR,TCEIL,EPSF90 15 LOGICAL TEQ,TNE,TGT,TGE,TLT,TLE 16!---Following are Fuzzy Comparison (arithmetic statement) Functions. 17! 18 TEQ(X,Y)=DABS(X-Y).LE.DMAX1(DABS(X),DABS(Y))*EPS3 19 TNE(X,Y)=.NOT.TEQ(X,Y) 20 TGT(X,Y)=(X-Y).GT.DMAX1(DABS(X),DABS(Y))*EPS3 21 TLE(X,Y)=.NOT.TGT(X,Y) 22 TLT(X,Y)=TLE(X,Y).AND.TNE(X,Y) 23 TGE(X,Y)=TGT(X,Y).OR.TEQ(X,Y) 24! 25!---Compute EPS for this computer. EPS is the smallest real number on 26! this architecture such that 1+EPS>1 and 1-EPS<1. 27! EPSILON(X) is a Fortran 90 built-in Intrinsic function. They should 28! be identically equal. 29! 30 EPS=D1MACH(NULL) 31 EPSF90=EPSILON(X) 32 IF(EPS.NE.EPSF90) THEN 33 WRITE(*,2)'EPS=',EPS,' .NE. EPSF90=',EPSF90 342 FORMAT(A,Z16,A,Z16) 35 ENDIF 36!---Accept a representation if exact, or one bit on either side. 37 EPS3=3.D0*EPS 38 WRITE(*,1) EPS,EPS, EPS3,EPS3 391 FORMAT(' EPS=',D16.8,2X,Z16, ', EPS3=',D16.8,2X,Z16) 40!---Illustrate Fuzzy Comparisons using EPS3. Any other magnitudes will 41! behave similarly. 42 Z=1.D0 43 I=49 44 X=1.D0/I 45 Y=X*I 46 WRITE(*,*) 'X=1.D0/',I,', Y=X*',I,', Z=1.D0' 47 WRITE(*,*) 'Y=',Y,' Z=',Z 48 WRITE(*,3) X,Y,Z 493 FORMAT(' X=',Z16,' Y=',Z16,' Z=',Z16) 50!---Floating-point Y is not identical (.EQ.) to floating-point Z. 51 IF(Y.EQ.Z) WRITE(*,*) 'Fuzzy Comparisons: Y=Z' 52 IF(Y.NE.Z) WRITE(*,*) 'Fuzzy Comparisons: Y<>Z' 53!---But Y is tolerantly (and algebraically) equal to Z. 54 IF(TEQ(Y,Z)) THEN 55 WRITE(*,*) 'but TEQ(Y,Z) is .TRUE.' 56 WRITE(*,*) 'That is, Y is computationally equal to Z.' 57 ENDIF 58 IF(TNE(Y,Z)) WRITE(*,*) 'and TNE(Y,Z) is .TRUE.' 59 WRITE(*,*) ' ' 60!---Evaluate Fuzzy FLOOR and CEILing Function values using a Comparison 61! Tolerance, CT, of EPS3. 62 X=0.11D0 63 Y=((X*11.D0)-X)-0.1D0 64 YFLOOR=TFLOOR(Y,EPS3) 65 YCEIL=TCEIL(Y,EPS3) 6655 Z=1.D0 67 WRITE(*,*) 'X=0.11D0, Y=X*11.D0-X-0.1D0, Z=1.D0' 68 WRITE(*,*) 'X=',X,' Y=',Y,' Z=',Z 69 WRITE(*,3) X,Y,Z 70!---Floating-point Y is not identical (.EQ.) to floating-point Z. 71 IF(Y.EQ.Z) WRITE(*,*) 'Fuzzy FLOOR/CEIL: Y=Z' 72 IF(Y.NE.Z) WRITE(*,*) 'Fuzzy FLOOR/CEIL: Y<>Z' 73 IF(TFLOOR(Y,EPS3).EQ.TCEIL(Y,EPS3).AND.TFLOOR(Y,EPS3).EQ.Z) THEN 74!---But Tolerant Floor/Ceil of Y is identical (and algebraically equal) 75! to Z. 76 WRITE(*,*) 'but TFLOOR(Y,EPS3)=TCEIL(Y,EPS3)=Z.' 77 WRITE(*,*) 'That is, TFLOOR/TCEIL return exact whole numbers.' 78 ENDIF 79 STOP 80 END 81 DOUBLE PRECISION FUNCTION D1MACH (IDUM) 82 INTEGER IDUM 83!======================================================================= 84! This routine computes the unit roundoff of the machine in double 85! precision. This is defined as the smallest positive machine real 86! number, EPS, such that (1.0D0+EPS > 1.0D0) & (1.D0-EPS < 1.D0). 87! This computation of EPS is the work of Alan C. Hindmarsh. 88! For computation of Machine Parameters also see: 89! W. J. Cody, "MACHAR: A subroutine to dynamically determine machine 90! parameters, " TOMS 14, December, 1988; or 91! Alan C. Hindmarsh at http://www.netlib.org/lapack/util/dlamch.f 92! or Werner W. Schulz at http://www.ozemail.com.au/~milleraj/ . 93! 94! This routine appears to give bit-for-bit the same results as 95! the Intrinsic function EPSILON(x) for x single or double precision. 96! hdk - 25 August 1999. 97!----------------------------------------------------------------------- 98 DOUBLE PRECISION EPS, COMP 99! EPS = 1.0D0 100!10 EPS = EPS*0.5D0 101! COMP = 1.0D0 + EPS 102! IF (COMP .NE. 1.0D0) GO TO 10 103! D1MACH = EPS*2.0D0 104 EPS = 1.0D0 105 COMP = 2.0D0 106 DO WHILE ( COMP .NE. 1.0D0 ) 107 EPS = EPS*0.5D0 108 COMP = 1.0D0 + EPS 109 ENDDO 110 D1MACH = EPS*2.0D0 111 RETURN 112 END 113 DOUBLE PRECISION FUNCTION TFLOOR(X,CT) 114!===========Tolerant FLOOR Function. 115! 116! C - is given as a double precision argument to be operated on. 117! it is assumed that X is represented with m mantissa bits. 118! CT - is given as a Comparison Tolerance such that 119! 0.lt.CT.le.3-Sqrt(5)/2. If the relative difference between 120! X and a whole number is less than CT, then TFLOOR is 121! returned as this whole number. By treating the 122! floating-point numbers as a finite ordered set note that 123! the heuristic eps=2.**(-(m-1)) and CT=3*eps causes 124! arguments of TFLOOR/TCEIL to be treated as whole numbers 125! if they are exactly whole numbers or are immediately 126! adjacent to whole number representations. Since EPS, the 127! "distance" between floating-point numbers on the unit 128! interval, and m, the number of bits in X's mantissa, exist 129! on every floating-point computer, TFLOOR/TCEIL are 130! consistently definable on every floating-point computer. 131! 132! For more information see the following references: 133! {1} P. E. Hagerty, "More on Fuzzy Floor and Ceiling," APL QUOTE 134! QUAD 8(4):20-24, June 1978. Note that TFLOOR=FL5 took five 135! years of refereed evolution (publication). 136! 137! {2} L. M. Breed, "Definitions for Fuzzy Floor and Ceiling", APL 138! QUOTE QUAD 8(3):16-23, March 1978. 139! 140! H. D. KNOBLE, Penn State University. 141!===================================================================== 142 DOUBLE PRECISION X,Q,RMAX,EPS5,CT,FLOOR,DINT 143!---------FLOOR(X) is the largest integer algegraically less than 144! or equal to X; that is, the unfuzzy Floor Function. 145 DINT(X)=X-DMOD(X,1.D0) 146 FLOOR(X)=DINT(X)-DMOD(2.D0+DSIGN(1.D0,X),3.D0) 147!---------Hagerty's FL5 Function follows... 148 Q=1.D0 149 IF(X.LT.0)Q=1.D0-CT 150 RMAX=Q/(2.D0-CT) 151 EPS5=CT/Q 152 TFLOOR=FLOOR(X+DMAX1(CT,DMIN1(RMAX,EPS5*DABS(1.D0+FLOOR(X))))) 153 IF(X.LE.0 .OR. (TFLOOR-X).LT.RMAX)RETURN 154 TFLOOR=TFLOOR-1.D0 155 RETURN 156 END 157 DOUBLE PRECISION FUNCTION TCEIL(X,CT) 158!==========Tolerant Ceiling Function. 159! See TFLOOR. 160 DOUBLE PRECISION X,CT,TFLOOR 161 TCEIL= -TFLOOR(-X,CT) 162 RETURN 163 END 164 DOUBLE PRECISION FUNCTION ROUND(X,CT) 165!=========Tolerant Round Function 166! See Knuth, Art of Computer Programming, Vol. 1, Problem 1.2.4-5. 167 DOUBLE PRECISION TFLOOR,X,CT 168 ROUND=TFLOOR(X+0.5D0,CT) 169 RETURN 170 END 171