0

I am writing a numerical code that needs to make extensive (and possibly fast) comparisons among double precision numbers. My solution to compare two numbers A and B involves shifting A to the left (or right) by an epsilon and checking whether the result is bigger (or smaller) than B. If so, the two doubles are the same. (Extra coding needs to be done for negative or zero numbers).

This is the comparing function:

#define S_
inline double s_l (double x){
    if(x>0){return 0.999999999*x;}
    else if(x<0){return 1.00000001*x;}
    else {return x-0.000000000001;}
}
inline double s_r (double x){
    if(x>0){return 1.00000001*x;}
    else if(x<0){return 0.999999999*x;}
    else{return x+0.000000000001;}
}
inline bool s_equal (double x,double y){
    if(x==y){return true;}
    else if(x<y && s_r(x)>y){return true;}
    else if(x>y && s_l(x)<y){return true;}
    else{return false;}
}
#endif

Since this is part of a MonteCarlo algorithm and s_equal(x,y) is called millions of times, I wonder if there is any better or faster way to code this, understandable at a simple level.

Bill the Lizard
  • 369,957
  • 201
  • 546
  • 842
Ale
  • 5
  • 1
  • 2
    I do something like abs( (x-y)/x ) < 1.0e-10 – brian beuning Apr 23 '13 at 23:33
  • "Almost equal" is an advanced technique and should not be lightly undertaken. For example if `a` almost equals `b` and `b` almost equals `c`, it does not follow that `a` almost equals `c`. This can lead to unanticipated complications. – Pete Becker Apr 24 '13 at 00:39
  • @brianbeuning Consider posting that as an answer. It's a *very* close match to the presented algorithm, and by removing all that if-branching it should be much faster than what's presented. – Drew Dormann Apr 24 '13 at 03:30

3 Answers3

0

If you're using the C++11, then you could use the new math library functions, such as:

bool isgreater(float x, float y)

More documentation on std::isgreater can be had here.

Otherwise, there's always is_equal in boost. Also, SO already has a bunch of related (not sure if same) questions such as the ones here, here and here.

Community
  • 1
  • 1
Bhargav Bhat
  • 8,939
  • 1
  • 17
  • 29
  • Thanks. I am not sure this solves my problem though. I tried to find a definition of isgreater or is_equal, and cannot see if they do epsilon comparison and if so, to which precision. Regarding other threads, I did search for them, but found out that most of the discussion is about "fixed" epsilon comparison with abs(x-y) – Ale Apr 24 '13 at 20:02
  • 1
    @Ale http://www.boost.org/doc/libs/1_36_0/libs/test/doc/html/utf/testing-tools/floating_point_comparison.html – Jeff Aug 06 '14 at 12:52
0

I do something like abs( (x-y)/x ) < 1.0e-10 .

You need to divide by x in case both values are huge or tiny.

brian beuning
  • 2,664
  • 13
  • 18
  • Thanks for the reply! I thought about this kind of implementation, but realized it can be dangerous: 1) when x is 0, it will give NaN 2) when y is 0, the lhs is equal to 1, no matter how x is close to 0. Including 'abs' in the solution seemed natural to me as well in order to get rid of many if-statements, but treating the numbers near 0 didn't seem so obvious. Did you think about this problem? How did you solve it? – Ale Apr 24 '13 at 19:42
  • abs( ((x == 0.0) ? y : ((y == 0.0) ? x : (x-y)/x)) ) < 1.0e-10 – brian beuning Apr 24 '13 at 23:35
0

I was surprised to find a significant speedup by avoiding all the double-precision math:

#define S_L(x) (x)+((x)<0?1024:-1024)
#define S_R(x) (x)+((x)<0?-1024:1024)
#define S_EQUAL(x,y) (S_L(x)<(y) && S_R(x)>(y))

double foo;                                                                      
double bar;                                                                      
long *pfoo;                                                                      
long *pbar;                                                                      

pfoo = (long*)&foo;                                                              
pbar = (long*)&bar;          

double result1 = S_R(*pfoo);
double result2 = S_L(*pbar);
bool result3 = S_EQUAL(*pfoo, *pbar);

(In testing, I operated on randomly-generated doubles between -1M and 1M, executing each operation 100M times with different input for each iteration. Each operation was timed in an independent loop, comparing system times - not wall times. Including loop overhead and generation of random numbers, this solution was about 25% faster.)

A word of warning: there are lots of dependencies here on your hardware, the range of your doubles, the behavior of your optimizer, etc., etc. Such pitfalls are commonplace when you start second-guessing your compiler. I was shocked to see how much faster this was for me, since I'd always been told that integer and floating point units are kept so separate on hardware that the transport of bits from one to the other always requires a hardware memory operation. Who knows how well this will work for you.

You will likely have to play with the magic numbers a bit (the 1024s) to get it to do about what you want it to - if it's even possible.

Sniggerfardimungus
  • 10,708
  • 10
  • 46
  • 91
  • Thanks. It is in an interesting suggestion, but scaring. The code has to be portable and compiled by anyone in the community. I'd rather sacrifice some speed to make something less prone to errors! – Ale Apr 24 '13 at 20:30