problem incrementing my variable

[Gary]

Here is my program, which does okay for awhile then messes up (output
on bottom):

#include <stdio.h>

main()
{
float p;

for ( p = 0; p <= 1; p = p + .005)
{
printf("p value is %f\n", p);
}

}

output
.....
p value is 0.695000
p value is 0.700000
p value is 0.705000
p value is 0.710000
p value is 0.714999
p value is 0.719999
.....

Thanks for insights.

 

[Ian Collins]

Here is my program, which does okay for awhile then messes up (output
on bottom):

#include <stdio.h>

main()
{
float p;

for ( p = 0; p <= 1; p = p + .005)
{
printf("p value is %f\n", p);
}

}

output
.....
p value is 0.695000
p value is 0.700000
p value is 0.705000
p value is 0.710000
p value is 0.714999
p value is 0.719999
.....

Thanks for insights.

Define messed up. Your output shows that 0.715 does not have an exact
float representation as a float. It may have one as a double, but what
are you trying to prove?

 

[Gary]

Gary wrote:

Define messed up.  Your output shows that 0.715 does not have an exact
float representation as a float.  It may have one as a double, but what
are you trying to prove?

Thanks double did work.

I am working with line segments defined by two ordered pairs (x1, y1)
and (x2, y2). I am trying to check if another ordered pair (x,y) is a
member of the line segment. I am using the p variable to increment
through the possible values on a resolution of .005 I don't know if
I said this right, because its been awhile since math, and I am trying
to automate some of my CAD work, so I am getting back into Math and
trying to learn C as well.

 

[dj3vande]

I am working with line segments defined by two ordered pairs (x1, y1)
and (x2, y2). I am trying to check if another ordered pair (x,y) is a
member of the line segment. I am using the p variable to increment
through the possible values on a resolution of .005 I don't know if
I said this right, because its been awhile since math, and I am trying
to automate some of my CAD work, so I am getting back into Math and
trying to learn C as well.

That's not the right way to do it.

It's a math problem, not a C problem, so I won't go into great depth
about it in comp.lang.c, but the summary is:
Every line (not line segment) is the set of points {(x0,y0)+t*(dx,dy)}
where (x0,y0) is point on the line, (dx,dy) is a direction vector, and
t is a real-valued parameter.
A line segment is a subset of these points with t constrained to be
inside a finite real interval.
This gives you a way to check directly whether a point (x,y) is on the
line, and it's a Rather Better way than trying to generate points on
the line and see if any of them match.

The details are an exercise in linear algebra.


ObC: <http://www.c-faq.com/fp/index.html>.


dave

 

[Kaz Kylheku]

Here is my program, which does okay for awhile then messes up (output
on bottom):

#include <stdio.h>

