Weird numerical behavior regarding zero

[matthias_k]

Hi,

I am currently implementing a solver for linear optimization problems
using the revised simplex method and I stumbled accross some strange
behavior regarding the treatment of the number 0.
I am not sure whether this is compiler or language related though.

The first problem is, that there are two representations of zero, namely
+0 and -0. Does IEEE for floating points allow this? Maybe I didn't
activate IEEE floating point numbers?
Anyway, this does have some very annoying consequences, since this
expression:
x / -0 yields -inf (x is positive).

Which brings me to my second question. Shouldn't division by zero make
my computer go up in flames or so? Because, on my system (Arch Linux 0.7
i686) with my compiler (g++ 3.4.3), dividing by zero is allowed.
In fact, it yields +inf when dividing by +0 and -inf when dividing by
-0, given that the nominator (is that the english term for "Zähler"?) is
positive of course.

Any comments/ideas?

Regards,
Matthias

 

[Andrew Koenig]

The first problem is, that there are two representations of zero, namely
+0 and -0. Does IEEE for floating points allow this? Maybe I didn't
activate IEEE floating point numbers?

IEEE floating point requires +0 and -0 as distinct representations.

Anyway, this does have some very annoying consequences, since this
expression:
x / -0 yields -inf (x is positive).

Yes, that is what IEEE floating point requires.

Which brings me to my second question. Shouldn't division by zero make my
computer go up in flames or so? Because, on my system (Arch Linux 0.7
i686) with my compiler (g++ 3.4.3), dividing by zero is allowed.

C++ says that what happens after division by zero is undefined.

IEEE says that the result is either a floating-point trap or an
appropriately signed infinity, depending on how you have set things up.

The behavior you have described so far is completely normal according to
IEEE.

 

[matthias_k]

IEEE floating point requires +0 and -0 as distinct representations.




Yes, that is what IEEE floating point requires.




C++ says that what happens after division by zero is undefined.

IEEE says that the result is either a floating-point trap or an
appropriately signed infinity, depending on how you have set things up.

The behavior you have described so far is completely normal according to
IEEE.

Thanks Andrew. That is terrible news.
Does that mean I have to multiply each -0 with -1 to get a sane
representation of zero?

 

[Dave Moore]

Hi,

I am currently implementing a solver for linear optimization problems
using the revised simplex method and I stumbled accross some strange
behavior regarding the treatment of the number 0.
I am not sure whether this is compiler or language related though.

The first problem is, that there are two representations of zero, namely
+0 and -0.

Depending on how you are obtaining your "signed zero" values, many different
things could be happening. You should use the constants defined in
std::numeric_limits against your values to figure out what is happening.
For example, if the problem values results from subtraction operations, you
have to check them against the machine accuracy (e.g.
std::numeric_limits<float>::epsilon()) to make sure they haven't lost all
significance. I guess this is what happened in your case, and the sign bit
is just staying set.


HTH,

Dave Moore

 

[matthias_k]

Depending on how you are obtaining your "signed zero" values, many different
things could be happening. You should use the constants defined in
std::numeric_limits against your values to figure out what is happening.
For example, if the problem values results from subtraction operations, you
have to check them against the machine accuracy (e.g.
std::numeric_limits<float>::epsilon()) to make sure they haven't lost all
significance. I guess this is what happened in your case, and the sign bit
is just staying set.


HTH,

Dave Moore

I am doing lots of subtraction, multiplication and division.
What do you mean by check against the machine accuracy? Does -0 imply
that it's not 0 but something really close to 0? I can't check that for
each arithmetical operation, that's too expensive.

 

[Pete Becker]

Does that mean I have to multiply each -0 with -1 to get a sane
representation of zero?

No, it means that you have to rethink what you mean by zero. <g>

The problem is that for negative values of x, 1/x is also negative. As x
approaches zero from below, 1/x becomes a larger and larger negative
value. The limit of 1/x as x approaches zero from below is negative
infinity. Similarly, for positive valus of y, 1/y is also positive. As y
approaches zero from above, 1/y becomes a larger and larger positive
value. The limit of 1/y as y approaches zero from above is positive
infinity. If you don't distinguish between negative and positive zeroes
you end up with the limits of both of those expressions being the same.
While it's possible to do math that way, it's far more intuitive to have
negative and positive infinity as two distinct values, and that in turn
means that you need two distinct zeroes, one positive and one negative.

 

[Jerry Coffin]

matthias_k wrote:

[ ... ]

The first problem is, that there are two representations of zero, namely
+0 and -0. Does IEEE for floating points allow this? Maybe I didn't
activate IEEE floating point numbers?

The C and C++ standards allow this. The IEEE floating point standards
_require_ it, if memory serves at all.

