Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

A function to compare two floating-point numbers #145

Open
masuday opened this issue Feb 6, 2020 · 9 comments
Open

A function to compare two floating-point numbers #145

masuday opened this issue Feb 6, 2020 · 9 comments
Labels
topic: mathematics linear algebra, sparse matrices, special functions, FFT, random numbers, statistics, ...

Comments

@masuday
Copy link

masuday commented Feb 6, 2020

Inspied by Julia, I would like to have a function to compare two floating-point numbers with a tolerance, useful for testing.

elemental function is_approximated_1(x) result(nearly_equal)
    real(dp),intent(in) :: x
    logical :: nearly_equal
    nearly_equal = abs(x) < d_tol
end function is_approximated_1

elemental function is_approximated_2(x,y) result(nearly_equal)
    real(dp),intent(in) :: x,y
    logical :: nearly_equal
    nearly_equal = abs(x-y) < d_tol
end function is_approximated_2

elemental function is_approximated_2_tol(x,y,tol) result(nearly_equal)
    real(dp),intent(in) :: x,y,tol
    logical :: nearly_equal
    nearly_equal = abs(x-y) < tol
end function is_approximated_2

Defining a general name:

interface is_approximated
   module procedure is_approximated_1, is_approximated_2, is_approximated_2_tol
end interface is_approximated

Usage:

real(real64) :: x(5),y(5)

print *,any(is_approximated(x,y))
@jvdp1
Copy link
Member

jvdp1 commented Feb 6, 2020

I would propose the following API:

is_approximated( x [, tol])
is_approximated( x, y, [, tol])

where tol is optional for both functions. A default tol can be determined with optval.
So, something like:

elemental function is_approximated_2_tol_dp(x,y,tol) result(nearly_equal)
    real(dp),intent(in) :: x,y
    real(dp), intent(in), optional :: tol
    logical :: nearly_equal
    nearly_equal = abs(x-y) < optval(tol, d_tol)
end function is_approximated_2_tol_dp

I would also think that we can drop is_approaximated(x [, tol]) because the user would do explicitely what the second API does implicitely.

This function could be used in the tests of stdlib. A bitwise check with (almost) the same API could also be implemented. See #142 for more discussion.

Finally I may envisage some rebuttals from the community regarding this function because it is a one-line function.

@masuday
Copy link
Author

masuday commented Feb 6, 2020 •
edited
Loading

I agree with your suggestion. A simple API is good enough.

I like this function because this comparison always occurs in any tests. A function with a clear purpose makes the code readable. For more precise comparisons, a bit-wise approach is preferred.

@aradi
Copy link
Member

aradi commented Feb 6, 2020

I'd suggest to go with something more general, which allows also relative error. I like quite a lot Numpy's isclose() function which enables to specify both, absolute and relative error. So the interface could be like

elemental function is_close(x, y, reltol, abstol) result(isclosel)
    real(dp),intent(in) :: x, y
    real(dp), intent(in), optional :: reltol, abstol
    logical :: isclose
    isclose = abs(x - y) < reltol * abs(y) + atol
end function is_close

Also, the name is_close() is much shorter than is_approximated().

@masuday
Copy link
Author

masuday commented Feb 6, 2020

@aradi Julia also uses a similar function. Personally, for a quick comparison, I am using the following function (as suggested by Tamas_Papp).

abs(x - y) < reltol * (1+abs(x)+abs(y))

@wclodius2
Copy link
Contributor

I suggest is_within.

@ivan-pi ivan-pi mentioned this issue Feb 26, 2021
Open
@ivan-pi ivan-pi mentioned this issue Aug 19, 2021
Merged
6 tasks
@awvwgk awvwgk added the topic: mathematics linear algebra, sparse matrices, special functions, FFT, random numbers, statistics, ... label Sep 18, 2021
@ivan-pi
Copy link
Member

ivan-pi commented Sep 27, 2021

In #488 @bilderbuchi mentions the Python math.isclose as another possible implementation. The choices behind the Python version (available with Python >= 3.5) are discussed in a Python Enhancement Proposal PEP 485 which also considers the numpy and Boost C++ implementations.

Tagging @zoziha, @milancurcic

@certik
Copy link
Member

certik commented Oct 29, 2024 •
edited
Loading

It looks like some good choices are:

  • Python: abs(a-b) <= max(rtol*max(abs(a), abs(b)), atol)
  • NumPy: abs(a-b) <= atol + rtol*abs(b)
  • Tamas Papp: abs(a-b) <= tol*(1 + abs(a) + abs(b))

They all seem fundamentally the same and all follow from the Python version using the relation max(abs(a), abs(b)) <= abs(a) + abs(b):

max(rtol*max(abs(a), abs(b)), atol) <= max(rtol*(abs(a) + abs(b)), atol) <= atol + rtol*(abs(a) + abs(b)) <= max(atol, rtol)*(1+abs(a)+abs(b))

So for example:

  • If Python is satisfied, then also Tamas Papp (but not necessarily vice versa) with tol = max(atol, rtol)
  • Assuming abs(b) >= abs(a), if Python is satisfied, then so is NumPy (but not vice versa)

Based on my experience, when comparing some results of numerical computation which happens to be around 0, it can easily be that all digits are wrong, just the order of magnitude (1e-16) is correct, so the rtol is meaningless, we want to use atol=1e-14 (let's say) to get a good comparison. When numbers become large, say 1e8, then the atol is meaningless, however the significant digits will usually be mostly correct, and the rtol measures exactly that (for example rtol=1e-10 means roughly 10 significant digits are correct). So we want to use abs(a-b) <= rtol*max(abs(a), abs(b)) for large numbers, and abs(a-b) <= atol for small numbers, and it turns out the expression max(rtol*max(abs(a), abs(b)), atol) combines the two in the correct way: for large numbers it is equal to rtol*max(abs(a), abs(b)), for small numbers it is equal to atol. The Tamas Papp's expression is a relaxation of it, but even that still works well: for small numbers 1 + abs(a) + abs(b) -> 1 and tol becomes atol. For large numbers 1+abs(a)+abs(b) ->2*max(abs(a), abs(b)) and 2*tol becomes rtol.

@arjenmarkus
Copy link
Member

arjenmarkus commented Nov 2, 2024 via email •
edited by certik
Loading

@Beliavsky
Copy link

A current project is https://github.com/suzuyuyuyu/fortran-approx

"This repository distributes approx.f90. It is a simple Fortran program that enables approximate comparisons, such as .eq., .ne., .lt., .le., .gt., and .ge.. The tolerance is set to 1.0e-7 by default, but it can be changed by the user—either by modifying the module or by passing an argument to the function.

In addition, custom operators such as .aeq., .ane., .alt., .ale., .agt., and .age. are defined. These operators perform comparisons that take the specified tolerance into account. The same tolerance value is used for these operators, and it can only be changed by modifying the module."

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
topic: mathematics linear algebra, sparse matrices, special functions, FFT, random numbers, statistics, ...
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
9 participants
@certik @arjenmarkus @aradi @masuday @jvdp1 @ivan-pi @awvwgk @Beliavsky @wclodius2