main()
{
float p;

You might want to use double for general-purpose floating calculations.

Note that most math functions in the C standard library pass and
return double.

Also note that when floating-point values of type float are passed as trailing
arguments to variadic functions, or as arguments to unprototyped functions,
they are promoted to double.

The type float has smaller limits on what its range and precision may be; it
may or may not be smaller than double, and thus is potentially useful as a
space-saving device in large arrays.

for ( p = 0; p <= 1; p = p + .005)

This is a very bad idea in general. The number .005 might not be exactly
representable in the floating point arithmetic your C programming system
is using (usually that of the underlying machine). So by doing these
repeated additions, you accumulate error.


You want:

double i;
const double step = 0.005;

for (i = 0; i < 1/step; i++) {
double p = i * step;

Here, i is still a floating-point number that is being stepped by repeated
additions. But small, positive integers are generally represented exactly in
floating point arithmetic. If you initialize an accumulator with with 0.0 and
keep adding 1.0, you don't expect errors until the number becomes too large to
be representable exactly in the mantissa of the floating-point format. Until
that point, the arithmetic will be exact.

With fractions, it's a different story, because some of them are not
representable all in a given-floating point format. The fraction 0.005, or
5/1000 or 1/200, cannot be represented in binary using a finite number of
digits.

0.005 in binary leads, I believe, to this expansion:
____________________
0.000000010100011110101110000

The digits with the bar over them repeat over and over, indefinitely,
exactly like the way the decimal expansion of 1/3 is 0.3333333...

Since your machine most likely uses a binary floating-point format, and not a
decimal representation, it simply truncates the above sequence of repeating,
and takes the truncated digits to be the mantissa. Truncation, of course,
produces a number which is no longer 0.005, only very close to it.

Now if you start adding these together, the error compounds. If 0.005 is not
really 0.005, then 0.005 + 0.005 is not really 0.01, and so on.

{
printf("p value is %f\n", p);

Here is an example of a float being converted to double
(trailing argument in variable-length argument list). That's why
it works; if the float was passed as float, it would be wrong, because the
printf conversion specifier %f extracts a value of type double.

}

}

output
....
p value is 0.695000
p value is 0.700000
p value is 0.705000
p value is 0.710000
p value is 0.714999

Ooops!

 

[Keith Thompson]

Here is my program, which does okay for awhile then messes up (output
on bottom):

#include <stdio.h>

main()
{
float p;

for ( p = 0; p <= 1; p = p + .005)
{
printf("p value is %f\n", p);
}

}

output
....
p value is 0.695000
p value is 0.700000
p value is 0.705000
p value is 0.710000
p value is 0.714999
p value is 0.719999
....

Thanks for insights.

Read section 14 of the comp.lang.c FAQ, <http://www.c-faq.com/>.

 

[Keith Thompson]

Define messed up. Your output shows that 0.715 does not have an exact
float representation as a float. It may have one as a double, but what
are you trying to prove?

0.715 almost certainly doesn't have an exact representation as float,
double, or long double. Typically only numbers that are integer
multiples of some power of two can be represented exactly. (This
might vary for some bizarre floating-point representations, and is
also subject to range and precision limits.)

 

[Kaz Kylheku]

0.715 almost certainly doesn't have an exact representation as float,
double, or long double.

The apparent 0.715 computed in the program most likely no longer even has the
representation that is obtained by writing the constant 0.715 in your program.
That 0.715 is the result of adding the approximation 0.005 to zero, 143 times.

 

[Keith Thompson]

Here is my program, which does okay for awhile then messes up (output
on bottom):

#include <stdio.h>

main()
{
float p;

for ( p = 0; p <= 1; p = p + .005)
{
printf("p value is %f\n", p);
}

}

output
....
p value is 0.695000
p value is 0.700000
p value is 0.705000
p value is 0.710000
p value is 0.714999
p value is 0.719999
....

Thanks for insights.

If you want to iterate over a range and expect to hit the end of the
range exactly, you need to use integers. (You can get away with using
floating-point in some restricted cases, with a sufficiently careful
choice of values, but we'll ignore that.)

Note that neither 0.715 nor 0.005 can be represented exactly. Since
the value you're adding to p each time isn't exactly 0.005, you have a
*cumulative* error. If the actual value is, say, 0.005-delta, then
for the value that you expect to be 0.715, the best you can hope for
is 0.715 - 143*delta.

Here's a version of your program the avoids the cumulative roundoff
error, by using an integer to control the iteration and computing each
floating-point value individually:

#include <stdio.h>
int main(void)
{
int i;

for (i = 0; i <= 200; i ++) {
double p = i / 200.0;
printf("%f\n", p);
}
return 0;
}

The output *looks* better, but you still have the problem that most of
the values of p still can't be represented exactly (the multiples of
0.125 can be). Change the "%f\n" format to, say, "%.50f\n", and
you'll see what I mean.

And using double, or even long double, rather than float can get you
closer to the exact value, but it can never get you all the way there.

 

[user923005]

Thanks double did work.

I am working with line segments defined by two ordered pairs (x1, y1)
and (x2, y2). I am trying to check if another ordered pair (x,y) is a
member of the line segment. I am using the p variable to increment
through the possible values on a resolution of .005    I don't know if
I said this right, because its been awhile since math, and I am trying
to automate some of my CAD work, so I am getting back into Math and
trying to learn C as well.

First, check that x is between x1 and x2
Second, check that y is between y1 and y2
If it fails either check, then the point is not on the segment.

The segment gives you the equation of a line.
Find the equation of the line.
Plug in the x component for the point into the line equation.
If the y value computed by the line equation is the y value you were
given (use fuzzy compare as recommended by the C-FAQ) then the point
is on the line.

 

[Bartc]

Thanks double did work.

I am working with line segments defined by two ordered pairs (x1, y1)
and (x2, y2). I am trying to check if another ordered pair (x,y) is a
member of the line segment.

You want to know if the point (x,y) is on the line (x1,y1) to (x2,y2)?

Whatever method you use (maybe comp.graphics.algorithms may be of help
here), and even with doubles, you may run into accuracy problems: how close
to the line or it's endpoints does (x,y) have to be, to be considered *on*
it?

If you imagine an XY plane with a fine grid corresponding to legal floating
point values (any X or Y only exists at an intersection), then a line drawn
from x1y1 to x2y2 quite likely will miss a lot of intersections on the way.
So you can't step smoothly along it, you can only follow the grid (and it
gets coarser away from 0,0).

