2

Hello I have some difficulties with the ceil function in c++ :

I've got a regular grid of points, and I need to perform an interpolation over it to compute the z-values of a set of points.

In order to do that, I need, for each computed point, to get the nearest points on the grid. I do it like this :

y1 = dy*floor(p.y/dy);
y2 = dy*ceil(p.y/dy);

where dy is the space beetween two points of the grid. (y1, p.y and dy are double) If I display the results using 

cout << static_cast<double>(p.y/dy) << ": " << y1 << ", " << y2 << endl;

I've got these strange results :

0: 0, 0
1: 0.1, 0.1
2: 0.2, 0.2
3: 0.3, 0.4

The three first results are OK but the last one is wrong and make an assertion fail.

I would like to know where did this strange error came and how to avoid it. Thanks.

I apologize for my English

EDIT

I call the function with dy = 0.1, but during the execution, it take the folowing value dy = 0.10000000000000001. p.y is initialized like that : 

 const uint N = round((x2 - x1) / dx2);
 const uint M = round((y2 - y1) / dy2);

 double p = persistence;
 double n = number_of_octaves;

 // generation of the points where the perlin noise is generated
 std::vector<Vertex3d> ret;
 std::vector<Vertex3d> dummy;
 for (uint i=0;i<=N;++i)
 {
        for (uint j=0;j<=M;++j)
        {
                ret.push_back({.x = i*dx2, .y=j*dy2, .z=0});
                dummy.push_back({.x = i*dx2, .y=j*dy2, .z=0});
        }
 }

where x1 = 0 and x2 = 1 (according gdb)

Pascal Cuoq
  • 79,187
  • 7
  • 161
  • 281
Alexandre Hoffmann
  • 341
  • 2
  • 9
  • .. And give us the values used in the calculations (i.e. p.y and dy) – Ed Heal Mar 08 '13 at 10:25
  • In addition to the excellent answers given below, you might want to see this post on floating point comparison as I'm sure you've got other bugs lurking about in your code. For example, your assert is probably expecting a very specific value. Don't do that! http://stackoverflow.com/questions/17333/most-effective-way-for-float-and-double-comparison – CadentOrange Mar 08 '13 at 10:27
  • When you mean the actual code, you mean the 62 lines of the function, the 48 lines of the function that call this function and the 9 of the main that call the function that call the function ? There is nothing beetween the 3 lines that I showed you, and all of these values are parameters. – Alexandre Hoffmann Mar 08 '13 at 11:02
  • In one of the comments, you state that `p.y` is 0.30000000000000004 and `dy` is 0.10000000000000001. Obviously `p.y` is more than 3 times `dy`, so `p.y/dy` is greater than three, so `ceil(p.y/dy)` is 4, so `dy*ceil(p.y/dy)` is (about) .4. You say “the last one is wrong”, but clearly all the math has been done correctly, and .4 is the nearest grid point above `p.y`. So your problem description is inadequate; it does not explain why you do not want .4 for this grid point. – Eric Postpischil Mar 08 '13 at 11:20
  • In another comment, you ask whether there is a way to assure that the value of `dy` is .1 rather than 0.10000000000000001. No, this is not possible as long as you use binary floating point. The value .1 is not exactly representable in binary floating point. One way to reduce such problems is to calculate in multiples of values that are exactly representable, such as multiples of 1 or of .125 or .0078125 other powers of two. That is, you might scale your points to match the binary values. However, since you have not shown the other arithmetic you are doing, there could be other problems. – Eric Postpischil Mar 08 '13 at 11:22
  • Your edit to the problem states “p.y is initialized like that:”, but the new code does not contain either “p.y” or “p”, so it does not initialize `p.y`. Your code does appear to be doing something with a member named `y` in some objects related to `ret` and `dummy`. Perhaps this is the value for `.y` you are referring to? – Eric Postpischil Mar 08 '13 at 11:33
  • Even if `p.y` is initialized with the expression you show for `.y` in `.y=j*dy2`, you have not shown the value of `dy2`. We can deduce `dy2` is not the double that results from “.1”, because the C or C++ source code `3*.1` yields the result 0.299999999999999988897769753748434595763683319091796875, not the 0.30000000000000004 that you have reported. If you want help, you **must** show exactly the problem. You must show **exactly and completely** how `p.y` is produced, including all expressions used to produce it and all values used in those expressions. – Eric Postpischil Mar 08 '13 at 11:38

3 Answers3

1

When you are using floating-point arithmetic, you generally should not rely on the result to be exactly an integer, because of imprecision in floating point calculations. You probably need to re-design your code so that it doesn't rely on this.

Here is a classic discussion of floating point arithmetic.

  • @EricPostpischil, you are right. My language was not precise enough. –  Mar 08 '13 at 11:38
1

Given your results I suppose dy is equal to 0.1.

Your results aren't wrong. If you have 0.3 < p.y < 0.4 then you will have 3 < p.y/dy < 4 so the ceil will be 4 and the floor will be 3.

Maybe you are confused because somewhere else in your code you set p.y to 0.3 or 0.4. You should know that floats are not that precise. It means that even if you set p.y = 0.3 it can have a value of 0.30000001 or something like that, causing your issue.

alestanis
  • 21,519
  • 4
  • 48
  • 67
  • in fact, dy = 0.10000000000000001, and p.y = 0.30000000000000004, but when I call my function, I give it dy = 0.1. Is there a way to assure that the value of dy is realy 0.1 ? – Alexandre Hoffmann Mar 08 '13 at 10:58
  • @AlexandreHoffmann Not with binary floating point numbers, since they cannot represent `0.1` exactly. – Christian Rau Mar 08 '13 at 11:43
-1

This looks simply like a case of rounding errors for floating point arithmetic.

It could be that the result in the last case is 3.00000000000001 so its ceiling is 4 not 3 (and then is multiplied by 0.1)

CashCow
  • 30,981
  • 5
  • 61
  • 92
  • This is not a rounding error in division. In IEEE-754 floating-point arithmetic, division of one floating-point number by another floating-point number will never produce a value greater than three unless the exact quotient is greater than three. – Eric Postpischil Mar 08 '13 at 11:25
  • I never said the rounding error happened during the division, it was already there and then got exposed when he did ceil() – CashCow Apr 11 '13 at 09:58
  • It is difficult to reconcile the statement “I never said the rounding error happened during the division” with the edit history that shows the original statement “This looks simply like a case of rounding errors for floating point division.” – Eric Postpischil Apr 11 '13 at 12:33