Anyway, this does have some very annoying consequences, since this
expression:
x / -0 yields -inf (x is positive).

Which brings me to my second question. Shouldn't division by zero make
my computer go up in flames or so?

According to C it gives undefined behavior -- which can include
returning a value it considers reasonable. C99 added a number of
functions for controlling how floating point is done. AFAIK, they
haven't officially been adopted into C++ yet, but your compiler may
include them anyway -- if you have an fenv.h, chances are you can use
its functions from C++ as well as from C. I could be wrong, but at
first glance I'd expect FE_OVERFLOW to apply to division by 0.

Because, on my system (Arch Linux 0.7
i686) with my compiler (g++ 3.4.3), dividing by zero is allowed.
In fact, it yields +inf when dividing by +0 and -inf when dividing by

-0, given that the nominator (is that the english term for "Zähler"?) is
positive of course.

"Numerator", is the word you seem to be looking for. In any case, as I
said above, the C++ standard simply makes division by 0 undefined
behavior, so anthing is possible.

The usual way to deal with negative zeros is not by multiplication, but
by addition -- a negative zero plus a positive zero gives a positive
zero, so to remove negative zeros, you add 0.0 to anything that might
have resulted in a negative zero.

 

[Karl Heinz Buchegger]

I am doing lots of subtraction, multiplication and division.
What do you mean by check against the machine accuracy? Does -0 imply
that it's not 0 but something really close to 0? I can't check that for
each arithmetical operation, that's too expensive.

Taking your comments into account I would say that lots
of surprises are waiting for you.

First of all, you really shouldn't divide by 0. No matter if
IEEE allows it or not. Simply don't do it. If you are bound
to divide by 0, then something strange has happend and your
program has to react on that.

What Dave means. How do you know that something is exactly
+0 or -0. In programmers real life you can't distinguish between
them, because you don't know if the sign is correct or not.
If eg. the result of a subtraction or addition happens to be
-0 then you never know if the sign is indeed '-' or if it
(mathematically) should be '+'. This is because of small round
of errors which happen all the time in floating point arithmetic.
So the best you can do is:
do the arithmetic and before you you perform the division,
check if you don't divide by a number that is small enough
to be assumed to be 0. If you figure out that you are going
to do that, consider it to be an error and react to it.
Do *not* depend on the result to be +inf or -inf, those
results will be wrong on most of the cases anyway.

In a computer 2 * 0.1 does *not* equal 0.2, but something that
is very close to 0.2 but not exactly 0.2
You need to take that into account with every floating point
operation you make.

And no: checking for division by 0.0 is *never* too expensive
no matter how many operations you have to perform. There is only
one exception: If you can proove that your division will never
be executed with something that is even close to 0.0. In all
other cases: Do the check, but do it correct:

if( fabs( some_number ) < some_epsilon ) {
// some_number is close to 0.0
// treat it as 0.0 and react to it
}
else
op = 5.0 / some_number;

As for the value of some_epsilon. It mostly depends on the
history of some_number, how many 'floating point debris'
has accumulated in it. This is dependent on
a) the number of operations
b) the range of the numbers participating in those operations.

Using a value like std::numeric_limits<float>::epsilon() isn't going
to help.

And by the way: Never use float, until you know what you
do, have the knowledge to fight that beast, are willing to fight
that beast and have a very, very, very, very good reason to use
it. In all other cases (that is 99.99% of all cases) use double.

 

[Pete Becker]

And no: checking for division by 0.0 is *never* too expensive
no matter how many operations you have to perform.

The reason for having both forms of zero, infinities, and NaNs is so
that you don't have to check before every operation. Once an error
occurs it propogates through the subsequent calculations, so you only
have to check at the end. The underlying assumption is that errors are
rare, so not checking for them makes complicated computations much faster.

 

[White Wolf]

No, it means that you have to rethink what you mean by zero. <g>

[SNIP]

Just as an (off-topic) sidenote to that: (probably) that definition of
floating point zero has kept a lot of people wondering for quite some time
in Hungary. A public digital thermometer was installed at the capital, and
it was sometimes showing -0 Celsius (Centigrade) and sometimes +0. People
came up with loads of explanations as to why that thermometer made a
difference between -0 and +0 degrees, why was it important. I guess that
Forest Gump was closest to the solution when he said: for no particular
reason.

 

[Karl Heinz Buchegger]

The reason for having both forms of zero, infinities, and NaNs is so
that you don't have to check before every operation. Once an error
occurs it propogates through the subsequent calculations, so you only
have to check at the end. The underlying assumption is that errors are
rare, so not checking for them makes complicated computations much faster.

