-1

I have an array of "x" values (a grid for a PDE solver) that when I pass to a function that fills another array based on these x-values, a particular expression involving one x-value does not evaluate properly. The x-values range is -1:1 at an increment of 0.0125, and at x = -0.5, and at x = 0.5, I need to handle these cases differently than the other values. However the block below fails to evaluate to TRUE for the point x = 0.5 (it's ok for x = -0.5). Here's a stripped down snippet of the problem block, with details to follow:

int N = 160;
double delta_x = 0.0125;
const double lims = 0.5 * delta_x;

for(int i = 0; i <= N; i++)
{                                
  if((x[i] < -0.5) || (x[i] > 0.5)) sol[i] = 0;
  else if( (abs(x[i] + 0.5) < lims) || (abs(x[i] - 0.5) < lims) )  sol[i] = 0.5;
  else sol[i] = 1;

  cout << setprecision(30) << "lims: " << lims << ", abs(x[i] - 0.5): " << abs(x[i] - 0.5) << endl;
  cout << "sol[" << i << "]: " << sol[i] << endl;                                                                                         
}

Here is the output for x = 0.5:

lims: 0.00625000000000000034694469519536, abs(x[i] - 0.5): 1.11022302462515654042363166809e-16
sol[120]: 0

So, it's looks like the expression in the if-statement should return TRUE when x = 0.5, even though it's not precisely 0.5 of course, since it is within the "lims" of the range. Any thoughts??

bcf
  • 2,104
  • 1
  • 24
  • 43
  • 1
    That sure looks like an off-by-one bug to me. Surely it should be < N instead of <= N? – Hans Passant Nov 08 '13 at 22:36
  • 1
    Just checking -- are you including `cmath.h`? `abs` is only overloaded for doubles in that, `stdlib` declares int-only `abs`. – Jongware Nov 08 '13 at 22:37
  • 1
    Please clarify what's not working as you expect. In the example output, 1e-16 is certainly within the tolerance `lims`. – MooseBoys Nov 08 '13 at 22:37
  • @HansPassant, The variable "N" is the number of x-values starting from 0, so the range of x-values are x_0:x_N for N+1 total array elements. Thanks – bcf Nov 08 '13 at 22:38
  • I don't have that many neurons left, but one of the five is permanantly wired to flash madly whenever I see "<=" instead of just "<" in a C/C++/Java for loop. Not that such code must be wrong, just likely that it is. – DarenW Nov 08 '13 at 22:39
  • @David: I don't understand. The output you posted clearly shows that the `if` **is** entered for `x = 0.5`. What made you think that it is not entered? On the contrary, I don't see any output for `x = -0.5`. Yet you claim that `x = -0.5` "works". So what is all that supposed to mean? – AnT stands with Russia Nov 08 '13 at 22:39
  • @Jongware My include line is #include -- Correct me if I'm wrong but I think that's the same as #include "cmath.h"? – bcf Nov 08 '13 at 22:41
  • The cause of this is most definitely that `abs(int)` is called. Just replace `abs` with `fabs`. – Paul Groke Nov 08 '13 at 22:42
  • @Paul Groke: If `abs(int)` was called, then how did he manage to obtain `1.11022302462515654042363166809e-16` output after sending the result of `abs` to `cout`??? – AnT stands with Russia Nov 08 '13 at 22:43
  • #include is just the C++ header of #include , so you are doing that right. Are you linking it as wel? – Arjen Nov 08 '13 at 22:43
  • I have a strong suspicion that the OP either missed the `e-16` suffix in the output or simply doesn't know what that `e-16` means. I see no other explanation for his claim that `x = 0.5` skips the `if`, when the output clearly shows that `if` is entered. – AnT stands with Russia Nov 08 '13 at 22:47
  • @AndreyT: Good question. In don't know. :-) I just know that his code works with `fabs`, and there is no "off by one" error. See http://ideone.com/8aMIUO – Paul Groke Nov 08 '13 at 22:59
  • Note that in both `if` statements you could remove two-thirds of the parentheses without changing the meaning. – Pete Becker Nov 09 '13 at 12:26

2 Answers2

3

Answer to the original question

The expression in the if statement did evaluate true. If it did not, then you would not have seen the output that you included in the question.

FWIW, a simpler test would be for

abs(abs(x[i]) - 0.5) < lims

Answer to the latest version of the question

if((x[i] < -0.5) || (x[i] > 0.5)) sol[i] = 0;
else if( (abs(x[i] + 0.5) < lims) || (abs(x[i] - 0.5) < lims) )  sol[i] = 0.5;
else sol[i] = 1;

You state that x[i] is close to 0.5, but does not set sol[i] to 0.5 and in fact sets it to 0. In which case the only sane conclusion is that x[i] > 0.5 and the first condition is is met:

if((x[i] < -0.5) || (x[i] > 0.5)) sol[i] = 0;

So you need to change the order of your tests:

if( (abs(x[i] + 0.5) < lims) || (abs(x[i] - 0.5) < lims) )  sol[i] = 0.5;
else if((x[i] < -0.5) || (x[i] > 0.5)) sol[i] = 0;
else sol[i] = 1;

And I would write it like this:

if (abs(abs(x[i]) - 0.5) < lims)
    sol[i] = 0.5;
else if ((x[i] < -0.5) || (x[i] > 0.5)) 
    sol[i] = 0;
else 
    sol[i] = 1;

Please in future be sure to include the right code in the question that you ask.

David Heffernan
  • 601,492
  • 42
  • 1,072
  • 1,490
  • 1
    Exactly. The output quoted in the question was generated for `x = 0.5`. Yet the OP claims that for `x = 0.5` the `if` is not entered. The question, as stated, is contradicting itself. – AnT stands with Russia Nov 08 '13 at 22:45
  • @AndreyT Thanks for noticing that, my mistake. I was trying to make it simpler but then it wasn't showing the error. I just edited it so that the error results. – bcf Nov 08 '13 at 22:49
  • @David Heffernon - Thanks! I just realized I need to simply check that problem if-statement first. – bcf Nov 08 '13 at 22:56
  • Leaving out the headers for brevity added a layer to the confusion as well. Mistaking `abs` for `fabs` is, alas, a common mistake -- and, as I just tested, `-Wall` on `printf("%f", abs(1.5));` only warns for the `%f`, not for the implicit conversion to int. – Jongware Nov 08 '13 at 22:59
1

Well, now after you corrected your code, the problem is obvious. For x[i] = 0.5 the control is intercepted by the very first condition

if((x[i] < -0.5) || (x[i] > 0.5)) sol[i] = 0;

It doesn't even get to your abs evaluations.

Floating-point number are not always precise. While it is true that 0.5 can be represented precisely, it only applies to situations when you assign 0.5 directly or use some simple and particularly stable evaluations to arrive at that value (like in case of 1.0/2). In more complex cases you might end up with imprecise 0.5 value (like in case of 0.1 * 5). And that is what's causing it to get intercepted by that first if. Your positive 0.5 is imprecise and just happens to be greater than the precise 0.5.

One way to fix it might be to put your approximate lim-based comparisons first and make "precise" comparisons follow 

if( (abs(x[i] + 0.5) < lims) || (abs(x[i] - 0.5) < lims) )  sol[i] = 0.5; 
else if((x[i] < -0.5) || (x[i] > 0.5)) sol[i] = 0;
else sol[i] = 1;
AnT stands with Russia
  • 312,472
  • 42
  • 525
  • 765