!**********************************************************************
! This file: http://ftp.aset.psu.edu/pub/ger/fortran/hdk/eps.f90
!  ROUTINE:   FUZZY FORTRAN OPERATORS
!  PURPOSE:   Illustrate Fortran implementation of the APL
!             tolerant comparisons, tolerant CEIL/FLOOR, and Tolerant
!             ROUND functions.
!  PLATFORM:  PC Windows/Linux Fortran compilers, F-Compiler, AIX XLF90
!  TO RUN:    F eps.f90 -o eps.exe
!             AIX: XLF90 eps.f -o eps.exe -qfloat=nomaf
!  CALLS:     none
!  AUTHOR:    H. D. Knoble <hdkLESS at SPAM psu dot edu>
!  DATE:      7 January 2004
!  REVISIONS:
!**********************************************************************
!
 module fuzzcomp
   implicit none
   integer, public, parameter :: PREC = selected_real_kind(15)
   public :: feq,fgt,fne,fle,flt,fge, tfloor,tceil,round
   real(kind=PREC), public :: one=1.0_PREC, fuzz=2*epsilon(one)

 contains

  function feq(x,y) result (test)
    real(kind=PREC), intent(in) :: x,y
    logical :: test
    test = abs(x-y) <= max(abs(x),abs(y))*fuzz
  end function feq

  function fgt(x,y) result (test)
    real(kind=PREC), intent(in) :: x,y
    logical :: test
    test = (x-y) > max(abs(x),abs(y))*fuzz
  end function fgt

  function fne(x,y) result (test)
    real(kind=PREC), intent(in) :: x,y
    logical :: test
    test = .not. feq(x,y)
  end function fne

  function fle(x,y) result (test)
    real(kind=PREC), intent(in) :: x,y
    logical :: test
    test = .not. fgt(x,y)
  end function fle

  function flt(x,y) result (test)
    real(kind=PREC), intent(in) :: x,y
    logical :: test
    test = fle(x,y) .and. fne(x,y)
  end function flt

  function fge(x,y) result (test)
    real(kind=PREC), intent(in) :: x,y
    logical :: test
    test = fgt(x,y) .or. feq(x,y)
  end function fge

  function tfloor(x) result(r)
 !===========Tolerant FLOOR Function.
 !
 !  x - is given as a double precision argument to be operated on.
 !      it is assumed that X is represented with m mantissa bits.
 !  fuzz - is given as a Comparison Tolerance such that
 !         0.lt.fuzz.le.3-Sqrt(5)/2. If the relative difference between
 !         X and a whole number is less than fuzz, then TFLOOR is
 !         returned as this whole number. By treating the
 !         floating-point numbers as a finite ordered set note that
 !         the heuristic eps=2.**(-(m-1)) and fuzz=2*eps  causes
 !         arguments of TFLOOR/TCEIL to be treated as whole numbers
 !         if they are exactly whole numbers or are immediately
 !         adjacent to whole number representations. Since EPS, the
 !         "distance" between floating-point numbers on the unit
 !         interval, and m, the number of bits in X's mantissa, exist
 !         on every floating-point computer, TFLOOR/TCEIL are
 !         consistently definable on every floating-point computer.
 !
 !         For more information see the following references:
 !    {1} P. E. Hagerty, "More on Fuzzy Floor and Ceiling," APL  QUOTE
 !        QUAD 8(4):20-24, June 1978. Note that TFLOOR=FL5 took five
 !        years of refereed evolution (publication).
 !
 !    {2} L. M. Breed, "Definitions for Fuzzy Floor and Ceiling",  APL
 !        QUOTE QUAD 8(3):16-23, March 1978.
 !
 !   H. D. KNOBLE, Penn State University - 7 January 2004
 !=====================================================================
       real (kind=PREC), intent(in) :: x
       real (kind=PREC) :: q,rmax,eps5,r
 !---------Hagerty's FL5 Function follows...
       q=one
       if( x<0 ) then
         q=one-fuzz
       endif
       rmax=q/(2.0_PREC-fuzz)
       eps5=fuzz/q
       r=floor(x+max(fuzz,min(rmax,eps5*abs(one+floor(x)))))
       if(.not.(x <=0 .or. (r-x) < rmax)) then
         r=r-one
       endif
  end function tfloor

  function tceil(x) result(r)
 !==========Tolerant Ceiling Function.
      real(kind=PREC), intent(in) :: x
      real(kind=PREC) :: r
      r= -tfloor(-x)
  end function tceil

  function round(x) result(r)
      real(kind=PREC), intent(in) :: x
      real(kind=PREC) :: r
 !=========Tolerant Round Function
 !  See Knuth, Art of Computer Programming, Vol. 1, Problem 1.2.4-5.
       r=tfloor(x+0.5_PREC)
  end function round
 end module fuzzcomp

 program fuzztest
! Example of Fuzzy comparison testing.
 use fuzzcomp
 integer :: i
 real(kind=PREC) :: x,y,z

! Example: x <> y but x feq y.
      x=one
      i=-49
      z=one/i
      y=z*i

! Fuzzy comparisons
      if (x /= y) then
       write(unit=*,fmt="(1X,A,2G25.17)") "x,y=",x,y
        print *,"Fuzzy Equality test: x /= y;"
      endif
      if (fgt(x,y)) then
        print *, "but x is fuzzily > y"
      else if(feq(x,y)) then
        print *, "but x is fuzzily = y"
      else
        print *, "but x is fuzzily < y"
      endif
      print *, " "

!---Evaluate Fuzzy FLOOR and CEILing Function values using a Comparison
!   Tolerance, (global) fuzz.
      x=0.11_prec
      y=((x*11.0_prec)-x)-0.1_prec+2*epsilon(x)
      z=one
      write(unit=*,fmt="(1X,A,2G25.17)") "Y,Z=",y,z
!---Floating-point Y /= to floating-point Z.
      if(y == z) then
        print *, "fuzzy FLOOR/CEIL test: Y==Z"
      else if(y /= z) then
        print *, "fuzzy FLOOR/CEIL test: Y/=Z"
      endif
      if(tfloor(y) == tceil(y).and.tfloor(y) == z) then
!---But Tolerant Floor/Ceil of Y is identical to Z.
        print *, "but TFLOOR(Y)=TCEIL(Y)==Z."
        print *, "That is, TFLOOR/TCEIL return exact whole numbers."
      else
        print *, "and TFLOOR(Y)=TCEIL(Y)/=Z."
        print *, "TFLOOR/TCEIL did NOT return exact whole numbers."
      endif

! Intel x86 output:
! x,y=   1.0000000000000000  0.99999999999999989
! Fuzzy Equality test: x /= y;
! but x is fuzzily = y
!
! Y,Z=  1.0000000000000004  1.0000000000000000
! fuzzy FLOOR/CEIL test: Y/=Z
! but TFLOOR(Y)=TCEIL(Y)==Z.
! That is, TFLOOR/TCEIL return exact whole numbers.

end program fuzztest