floating point problems

[G Patel]

Hello,

I'm having trouble with floating point problems.

I'm trying to write a function that calculates the distance between two
cartesian points (integer coordinates). I have the function working,
but it fails on one test case.

Here is my function:

double
distance(int x1, int y1, int x2, int y2)
{
double deltax = (double)x1 - x2;
double deltay = (double)y1 - y2;
return sqrt(deltax*del tax + deltay*deltay);
}


It works fine for most cases except that if deltax is really small
compared to deltay (or vice versa), the addition inside the sqrt
argument causes the smaller square value to be lost (because of adding
a large float number with a small float number)).

example: deltax = 1 deltay = 1999998

What this distance function returns is exactly what would be returned
if only one of the dimensions were 1999998.

This is crucial because I need the distance formula to work within a
certain tolerance so I can check to see if triangles are right or
oblique.

Does anyone know any tricks that I can use to preserve the "effect" of
the small deltax in the distance?

Thanks

[P.J. Plauger]

Re: floating point problems

"G Patel" <gaya.patel@gma il.com> wrote in message
news:1138912880 .138182.67920@g 49g2000cwa.goog legroups.com...
[color=blue]
> I'm having trouble with floating point problems.
>
> I'm trying to write a function that calculates the distance between two
> cartesian points (integer coordinates). I have the function working,
> but it fails on one test case.
>
> Here is my function:
>
> double
> distance(int x1, int y1, int x2, int y2)
> {
> double deltax = (double)x1 - x2;
> double deltay = (double)y1 - y2;
> return sqrt(deltax*del tax + deltay*deltay);
> }
>
>
> It works fine for most cases except that if deltax is really small
> compared to deltay (or vice versa), the addition inside the sqrt
> argument causes the smaller square value to be lost (because of adding
> a large float number with a small float number)).
>
> example: deltax = 1 deltay = 1999998
>
> What this distance function returns is exactly what would be returned
> if only one of the dimensions were 1999998.
>
> This is crucial because I need the distance formula to work within a
> certain tolerance so I can check to see if triangles are right or
> oblique.
>
> Does anyone know any tricks that I can use to preserve the "effect" of
> the small deltax in the distance?[/color]

It ain't easy. What you need is a well written hypot function
from C99. Barring that, take a look at the (template) function
_Fabs in the header xcomplex in our Standard C++ library. It
computes the magnitude of a complex number, which is the same
problem as you have. Reading it is tough sledding, but it's
the minimum you have to do to preserve all those bits you need.

Alternatively, I suspect there's some way to recast your problem
to an atan2 call, if all you really care about is right-angle
detection.

HTH,

P.J. Plauger
Dinkumware, Ltd.

http://www.dinkumware.com

[Rod Pemberton]

Re: floating point problems


"P.J. Plauger" <pjp@dinkumware .com> wrote in message
news:pJmdnQBf54 mo7n_enZ2dnUVZ_ sSdnZ2d@giganew s.com...[color=blue]
> "G Patel" <gaya.patel@gma il.com> wrote in message
> news:1138912880 .138182.67920@g 49g2000cwa.goog legroups.com...
>[color=green]
> > I'm having trouble with floating point problems.
> >
> > I'm trying to write a function that calculates the distance between two
> > cartesian points (integer coordinates). I have the function working,
> > but it fails on one test case.
> >
> > Here is my function:
> >
> > double
> > distance(int x1, int y1, int x2, int y2)
> > {
> > double deltax = (double)x1 - x2;
> > double deltay = (double)y1 - y2;
> > return sqrt(deltax*del tax + deltay*deltay);
> > }
> >
> >
> > It works fine for most cases except that if deltax is really small
> > compared to deltay (or vice versa), the addition inside the sqrt
> > argument causes the smaller square value to be lost (because of adding
> > a large float number with a small float number)).
> >
> > example: deltax = 1 deltay = 1999998
> >
> > What this distance function returns is exactly what would be returned
> > if only one of the dimensions were 1999998.
> >
> > This is crucial because I need the distance formula to work within a
> > certain tolerance so I can check to see if triangles are right or
> > oblique.
> >
> > Does anyone know any tricks that I can use to preserve the "effect" of
> > the small deltax in the distance?[/color]
>
> It ain't easy. What you need is a well written hypot function
> from C99.[/color]

