Tolerant comparison of real numbers

[arjenmarkus]

I have the source for a set of comparison functions and other functions such as floor and ceil that take into account that real numbers are not as well-behaved as, say, integers, by taking a margin into account. The source comes from Skip Knoble, once quite active on comp.lang.fortran, and builds on work by Hindmarsh and others.

Is this something useful for stdlib? The source needs to be updated to the current standard, but I can do that :).

[epagone]

Personally, I would be very interested in that. The naive approach that I currently use in my codes is the following:

elemental logical function eqreal(real1, real2)
    real(wp), intent(in) :: real1, real2
    eqreal = abs(real1 - real2) < spacing( max( abs(real1), abs(real2) ) )
end function eqreal

I would be curious to learn what a proper approach would be.

[arjenmarkus]

I am not sure if there is a proper approach ;). The best I can think of is a well-chosen tolerance and that may actually be dependent on the application. The code I was referring predates Fortran 90, by the way, but it is definitely worth considering, if only because it will provide some common ground based on much prior work.

[ivan-pi]

I think this is a partial duplicate of https://github.com/fortran-lang/stdlib/issues/145.

We can keep this issue open for the well-behaved floor and ceil functions, but I would suggest changing the title in that case.

[arjenmarkus]

Indeed, the collection of functions I have includes the ones described there.

Op vr 26 feb. 2021 om 17:35 schreef Ivan Pribec [email protected]:

I think this is a partial duplicate of #145 https://github.com/fortran-lang/stdlib/issues/145.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/fortran-lang/stdlib/issues/327#issuecomment-786754837, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN6YRYLEJBQDN2OKVMV6RTTA7ETHANCNFSM4YIL5PBA .

[jvdp1]

Is this something useful for stdlib?

I think it is indeed something for stdlib. It could even be used inside stdlib for the different tests with real values (and for which we already got multiple issues/discussions).

[arjenmarkus]

I attach the original source code eps_f90.txt - I have not take the opportunity yet to convert it to a modern package, so it is really old-school :). However, the output from the program may be of interest to give better insight:

EPS= 0.22204460D-15 3CB0000000000000, EPS3= 0.66613381D-15 3CC8000000000000 X=1.D0/ 49 , Y=X* 49 , Z=1.D0 Y= 0.99999999999999989 Z= 1.0000000000000000
X=3F94E5E0A72F0539 Y=3FEFFFFFFFFFFFFF Z=3FF0000000000000 Fuzzy Comparisons: Y<>Z but TEQ(Y,Z) is .TRUE. That is, Y is computationally equal to Z.

X=0.11D0, Y=X*11.D0-X-0.1D0, Z=1.D0 X= 0.11000000000000000 Y= 0.99999999999999989 Z= 1.0000000000000000
X=3FBC28F5C28F5C29 Y=3FEFFFFFFFFFFFFF Z=3FF0000000000000 Fuzzy FLOOR/CEIL: Y<>Z but TFLOOR(Y,EPS3)=TCEIL(Y,EPS3)=Z. That is, TFLOOR/TCEIL return exact whole numbers.

[certik]

Some feedback from the Fortran call:

• allow to set the tolerance in an argument
• NumPy has similar functionality, perhaps other languages also, we should explore prior art (https://numpy.org/doc/stable/reference/generated/numpy.allclose.html)
• make it work for arrays

[awvwgk]

What about an API similar to pytest.approx? This allows to create an object, which holds the tolerance and overloads the usual comparison operators in a convenient way.

[ghost]

It's called isClose in the D language.

Interface isClose(lhs, rhs, maxRelDiff, maxAbsDiff)

Documentation https://dlang.org/phobos/std_math.html#.isClose

Implementation https://github.com/dlang/phobos/blob/stable/std/math.d#L8847

[ghost]

GNU GSL also has this functionality in the function gsl_fcmp.

Interface int gsl_fcmp (const double x1, const double x2, const double epsilon);

Documentation https://www.gnu.org/software/gsl/doc/html/math.html#approximate-comparison-of-floating-point-numbers

Implementation (License: GPL) https://git.savannah.gnu.org/cgit/gsl.git/tree/sys/fcmp.c

[dalon-work]

This blog post discuss 4 different ways to compare floating-point numbers: https://codingnest.com/the-little-things-comparing-floating-point-numbers/