# Function for comparing floating point numbers

I am trying to make a function for floating point numbers comparison.
The goal is to "evaluate" all operators in single function: >, >=, =, <=, <. If X>Y, then F>0. If X==Y, F==0. If X<=Y, F<=0 etc (where F is returned value).
I took a piece of code from numeric limits. First thing I changed, in scaling of epsilon, fabs(x+y) to fabs(x)+fabs(y) (because if x=-y, we would be comparing diff to 0; which probably would work fine, but sounds kinda stupid). I've also seen an approach using max(fabs(x), fabs(y)), can we say that one of them is better than the others?
Third parameter, DecadesOfInAccuracy, means how many digits above epsilon() should be ignored.

#include <limits>
#include <iostream>

template<class T>
typename std::enable_if<!std::numeric_limits<T>::is_integer, int>::type
cmp(T x, T y, unsigned int doia = 0)
{
std::cout << "X: " << x << std::endl;
std::cout << "Y: " << y << std::endl;
std::cout << "D: " << fabs(x - y) << std::endl;
std::cout << "M: " << std::numeric_limits<double>::min() * pow(10, doia+1) << std::endl;
std::cout << "E: " << std::numeric_limits<T>::epsilon() * (fabs(x) + fabs(y)) * pow(10, doia + 1) << std::endl;

if (fabs(x - y) < std::numeric_limits<double>::min() * pow(10, doia + 1))
return 0;
bool ltoeq = (x - y) <= std::numeric_limits<T>::epsilon() * (fabs(x)+fabs(y)) * pow(10, doia + 1);
bool gtoeq = (y - x) <= std::numeric_limits<T>::epsilon() * (fabs(x) + fabs(y)) * pow(10, doia + 1);
std::cout << "L: " << ltoeq << std::endl;
std::cout << "G: " << gtoeq << std::endl;
return gtoeq - ltoeq;
}

int main()
{
std::cout.precision(std::numeric_limits<double>::max_digits10 + 5);
double x = 0.0001000000000000001;
double y = 0.0001;

int res = cmp(x, y, 0);
std::cout << res << std::endl << std::endl;
if (res > 0)
std::cout << "X>Y" << std::endl;
if (res >= 0)
std::cout << "X>=Y" << std::endl;
if (res == 0)
std::cout << "X==Y" << std::endl;
if (res <= 0)
std::cout << "X<=Y" << std::endl;
if (res < 0)
std::cout << "X<Y" << std::endl;
return 0;
}


It "correctly" returns equality for pairs like 0.0001000000000000001; 0.0001, 1.0-4*0.2; 0.2, 0; -0, 100000000000.0001; 100000000000.
Is there something I forgot to consider, or maybe something can be simplified? (of course except couts, which will be removed :P).
I am using VS2019 if that matters.

• Welcome to code review where we review working code to provide suggestions on how the code can be improved. What is the expected output? Is this code working as expected? Exactly how many relations to expect to be output at the end. If this code is not working as expected then it is off-topic for the code review site. Dec 26, 2019 at 14:20
• It works as expected for these examples. I need a way to test much more cases or someone to evaluate if the code will always work. Expected output is on the beginning. Positive if x is definitely smaller than y, 0 for equal and negative if x is definitely greater than y. Dec 26, 2019 at 14:30
• Sorry for being obvious but to add more test cases make a struct with test data and then have an array or vector of those structs. Dec 26, 2019 at 14:36
• @pacmaninbw I did some testing, function works perfectly in range e308 => ~e-294 (for max precision). You can see the output here: github.com/herhor67/C-C-/tree/master/float_comp Unfortunately it can't compare to 0.0, which is most desired for me. Maybe I should add a special case? Dec 26, 2019 at 21:56
• Now it seems that the thread would fit better at stack overflow... Dec 26, 2019 at 23:24

Is there something I forgot to consider

Overflow

fabs(x) + fabs(y) is prone to overflow, even if mathematically epsilon() * (fabs(x) + fabs(y)) * pow(10, doia + 1) is representable as a T, resulting in incorrect results.

I'd expect cmp(T x, T y, unsigned int doia = 0) to work over the entire range of x,y.

// alternative
(fabs(x)/2 + fabs(y)/2) * (epsilon() * pow(10, doia + 1) * 2)


Not-a-number

The generation of gtoeq, ltoeq are 0 when either/both x or y are a not-a-number. Thus code returns 0, incorrectly implying near equality. Perhaps return something to indicate more that 3 conditions like a 4-bit (==,>,<,not comparable) or a double -1,0,1,NAN.

Precision

The + 1 in doia + 1 does not well allow code to just use epsilon, but minimally must use 10 * epsilon. Consider dropping the + 1.

Rather than an integer indicating some power-of-ten * epsilon, consider simply using the desired precision or n * epsilon(). For me I'd rather use a power-of-2.