I'm not sure if it's C99, but there is a hypot function in Sun Microsystem's
Freely Distributable LIBM:

http://www.netlib.org/fdlibm/

Rod Pemberton

[G Patel]

Re: floating point problems


Rod Pemberton wrote:[color=blue]
> "P.J. Plauger" <pjp@dinkumware .com> wrote in message
> news:pJmdnQBf54 mo7n_enZ2dnUVZ_ sSdnZ2d@giganew s.com...[color=green]
> > "G Patel" <gaya.patel@gma il.com> wrote in message
> > news:1138912880 .138182.67920@g 49g2000cwa.goog legroups.com...
> >[color=darkred]
> > > I'm having trouble with floating point problems.
> > >
> > > I'm trying to write a function that calculates the distance between two
> > > cartesian points (integer coordinates). I have the function working,
> > > but it fails on one test case.
> > >
> > > Here is my function:
> > >
> > > double
> > > distance(int x1, int y1, int x2, int y2)
> > > {
> > > double deltax = (double)x1 - x2;
> > > double deltay = (double)y1 - y2;
> > > return sqrt(deltax*del tax + deltay*deltay);
> > > }
> > >
> > >
> > > It works fine for most cases except that if deltax is really small
> > > compared to deltay (or vice versa), the addition inside the sqrt
> > > argument causes the smaller square value to be lost (because of adding
> > > a large float number with a small float number)).
> > >
> > > example: deltax = 1 deltay = 1999998
> > >
> > > What this distance function returns is exactly what would be returned
> > > if only one of the dimensions were 1999998.
> > >
> > > This is crucial because I need the distance formula to work within a
> > > certain tolerance so I can check to see if triangles are right or
> > > oblique.
> > >
> > > Does anyone know any tricks that I can use to preserve the "effect" of
> > > the small deltax in the distance?[/color]
> >
> > It ain't easy. What you need is a well written hypot function
> > from C99.[/color]
>
> I'm not sure if it's C99, but there is a hypot function in Sun Microsystem's
> Freely Distributable LIBM:
> http://www.netlib.org/fdlibm/
>[/color]


Thanks.

Now I still have a bug happening near the end of my program. When I do
the test for right triangle, I take the difference between (hyp^2) and
the sum of the other two sides squared. Sometimes this results in the
floating point value -0.00000 (i checked by printing it out). When I
compare the result with 0.0 it fails.

Any idea what this -0.0000 ? How can I compare that to 0?

Thanks

[serrand]

Re: floating point problems

G Patel wrote:[color=blue]
> Hello,
>
> I'm having trouble with floating point problems.
>
> I'm trying to write a function that calculates the distance between two
> cartesian points (integer coordinates). I have the function working,
> but it fails on one test case.
>
> Here is my function:
>
> double
> distance(int x1, int y1, int x2, int y2)
> {
> double deltax = (double)x1 - x2;
> double deltay = (double)y1 - y2;
> return sqrt(deltax*del tax + deltay*deltay);
> }
>
>
> It works fine for most cases except that if deltax is really small
> compared to deltay (or vice versa), the addition inside the sqrt
> argument causes the smaller square value to be lost (because of adding
> a large float number with a small float number)).
>
> example: deltax = 1 deltay = 1999998
>
> What this distance function returns is exactly what would be returned
> if only one of the dimensions were 19 999 98.
>
> This is crucial because I need the distance formula to work within a
> certain tolerance so I can check to see if triangles are right or
> oblique.
>
> Does anyone know any tricks that I can use to preserve the "effect" of
> the small deltax in the distance?
>
> Thanks
>[/color]