I am using the p variable to increment
through the possible values on a resolution of .005

A double will still have problems but won't show up as easily. Things are
still workable, just remember floating point is rarely exact.

 

[CBFalconer]

Here is my program, which does okay for awhile then messes up
(output on bottom):

#include <stdio.h>
main() {
float p;

for ( p = 0; p <= 1; p = p + .005) {
printf("p value is %f\n", p);
}
}

output
....
p value is 0.695000
p value is 0.700000
p value is 0.705000
p value is 0.710000
p value is 0.714999
p value is 0.719999
....

0.005 cannot be expressed exactly as a float (or a double). You
can reduce the dependance by specifying that only 3 decimal places
be printed in the printf format string.

Note that printf expects to receive a double. You should so cast
the argument. This is one of the few places that a C cast is
needed.

 

[Ian Collins]

Note that printf expects to receive a double. You should so cast
the argument. This is one of the few places that a C cast is
needed.

Isn't that taken care of by default argument promotion?

 

[user923005]

You want to know if the point (x,y) is on the line (x1,y1) to (x2,y2)?

Whatever method you use (maybe comp.graphics.algorithms may be of help
here), and even with doubles, you may run into accuracy problems: how close
to the line or it's endpoints does (x,y) have to be, to be considered *on*
it?

The CGA FAQ won't help. This problem is considered too elementary for
a listing, I am sure.

If you imagine an XY plane with a fine grid corresponding to legal floating
point values (any X or Y only exists at an intersection), then a line drawn
from x1y1 to x2y2 quite likely will miss a lot of intersections on the way.
So you can't step smoothly along it, you can only follow the grid (and it
gets coarser away from 0,0).

It's the wrong way to solve it anyway.

A double will still have problems but won't show up as easily. Things are
still workable, just remember floating point is rarely exact.

Maybe something like this:

#include <float.h>
#include <math.h>
/* 'Close enough' comparison: */
int double_compare(double d1, double d2)
{
if (d1 > d2)
if ((d1 - d2) < fabs(d1 * DBL_EPSILON))
return 0;
else
return 1;
if (d1 < d2)
if ((d2 - d1) < fabs(d2 * DBL_EPSILON))
return 0;
else
return -1;
return 0;
}

/* Max between two doubles: */
double dmax(double a, double b)
{
if (a > b)
return a;
return b;
}

/* Min between two doubles: */
double dmin(double a, double b)
{
if (a < b)
return a;
return b;
}

/* Return 1 if point is on line segment, zero otherwise: */
int point_on_segment_2D(double x1, double y1, double x2, double y2,
double x, double y)
{
double x_1 = dmin(x1, x2);
double x_2 = dmax(x1, x2);
double y_1 = dmin(y1, y2);
double y_2 = dmax(y1, y2);
/* First, we check to see if the point is inside of the bounding
box: */
if ((x >= x_1) && (x <= x_2)) {
if ((y >= y_1) && (y <= y_2)) {
/* If the point is in the box, we check to see if it is on
the line: */
int result = double_compare(((x2 - x1) * (y - y1)), ((x -
x1) * (y2 - y1)));
return result == 0;
}
}
return 0;
}

#ifdef UNIT_TEST
#include <stdio.h>
int main(void)
{
double x1 = 0,
y1 = 0;
double x2 = 5,
y2 = 5;
double x,
y;
/* First, some points of which half are on the segment: */
for (x = -5, y = -5; x <= 5; x += 0.5, y += 0.5) {
if (point_on_segment_2D(x1, y1, x2, y2, x, y)) {
printf("(%g,%g) is on the segment\n", x, y);
}
}

/* Now, some points where half are close to, but NOT on the
segment: */
for (x = -5, y = -5; x <= 5; x += 0.5, y += 0.4999) {
if (point_on_segment_2D(x1, y1, x2, y2, x, y)) {
printf("(%g,%g) is on the segment\n", x, y);
}
}

return 0;
}
/* Possible output:
(0,0) is on the segment
(0.5,0.5) is on the segment
(1,1) is on the segment
(1.5,1.5) is on the segment
(2,2) is on the segment
(2.5,2.5) is on the segment
(3,3) is on the segment
(3.5,3.5) is on the segment
(4,4) is on the segment
(4.5,4.5) is on the segment
(5,5) is on the segment
*/
#endif

 