Understood.
AS for me (working in computer graphics): Attempting to divide
by 0.0 almost always means that some special case occours that
needs to be handled seperately. Eg. The formula for computing
the intersection point of 2 line segments contains a division.
Attempting to divide by 0.0 means that those lines are parallel
and no intersection is possible. So for me, premature testing is
way simpler, because it means that I don't have to figure out
afterwards that something special has occoured and the calculated
point is meaningless.
Usually in my area those 0.0-tests are indeed rare, once I know that
the parameters are in allowed limits I know for sure that the division
will work and omit further tests.

But others milage may (and will) vary, of course.

The worst case however is: The program crashes and the programmer,
when asked why he did not check for this case, answers: "Because I
hunted for speed." That's why I said: "*never* too expensive"

 

[Pete Becker]

Just as an (off-topic) sidenote to that: (probably) that definition of
floating point zero has kept a lot of people wondering for quite some time
in Hungary. A public digital thermometer was installed at the capital, and
it was sometimes showing -0 Celsius (Centigrade) and sometimes +0. People
came up with loads of explanations as to why that thermometer made a
difference between -0 and +0 degrees, why was it important. I guess that
Forest Gump was closest to the solution when he said: for no particular
reason.

Actually, in that case I'll bet it's rounding a (small) non-zero value
to zero and keeping the sign. That's a coding error.

 

[Matthias]

Actually, in that case I'll bet it's rounding a (small) non-zero value
to zero and keeping the sign. That's a coding error.

I have now introduced a function called check_zero:

template< typename FPT >
bool check_zero( FPT& value )
{
if( fabs(value) < EPSILON )
{
value = 0;
return true;
}
return false;
}


I hope this does the job.

 

[Lionel B]

Hi,

I am currently implementing a solver for linear optimization problems
using the revised simplex method and I stumbled accross some strange
behavior regarding the treatment of the number 0.
I am not sure whether this is compiler or language related though.

The first problem is, that there are two representations of zero, namely
+0 and -0. Does IEEE for floating points allow this? Maybe I didn't
activate IEEE floating point numbers?
Anyway, this does have some very annoying consequences, since this
expression:
x / -0 yields -inf (x is positive).

Which brings me to my second question. Shouldn't division by zero make
my computer go up in flames or so? Because, on my system (Arch Linux 0.7
i686) with my compiler (g++ 3.4.3), dividing by zero is allowed.
In fact, it yields +inf when dividing by +0 and -inf when dividing by
-0, given that the nominator (is that the english term for "Zähler"?) is
positive of course.

Any comments/ideas?

The following (admittedly rather scary) document should tell you more than
you wanted to know about the above (and many other) floating point issues:
"What Every Computer Scientist Should Know About Floating-Point Arithmetic".
Ubiquitous on the net, here's one link to a pdf:
http://docs-pdf.sun.com/800-7895/800-7895.pdf

Regards,

 

[Pete Becker]

I have now introduced a function called check_zero:

template< typename FPT >
bool check_zero( FPT& value )
{
if( fabs(value) < EPSILON )
{
value = 0;
return true;
}
return false;
}


I hope this does the job.

It will do whatever it does. If that's what you want, go for it. But a
better name might be "force_small_values_to_zero". Close to zero is not
the same as zero. You haven't shown the definition of EPSILON, so it's
hard to see what this really does, but in general the potential error in
a value depends on how that value was calculated. Slamming all values
within a fixed range to zero is rarely appropriate.

 

[Matthias]

It will do whatever it does. If that's what you want, go for it. But a
better name might be "force_small_values_to_zero". Close to zero is not
the same as zero. You haven't shown the definition of EPSILON, so it's
hard to see what this really does, but in general the potential error in
a value depends on how that value was calculated. Slamming all values
within a fixed range to zero is rarely appropriate.

I have set EPSILON to 1.0e-14 for now. This value is larger than my
system's epsilon, but I had an example input where the objective
function value of the linear problem was something like 1.68e-15 (so it
was zero...) and because it was not EXACTLY zero my program went into
another iteration just to find out that it can't continue.

As to the rounding erros... What else can I check besides tiny values
around zero? If I would know what the correct result of a floating point
operation would be, I wouldn't need to calculate it...
So how could I check if I lost any precision after a calculation?

 

[Pete Becker]

I have set EPSILON to 1.0e-14 for now. This value is larger than my
system's epsilon, but I had an example input where the objective
function value of the linear problem was something like 1.68e-15 (so it
was zero...) and because it was not EXACTLY zero my program went into
another iteration just to find out that it can't continue.

Just a note: the value of epsilon depends on the type of the floating
point value that you're using.

As to the rounding erros... What else can I check besides tiny values
around zero? If I would know what the correct result of a floating point
operation would be, I wouldn't need to calculate it...
So how could I check if I lost any precision after a calculation?