don't rely on the printing... with
printf ("dist ((%d, %d), (%d, %d)) == %14.8f\n", x1, y1, x2, y2, d);
i'va got the rigth result
dist ((0, 0), (1, 1999998)) == 1999998.0000002 5

Xavier

[Christian Bau]

Re: floating point problems

In article <1138919097.179 821.174410@g49g 2000cwa.googleg roups.com>,
"G Patel" <gaya.patel@gma il.com> wrote:
[color=blue]
> Now I still have a bug happening near the end of my program. When I do
> the test for right triangle, I take the difference between (hyp^2) and
> the sum of the other two sides squared. Sometimes this results in the
> floating point value -0.00000 (i checked by printing it out). When I
> compare the result with 0.0 it fails.[/color]

It looks like you are calculating a square root, and then you take the
result and square it immediately. That is a bit pointless, isn't it?
Just takes code and execution time and adds a bit of rounding error.
Better write a function that returns the square of the distance between
two points.

[Eric Sosman]

Re: floating point problems



G Patel wrote On 02/02/06 18:01,:[color=blue]
> [...]
>
> Thanks.
>
> Now I still have a bug happening near the end of my program. When I do
> the test for right triangle, I take the difference between (hyp^2) and
> the sum of the other two sides squared. Sometimes this results in the
> floating point value -0.00000 (i checked by printing it out). When I
> compare the result with 0.0 it fails.
>
> Any idea what this -0.0000 ? How can I compare that to 0?[/color]

It's a negative value whose magnitude is less
than half a ten-thousandth -- that is, it't a small
negative number that looks like zero when rounded
to four decimal places. Perhaps it's -0.00000312
or some such.

When comparing floating-point numbers, it is
usually a poor idea to check for exact equality.
F-P numbers should be regarded as approximations,
so rather than test for "equal" you should usually
test for "close enough." You should probably do
something like

if (fabs(result) < TOLERANCE) {
/* it's a right triangle, as near as
I can tell */
}
else {
/* it's definitely not a right triangle */
}

The test itself is simple, but choosing a good value
for TOLERANCE can be extraordinarily difficult.

See also Questions 14.1, 14.4, and 14.6 in the
comp.lang.c Frequently Asked Questions (FAQ) list at

http://www.c-faq.com/

for more information on floating-point comparisons.

--
Eric.Sosman@sun .com

[Keith Thompson]

Re: floating point problems

"G Patel" <gaya.patel@gma il.com> writes:
[...][color=blue]
> Now I still have a bug happening near the end of my program. When I do
> the test for right triangle, I take the difference between (hyp^2) and
> the sum of the other two sides squared. Sometimes this results in the
> floating point value -0.00000 (i checked by printing it out). When I
> compare the result with 0.0 it fails.
>
> Any idea what this -0.0000 ? How can I compare that to 0?[/color]

How did you print it? If you used printf's "%f" format, anything
sufficiently close to 0.0 will appear as either -0.000000 or 0.000000.
Use the "%e" or "%g" format to get a better idea of the actual value,
and specify extra precision to get an even better idea:

#include <stdio.h>
int main(void)
{
double x = -1.2345e-67;
double y = +1.2345e-67;
printf("x = %f, y = %f\n", x, y);
printf("x = %g, y = %g\n", x, y);
printf("x = %.20g, y = %.20g\n", x, y);
return 0;
}

Floating-point numbers are inexact. See section 14 of the
comp.lang.c FAQ, <http://www.c-faq.com/>.

--
Keith Thompson (The_Other_Keit h) kst-u@mib.org <http://www.ghoti.net/~kst>
San Diego Supercomputer Center <*> <http://users.sdsc.edu/~kst>
We must do something. This is something. Therefore, we must do this.

[Mark McIntyre]

Re: floating point problems