[Kaz Kylheku]

Isn't that taken care of by default argument promotion?

Yes, but that would require Chucky to read a good reference manual.

He's only been here a few months; still learning.

 

[Keith Thompson]

float p;

for ( p = 0; p <= 1; p = p + .005) {
printf("p value is %f\n", p);
}

[...]
0.005 cannot be expressed exactly as a float (or a double). You
can reduce the dependance by specifying that only 3 decimal places
be printed in the printf format string.

It reduces it, but doesn't eliminate it; the number of iterations is
still likely to vary, depending on whether the value of p on the 201st
iteration is just below, or just above, 1.0.

 

[Bartc]

The CGA FAQ won't help. This problem is considered too elementary for
a listing, I am sure.

I was thinking of a different approach, perhaps finding the point on the
line nearest x,y and deciding from the distance from x,y, whether it can be
on the line. The c.g.a group could be kept in mind by the OP for other
stuff,

Maybe something like this:

#include <float.h>
#include <math.h>
/* 'Close enough' comparison: */
int double_compare(double d1, double d2)
{
if (d1 > d2)
if ((d1 - d2) < fabs(d1 * DBL_EPSILON))

DBL_EPSILON is the smallest possible step for a floating point value? Within
the context of CAD, I would use a tolerance perhaps a million times bigger.
At least.

/* Return 1 if point is on line segment, zero otherwise: */ ....
/* First, we check to see if the point is inside of the bounding
box: */
if ((x >= x_1) && (x <= x_2)) {
if ((y >= y_1) && (y <= y_2)) {

I would use your 'fuzzy compare' for these too.

#ifdef UNIT_TEST
#include <stdio.h>
int main(void)
{
double x1 = 0,
y1 = 0;
double x2 = 5,
y2 = 5;
double x,
y;

And I would use a line that wasn't at 45 degrees for testing...

 

[Keith Thompson]

user923005 wrote: [...]

Maybe something like this:

#include <float.h>
#include <math.h>
/* 'Close enough' comparison: */
int double_compare(double d1, double d2)
{
if (d1 > d2)
if ((d1 - d2) < fabs(d1 * DBL_EPSILON))

DBL_EPSILON is the smallest possible step for a floating point value? Within
the context of CAD, I would use a tolerance perhaps a million times bigger.
At least.

DBL_EPSILON is "the difference between 1 and the least value greater
than 1 that is representable in the given floating point type" (double
in this case). The actual smallest step varies tremendously depending
on the magnitude of the value -- which is why it's called *floating*
point. DBL_EPSILON is far too coarse for values near zero, and
unrepresentably fine for large values.

I haven't used fuzzy floating-point equality comparisons very often,
but it seems to me that the problem is more difficult problem than the
way it's usually presented. You need to ask yourself just what it
means for two values to be "equal", and how the values you're trying
to compare were computed; the latter can affect how much of the
difference is real, and how much is due to rounding errors. Both of
these things can vary greatly depending on the application. It's not
just about comparing numbers; you have to understand what they *mean*.

 

[jameskuyper]

user923005 wrote:
....

#include <float.h>
#include <math.h>
/* 'Close enough' comparison: */
int double_compare(double d1, double d2)
{
if (d1 > d2)
if ((d1 - d2) < fabs(d1 * DBL_EPSILON))
return 0;
else
return 1;
if (d1 < d2)
if ((d2 - d1) < fabs(d2 * DBL_EPSILON))
return 0;
else
return -1;
return 0;
}

Note: double_compare(0.0, -DBL_MAX) returns 0.

You should be comparing fabs(d1) with fabs(d2) in order to determine
what scaling factor to use.

 

[Bartc]

user923005 wrote:

DBL_EPSILON is the smallest possible step for a floating point value?
Within
the context of CAD, I would use a tolerance perhaps a million times
bigger.
At least.

I missed the fact that DBL_EPSILON is scaled before use.

Even so, using a value a million times bigger means that on a line 1000km
long (from the origin), a point needs to be within 0.2mm to be considered on
the line.