<g> Yes, that's the problem. People get PhD's in this sort of thing. I
don't know much about it, but in general it's a matter of careful
algorithm design and tweaking. For example, if you're solving for the
roots of quadratic equations you have the usual quadratic formula:

x1 = (-b + sqrt(b^2 - 4ac)) / 2a

and

x2 = (-b - sqrt(b^2 - 4ac)) / 2a

But you never compute the two roots this way. If b is negative you
calculate x1 first. You don't use the formula for x2, because it ends up
subtracting two values (-b is positive, and you'd subtract the result of
sqrt). That avoids magnifying roundoff errors. From x1 you calculate x2
with a different formula (which escapes me at the moment). Similarly, if
b is positive you calculate x2 first, again avoiding subtraction (factor
out the -1), and calculate x1 from that result.

You might want to look at www.petebecker.com/journeymansshop.html for a
more detailed discussion of some of this stuff.

 

[Karl Heinz Buchegger]

As to the rounding erros... What else can I check besides tiny values
around zero? If I would know what the correct result of a floating point
operation would be, I wouldn't need to calculate it...
So how could I check if I lost any precision after a calculation?

The problem with this type of things is that there is almost
no 'one size fits all' rule. Each problem has to be addressed
indivdually.

For another example (taken from "Geometric and Solid Modeling, Christoph Hoffmann" )

Consider implementing a test of whether two points in the plane
are equal. Specifically assume that the point u is the intersection
of the pair of lines (L1,L2), and that the point v is the intersection
of the lines (L3,L4). The line equations are the input to the following
algorithm:

1. Compute the coordinates of u.
2. By substituting into the line equations L3 and L4, conclude that
u == v if both L3(u) and L4(u) are smaller then some tolerance.

Intuitively, this algorithm ought to be equivalent to a second version in
which the roles of u and v are reversed. Lets see if this is true.

(Hint: a1 denotes 'a' with a subcscript of 1)

1. The intersectoin (ux, uy) of the lines

a1*x + b1*y + c1 = 0
a2*x + b2*y + c2 = 0

is computed as:

D = a1*b2 - a2*b1
ux = (b1*c2 - b2*c1) / D
uy = (a2*c1 - a1*c2) / D

2. The point (ux,uy) is assumed to lie on the line a*x + b*y + c = 0 if the
distance is small; that is if |a * ux + b * uy + c| < eps * sqrt(a*a + b*b)

We assume eps to be 1E-10, a reasonable bound for double precision. We ask whether u
and v are incident using 2 different methods:

a. Compute the coordinates of u; conclude that u == v iff u is on both L3 and L4.
b. Compute the coordinates of v; comclude that u == v iff v is on both L1 and L2.

The line coefficients follow. Since pow(2,-23) ~ 1E-7, they differ from 1 and 0
by amounts that are several orders of magnitude larger then eps.

a1 = -1 b1 = 1 c1 = 0
a2 = -(1+pow(2,-23)) b2 = 1-pow(2,-23) c2 = pow(2,-22)
a3 = 1 b3 = 0 c3 = -(1+pow(2,-15))
a4 = 0 b4 = 1 c4 = -(1+pow(2,-15))

These coefficients can be represented *exactly* in double precision. The coordinates
of the points are now computed to be

u = ( 1.0, 1.0 )
v = ( 1.000030517578125, 1.000030517578125)

They are both exact. Moreover since ai*ai + bi*bi is aproximately between
1 and 2, the evaluation of the line equations after substituting the point coordinates
yields an error that can be compared directly with eps. We obtain the values

L3(u) ~ -3E-5 > eps
L4(u) ~ -3E-5 > eps

from which we conclude that u connot be incident to v, since it is too far from
the lines L3, L4 whose intersection is v.

But we also obtain

L1(v) = 0 < eps
L2(v) ~ -7E-12 < eps

from which we must conclude that v is incident to u, since it lies extremely close
to both lines.
Therefore, although they ask the same geometric question, the two computations
yield contradictory results.

So what can one do about this.
Clearly: Don't use that algorithm. One could choose to do:

1. Compute the coordinates of u and v, by intersecting the respective lines.
2. If the euclidian distance between u and v is smaller then eps, decide u == v;
otherwise decide u != v

(Note that this computation is more 'expensive' then the previous one. This time
both points and an euclidian distance needs to be computed. In the previous
algorithm only one point is computed and the testing does not involve a square
root.)

This method is symetric, but it does not exhibit transitivity. Specifically,
we choose 3 points u, v and w such that u is incident to v, v is incident to w,
*but* u is not incident to w. We assume an eps of 1E-10

u = ( 0, 0 ) v = ( 0, 0.8E-10) w = ( 0, 1.6E-10 )

Clearly the distance between the adjacent pairs is less then eps, but the distance
between u and w is greater then eps. So what is the correct answer? Are all those
points equal or are they not?