On 2 Feb 2006 15:01:21 -0800, in comp.lang.c , "G Patel"
<gaya.patel@gma il.com> wrote:

[color=blue]
>Any idea what this -0.0000 ? How can I compare that to 0?[/color]

Welcome to the Real World (tm), where maths is not exact.

Due to how computers store numbers, reals cannot generally be stored
precisely. Thus your maths has resulted in a number nearly (but not
quite) zero, lets say 1e-16. Google for "goldberg floating point" for
more detail.

As I'm sure you realise this makes it impossible in general terms to
compare floating point numbers. What you do instead is define some
EPSILON which is 'close enough' and compare the difference to that.

Mark McIntyre
--
"Debugging is twice as hard as writing the code in the first place.
Therefore, if you write the code as cleverly as possible, you are,
by definition, not smart enough to debug it."
--Brian Kernighan

----== Posted via Newsfeeds.Com - Unlimited-Unrestricted-Secure Usenet News==----
http://www.newsfeeds.com The #1 Newsgroup Service in the World! 120,000+ Newsgroups
----= East and West-Coast Server Farms - Total Privacy via Encryption =----

[Gregory Pietsch]

Re: floating point problems


G Patel wrote:[color=blue]
> Hello,
>
> I'm having trouble with floating point problems.
>
> I'm trying to write a function that calculates the distance between two
> cartesian points (integer coordinates). I have the function working,
> but it fails on one test case.
>
> Here is my function:
>
> double
> distance(int x1, int y1, int x2, int y2)
> {
> double deltax = (double)x1 - x2;
> double deltay = (double)y1 - y2;
> return sqrt(deltax*del tax + deltay*deltay);
> }
>
>
> It works fine for most cases except that if deltax is really small
> compared to deltay (or vice versa), the addition inside the sqrt
> argument causes the smaller square value to be lost (because of adding
> a large float number with a small float number)).
>
> example: deltax = 1 deltay = 1999998
>
> What this distance function returns is exactly what would be returned
> if only one of the dimensions were 1999998.
>
> This is crucial because I need the distance formula to work within a
> certain tolerance so I can check to see if triangles are right or
> oblique.
>
> Does anyone know any tricks that I can use to preserve the "effect" of
> the small deltax in the distance?[/color]

Let's see ... what you need looks something like this:

#include <errno.h>
#include <math.h>

double
hypot(double x, double y)
{
int errx = fpclassify (x);
int erry = fpclassify (y);
double fabsx, fabsy, z;

/* Inf takes precedence over NaN ... strange */
if (errx == FP_INFINITE || erry == FP_INFINITE)
{
errno = ERANGE;
return HUGE_VAL;
}
else if (errx == FP_NAN || erry == FP_NAN)
{
errno = EDOM;
return (errx == FP_NAN) ? x : y;
}
fabsx = x < 0.0 ? -x : x;
fabsy = y < 0.0 ? -y : y;
if (errx == FP_ZERO || erry == FP_ZERO)
{
if (errx == FP_ZERO && erry == FP_ZERO)
return +0.0;
else
return (errx == FP_ZERO) ? fabsy : fabsx;
}
else
{
if (fabsx > fabsy)
{
z = fabsy / fabsx;
return fabsx * sqrt (1.0 + z * z);
}
else
{
z = fabsx / fabsy;
return fabsy * sqrt (1.0 + z * z);
}
}
}
[color=blue]
>
> Thanks[/color]

HTH, HAND, Gregory Pietsch

[G Patel]

Re: floating point problems


Christian Bau wrote:[color=blue]
> In article <1138919097.179 821.174410@g49g 2000cwa.googleg roups.com>,
> "G Patel" <gaya.patel@gma il.com> wrote:
>[color=green]
> > Now I still have a bug happening near the end of my program. When I do
> > the test for right triangle, I take the difference between (hyp^2) and
> > the sum of the other two sides squared. Sometimes this results in the
> > floating point value -0.00000 (i checked by printing it out). When I
> > compare the result with 0.0 it fails.[/color]
>
> It looks like you are calculating a square root, and then you take the
> result and square it immediately. That is a bit pointless, isn't it?
> Just takes code and execution time and adds a bit of rounding error.
> Better write a function that returns the square of the distance between
> two points.[/color]

