Skip to content

Instantly share code, notes, and snippets.

@certik certik/a.f90
Created Sep 24, 2012

Embed
What would you like to do?
Donald Knuth has a proposal for a floating point comparison algorithm
logical elemental function fequal(a, b, eps) result(r)
! Returns .true. if a==b (withing "eps"), otherwise .false.
real(dp), intent(in) :: a, b, eps
integer :: e
if (abs(a) > abs(b)) then
e = exponent(a)
else
e = exponent(b)
end if
if (e < minexponent(1._dp)) e = minexponent(1._dp)
r = abs(a-b) < eps * 2._dp**e
end function
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.