Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Next revision
Previous revision
floating-point-equality [2015/03/10 01:54]
awf created
floating-point-equality [2018/09/05 15:14] (current)
awf
Line 1: Line 1:
-==== Is my floating point number equal to zero? ====+===== Is my floating point number equal to zero? =====
  
 One often sees code like the following, embodying the feeling programmers have that "you should never compare a float to zero". One often sees code like the following, embodying the feeling programmers have that "you should never compare a float to zero".
Line 13: Line 13:
 </​code>​ </​code>​
  
-But what are you trying to avoid here?   Is it the divide by zero exception?  Then avoid that:+But what are you trying to avoid here?   Is it the divide by zero?  Then avoid that...
 <code C++> <code C++>
 double func(double a, double b) { double func(double a, double b) {
Line 23: Line 23:
 } }
 </​code>​ </​code>​
-"​But",​ you say, "what if diff is really small"? ​ You mean like 1e-128? ​ Then return 1e+128. ​ Floating point can deal with it.   If you didn't need that dynamic range, you'ld be using floats.  ​And they would still deal fine with it  If you're worried about floating underflow, then compare ​to a constant ​you named MAX_OF_FLT_MIN_AND_RECIPROCAL_OF_FLT_MAX.+"​But",​ you say, "what if diff is really small"? ​ You mean like 1e-128? ​ Then this function will return 1e+128.  ​ 
 +Floating point can deal with it.   If you didn't need that dynamic range, you'ld be using ''​float''​ instead of ''​double''​.  ​ 
 + 
 +OK, what about 1e-309? ​ Aha, different issue (denormal floats)  
 +Now you won't get a divide by zeromost likely you'll return Inf, and... your code might run very slowly. 
 +Aha! No-one wants slow code, right? ​ OK, but then why the "​if"​ in the first place? ​  
 +That's slowing down every call to the function. ​  
 +If your workloads mostly don't call the function with nearly equal very small numbers,  
 +then you can go lot faster in the common case by not checking. 
 + 
 +"Ah, but I want the additional error checking"​. ​ OK, fine.  I'm with you.  But is this the right place? ​  
 +Why did the caller call the routine with nearly equal very small numbers? ​  
 +Is that physically likely in whatever real-world problem you're trying to solve? ​  
 +If not, one should probably protect at the call site with an assert.
  
 Or maybe it's a precondition of the function call that a not be equal to b, and it's a coding bug if not satisfied. Or maybe it's a precondition of the function call that a not be equal to b, and it's a coding bug if not satisfied.
 <code C++> <code C++>
 double func(double a, double b) { double func(double a, double b) {
-   ​assert(a != b);+   ​assert(a ​- b != 0);
        
    ​return 1/(a-b);    ​return 1/(a-b);
 } }
 </​code>​ </​code>​
-Don't forget to [[assert-always|leave your asserts on in all builds]] unless the profiler tells you otherwise.+Don't forget to [[assert-always|leave your asserts on in all builds]] unless the profiler tells you otherwise.   
 +But hang on, you're on a platform which will give you line number here when you get the divby0 right? ​   
 +So the assert is superfluous -- just enable exceptions on divide by zero, underflow, and overflow.
  
-=== But my case is more subtle than that ===+And there'​s also a very simple case. If I stored 0.3f somewhere, I expect to get it back.  It's irrelevant that 0.3f doesn'​t have an exact binary representation -- whatever it translates to, I should get exactly that back 
 +<code C++> 
 +void my_unit_test() { 
 +  my_matrix<​2,​2,​float>​ a = {3.0f, 1.1f, 2.2f, 4.0f}; 
 +  TEST_ASSERT_EQ_APPROX(a(0,​0),​ 3.0f, 1e-8f); // No, no, no.  There is no reason for the stored value to be different. 
 +  TEST_ASSERT_EQ(a(0,​0),​ 3.0f); // I gave you whatever 3.0f is.  I want it back. 
 +
 +</​code>​ 
 +You may worry about excess precision on x86, but there'​s no computation going on here, so we can and should ask for exact equality. 
 + 
 +==== But my case is more subtle than that ====
  
 <code C++> <code C++>
Line 46: Line 71:
 </​code>​ </​code>​
  
-Even here, you can think about what a better value to compare to might be.   First note that sin(x)/x works fine for all normalized ​floats except zero, so you can avoid a dependency on fabs and a definition of epsilon by just comparing to zero.+Even here, you can think about what a better value to compare to might be.   First note that sin(x)/x works fine for all floats except zero, so you can avoid a dependency on fabs and a definition of epsilon by just comparing to zero.
  
 But perhaps the profiler has told you that a big chunk of program time is spent in computing sin() and dividing, and furthermore you know that most times it's called with values near zero.  Well then, let's do what we do with any special function that goes slow: chop up the domain and special case.   ​It'​s true that sin(x)/x is exactly floating point 1.0 when x is small. ​ In fact, compute the taylor series, see that But perhaps the profiler has told you that a big chunk of program time is spent in computing sin() and dividing, and furthermore you know that most times it's called with values near zero.  Well then, let's do what we do with any special function that goes slow: chop up the domain and special case.   ​It'​s true that sin(x)/x is exactly floating point 1.0 when x is small. ​ In fact, compute the taylor series, see that
Line 62: Line 87:
 } }
 </​code>​ </​code>​
 +
 +It's a good idea then to write a blog post about the issue, because Toby Sharp might comment (see below) and point out that there is a much better choice:
 +
 +> "... since we need some kind of test for x=0 anyway, I don't see why we shouldn'​t take the opportunity to improve accuracy in this vicinity. And yes, it only costs 2 multiplies which is less than using the sin in those cases. So we are improving both performance and accuracy in the range (0, threshold). threshold is chosen to be the largest value such that the quadratic approximation is as good as the direct evaluation."​
 +
 +<code C++>
 +float sinc(float x) { 
 +  // Toby's threshold, determined through binary search as the best value for minimizing absolute error
 +  if (std::​abs(x) < 0.0406015441f)
 +      return 1.0f - (x * x) * (1.0f / 6.0f);
 +  else
 +      return std::sin(x) / x;
 +
 +</​code>​
 +
  
 === And templates? === === And templates? ===
Line 70: Line 110:
 Also, compares to zero may be faster on your platform. ​  And you don't have to define epsilon somewhere. Also, compares to zero may be faster on your platform. ​  And you don't have to define epsilon somewhere.
  
 +==== And more... ====
 +
 +
 +<code C++>
 +double sincminusone(double x) {
 +   if (fabs(x) < 1e-6)
 +     ​return 0; // But x^2/6 gives the right answer, not just truncating to the nearest zero.
 +   else
 +     ​return sin(x)/x-1; // Or try (sin(x)-x)/​x and compare to zero instead...
 +}
 +</​code>​
 +
 +See also
 +* https://​stackoverflow.com/​questions/​51134021/​floating-point-equality
 +
 +~~DISCUSSION~~