My hands are tied because I have to use the distance function's
distances inside my right triangle function.

[Christian Bau]

Re: floating point problems

In article <1138924312.688 876.80610@o13g2 000cwo.googlegr oups.com>,
"Gregory Pietsch" <GKP1@flash.net > wrote:
[color=blue]
> double
> hypot(double x, double y)
> {
> int errx = fpclassify (x);
> int erry = fpclassify (y);
> double fabsx, fabsy, z;
>
> /* Inf takes precedence over NaN ... strange */
> if (errx == FP_INFINITE || erry == FP_INFINITE)
> {
> errno = ERANGE;
> return HUGE_VAL;
> }
> else if (errx == FP_NAN || erry == FP_NAN)
> {
> errno = EDOM;
> return (errx == FP_NAN) ? x : y;
> }
> fabsx = x < 0.0 ? -x : x;
> fabsy = y < 0.0 ? -y : y;
> if (errx == FP_ZERO || erry == FP_ZERO)
> {
> if (errx == FP_ZERO && erry == FP_ZERO)
> return +0.0;
> else
> return (errx == FP_ZERO) ? fabsy : fabsx;
> }
> else
> {
> if (fabsx > fabsy)
> {
> z = fabsy / fabsx;
> return fabsx * sqrt (1.0 + z * z);
> }
> else
> {
> z = fabsx / fabsy;
> return fabsy * sqrt (1.0 + z * z);
> }
> }
> }[/color]

This implementation has the rather unfortunate property that there are
values x1, x2, y1, y2 such that

0 <= x1 <= x2
0 <= y1 <= y2

but

hypot (x1, y1) > hypot (x2, y2)

To see how this happens: Choose y1 = y2 = 1. Start with x = 1.99. The
function will calculate z = 1.0 / x, then tmp = sqrt (1.0 + z*z), then
finally return x * tmp. z will be around 0.5, tmp will be about 1.12.

Let u be the difference between x and the next larger floating point
number. If you do the same calculation with x+u, x+2u, x+3u etc. instead
of x, the value of z will eventually become smaller by u, so we can find
values x1, x2 such that x2 = x1 + u and tmp (x2) = tmp (x1) - u. In that
case the product fabsx * sqrt (1.0 + z * z), where sqrt (1.0 + z * z)
was calculated with double precision, will actually become smaller. Due
to the rounding of the result, the result will not always become
smaller, but sometimes it will become smaller when x is made larger.

I also suspect that you will not get the correctly rounded result in
some trivial cases like hypot (5.0, 12.0) which should be exactly equal
to 13.0.

You get a better and faster result as follows:

1. Calculate z = x*x + y*y. Make sure that the implementation doesn't
use fused multiply-add instructions for this; otherwise it can happen
that hypot (x, y) != hypot (y, x) which would be unfortunate.

2. Check that the result of z did not overflow, but also check that it
is large enough to guarantee that either none of the products x*x and
y*y did underflow, or that the other product is so large that it doesn't
make a difference. If both conditions are true, return sqrt (z).

3. If the result of z did overflow, you can find a constant c which is a
power of two so that c * sqrt ((x / c) * (x / c) + (y / c) * (y / c))
will give the correctly rounded result. The same can be done with a
different constant c' if the result of z was too small.

[Christian Bau]

Re: floating point problems

In article <1138925921.972 033.243230@g44g 2000cwa.googleg roups.com>,
"G Patel" <gaya.patel@gma il.com> wrote:
[color=blue]
> Christian Bau wrote:[color=green]
> > In article <1138919097.179 821.174410@g49g 2000cwa.googleg roups.com>,
> > "G Patel" <gaya.patel@gma il.com> wrote:
> >[color=darkred]
> > > Now I still have a bug happening near the end of my program. When I do
> > > the test for right triangle, I take the difference between (hyp^2) and
> > > the sum of the other two sides squared. Sometimes this results in the
> > > floating point value -0.00000 (i checked by printing it out). When I
> > > compare the result with 0.0 it fails.[/color]
> >
> > It looks like you are calculating a square root, and then you take the
> > result and square it immediately. That is a bit pointless, isn't it?
> > Just takes code and execution time and adds a bit of rounding error.
> > Better write a function that returns the square of the distance between
> > two points.[/color]
>
> My hands are tied because I have to use the distance function's
> distances inside my right triangle function.[/color]

If it is a homework problem, tell your instructor that his method is
stupid. If it is not a homework problem, your hands are not tied.

If you are given three points (x0, y0), (x1, y1), and (x2, y2), the
easiest way to determine whether there is a right angle at (x0, y0) is
the following:

Subtract (x0, y0) from (x1, y1) and (x2, y2), so the whole triangle is
moved to bring the first point to (0, 0)

x1 -= x0; y1 -= y0;
x2 -= x0; y2 -= y1;

Then all you do is to check whether

x1*y1 == - x2*y2

If they are the same, you have a right angle. If they are not the same,
you don't have a right angle. That's all.

[G Patel]

Re: floating point problems


Christian Bau wrote:[color=blue]
> In article <1138925921.972 033.243230@g44g 2000cwa.googleg roups.com>,
> "G Patel" <gaya.patel@gma il.com> wrote:
>[color=green]
> > Christian Bau wrote:[color=darkred]
> > > In article <1138919097.179 821.174410@g49g 2000cwa.googleg roups.com>,
> > > "G Patel" <gaya.patel@gma il.com> wrote:
> > >
> > > > Now I still have a bug happening near the end of my program. When I do
> > > > the test for right triangle, I take the difference between (hyp^2) and
> > > > the sum of the other two sides squared. Sometimes this results in the
> > > > floating point value -0.00000 (i checked by printing it out). When I
> > > > compare the result with 0.0 it fails.
> > >
> > > It looks like you are calculating a square root, and then you take the
> > > result and square it immediately. That is a bit pointless, isn't it?
> > > Just takes code and execution time and adds a bit of rounding error.
> > > Better write a function that returns the square of the distance between
> > > two points.[/color]
> >
> > My hands are tied because I have to use the distance function's
> > distances inside my right triangle function.[/color]
>
> If it is a homework problem, tell your instructor that his method is
> stupid. If it is not a homework problem, your hands are not tied.
>[/color]

It is a school assignment and I have the program working except that
our programs are testing automatically by some program.

He lets us know 8 of the cases, the rest are hidden. 2 of the cases
seem impossible to implement in computers (and yes, we have to sqrt in
a distance function, and then resquare it in the actualy "figuring out"
function).

One of the cases has a triangle with a side that is of length sqrt(41),
which is irrational. So when the other function gets it and resquares
it, accuracy is lost. Yet his testing program wants an exactly
solution (no approximate one allowed for this case).

The other case has a triangle with 2 really large sides and one very
small sides. He wants us to get it within some tolerance, but i can't
seem to "keep" the information from the small side to "stick" within
the sqrt function.

I am going to email him and tell him that his cases are ridiculous (in
a nice way).

Thanks.

[Coos Haak]

Re: floating point problems

Op 2 Feb 2006 16:58:41 -0800 schreef G Patel:
[color=blue]
> One of the cases has a triangle with a side that is of length sqrt(41),
> which is irrational. So when the other function gets it and resquares
> it, accuracy is lost. Yet his testing program wants an exactly
> solution (no approximate one allowed for this case).[/color]

<OT>
This is a rectacngular triangle with sides 4, 5 and sqrt(41).
See your Math book!
</OT>
--
Coos