Float comparison

[Keith Thompson]

Keith Thompson wrote:

... snip ...

And if I write this:

double x = 1.0;
if (x < 1.0) puts("Oops!");
if (x > 1.0) puts("Oops!");

do I have any grounds for complaint if the program prints
"Oops!"? If the stored value represents a range, I would
think the answer would be no.

x is a double. So is the expression (1.0). So they represent
the same range.

[...]

I think you're saying that the above code may not print "Oops!",
but you didn't actually say so. Can you confirm that?

As written, yes. With the modification I wrote, no.

In the following, I will ignore NaNs and infinities.

So in your model, (x < y) must yield 0 if both x and y represent
the same range, even though some values in the range of x are less
than some values in the range of y. Presumably (x < y) yields 1 if
and only if each value in the range of x is less than each value
in the range of y. And (x == y) if and only if both ranges are
exactly the same. Is that correct?

No. If x and y contain the identical ranges, there are no values
in the x range outside the y range, and vice-versa. So your
statement is impossible. The only thing that controls the ranges
is the value in x (or y). The actual size of the ranges is FP
system dependant.

I believe you have misunderstood what I wrote.

Let's suppose that, in accordance with your model, the stored double
value 1.0 represents the range of real numbers 1.0-EPS .. 1.0+EPS
assuming for simplicity that the range is symmetric). Then given:
double x = 1.0;
double y = 1.0;
both x and y represent identical ranges. Each of these identical
ranges contains infinitely many real numbers. Some of the numbers
within x's range are less that some of the numbers within y's range;
for example, 1.0-EPS (which is in x's range) is less than 1.0+EPS
(which is in y's range). Nevertheless, in your model, (x < y) must
yield 0.

And just what do you mean by "the value in x (or y)"? Are you
implying that x has just a single value?

Incidentally, in my article there were 5 paragraphs following what you
quoted, full of questions directed to you in an attempt to understand
your conception of C's floating-point model. You have snipped those
paragraphs without so much as a "[...]". My assumption is that you
were unable to answer my questions because your floating-point model
is actually inconsistent. Either that or you didn't bother to read
them.

 

[Keith Thompson]

True, but we are talking about FP representation (and int
representation), which have other limitations, including maximum
and minimum values. ints don't have a range associated with each
value, the adjacent integer value is always available. Reals don't
have an adjacent value, but rationals do (they are countable). The
FP system is built on rationals, but cannot express the adjacent
rational. Neither the integer nor the FP system can handle the
infinities that exist in arithmetic.

Yes, rationals are countable, but no, they do not have adjacent
values. Given any two rational numbers x and y, the rational number
(x+y)/2 is between them and distinct from both of them. (I'm
referring to mathematical rational numbers here, not to any C
representation of them.)

However floating-point types do have adjacent values. For example,
1.0 and 1.0+DBL_EPSILON are adjacent values of type double; there are
no double values between them.

 

[Guest]

(e-mail address removed) said:


It's a poor analogy.

oddly, I disagree

If I were to argue that the Standard requires
int to be exactly 16 bits (which, incidentally, is obviously not
the case), and you were to argue otherwise, would you call it a
draw if each of us stated our case "a few times"?

no. but if you kept on saying it and apparently were not
reading my rebuttals I'd give up.

You come from the view-point that people are fundametally rational
and if you keep working at it they will evantually grok you argument.
I consider people to be fundamentally irrational and if you don't,
quickly, get that brief flash of rationality you've had it.

Man is a rationlising animal. We form opinions and then make up
reasons for the opinions. You get less of this in technical news
groups but scratch a C programmer and you'll find an underlying
human being.

In a chess
context, the argument in this thread bears a closer parallel to "I
can use en passant even if other moves have happened since the
opponent advanced his pawn two places", "no you can't and here's
why", "yes I can", "no you can't", ad nauseam.


Sure. So does a killfile. So does skiing in Switzerland. So does
spacewalking. And, like those, it doesn't actually move the
argument any further forward.

sometimes it *can't* be moved forwards

Unfortunately, there doesn't seem to
be a way to move the argument forward. Chuck is clearly wrong, and
his opponents are clearly right for the reasons they have stated,
and his inability to accept this is his problem, not theirs.

Oddly enough, however, I think that for once Chuck is at least
/almost/ right. He's wrong about *_EPSILON, but I can understand
his point about the limitations of floating-point numbers. (I think
we all understand this, actually.) His problem is not so much in
what he's thinking as in the words he's using to think it.

so much for Worf's Hypothesis...

 

[Keith Thompson]

I still can't translate that into being a binary predicate, alas.

I think the claim is something like this:

Given a floating-point representation with a 40-bit mantissa, an
attempt to find a solution to a system of linear equations involving
more than about 40 variables will involve so many rounding errors in
intermediate results that the cumulative error in the final result
will be so large as to make that result meaningless.

I don't know whether there's any direct relationship between the two
occurrences of 40 (bits vs. variables). I have no informed opinion on
whether the claim is correct, other than an assumption that von
Neumann wouldn't have gotten something like this wrong.

 

[Keith Thompson]

You snipped the whole explanation. I recommend you don't snip in
the middle of paragraphs.

What explanation?

Here's the full quote:
| Try thinking about it. For example, if I write:
|
| float x = 1.0/3.0;
|
| you are claiming the value stored is 1/3. But if you examine the
| value stored you will not find that. You will find something in
| the range x*(1+FLT_EPSILON) and x*(1-FLT_EPSILON). That range is
| system dependant. We won't move that range very much if we
| substitute the desired value, 1/3, for the value found in x.
| Things may be worse, but they won't be better.

Assuming that the "." in "1/3." was intended to be a period, not a
decimal point, your statement "you are claiming the value stored is
1/3" is simply false. Nothing that follows justifies or explains it.

Assuming binary floating-point, the mathematical value 1/3 is not
representable in type float. The value stored is not 1/3; it's a
representable value close to the real number 1/3.

 

[Guest]

(e-mail address removed) writes:

I suspect the spelling "Heisenburg" has arisen as a superposition of
states in which the eigenstates are the strings "Heisenberg" and
"Heisenbug".  Measuring this spelling in comp.lang.c will probably cause
a collapse to the state |"Heisenbug">.  Measuring it in a group devoted
to discussion of the German physicist Werner Heisenberg would probably
have a different outcome.  Of course, if you are a Many Worlds devotee
(or subscribe to both newsgroups), both could happen.  

Quantum Boggum Sort:
Q1. use a source of quantum noise (eg. radioactive decay) to
randomly permutate an array.
Q2. if the array is not ordered, destroy the universe (*)
Q3. if you reached this step your universe has sorted the array
in O(n) time.
(*) [100] this is left as an exercise

 

[James Kuyper]

Yes, rationals are countable, but no, they do not have adjacent
values.

The fact that they are countable means that they can be arranged in a
sequence that allows them to be mapped one-to-one with the integers.
It's therefore meaningful, within the context of such a mapping, to talk
about "adjacent" rationals. However, there are infinitely many different
such sequences, and the meaning of "adjacent" would be different in each
one. One of the few things all such sequence have in common is that they
cannot be in order by the values of the rational numbers. Those two
facts makes utter nonsense of the conclusions CBFalconer reaches from
the countability of rationals.

 

[Chris Dollin]

Keith Thompson said:



(3/2 + 3/2) / 2 = 3/2

ITYM any two /different/ rational numbers.

If x and y are equal, they're not "two rational numbers x and y";
the /different/ is encoded in the /two/.

 

[James Kuyper]

If x and y are equal, they're not "two rational numbers x and y";
the /different/ is encoded in the /two/.

If that were the case, then the phrase "two identical numbers" would
count as a contradiction in terms.

 

[Dik T. Winter]

> Keith Thompson wrote: ....
>
> Oh? I have a nuclear counter and a sensor. Over a given period T
> I detect 10,000 counts of events of more than 1 keV energy. I
> compute the expected error in the counts as sqrt(count), or 100 and
> store that in cterror. Now you claim that 'counts' is exact, and
> represents the radiation strength at the time of testing? You are
> telling me to ignore cterror? I haven't even considered the time
> resolution of the detection system.

So you think the calculation of sqrt(count) is exact? I would say the
error is in the calculation and *not* anywhere else.

 

[Dik T. Winter]

>
> I think the claim is something like this:
>
> Given a floating-point representation with a 40-bit mantissa, an
> attempt to find a solution to a system of linear equations involving
> more than about 40 variables will involve so many rounding errors in
> intermediate results that the cumulative error in the final result
> will be so large as to make that result meaningless.
>
> I don't know whether there's any direct relationship between the two
> occurrences of 40 (bits vs. variables). I have no informed opinion on
> whether the claim is correct, other than an assumption that von
> Neumann wouldn't have gotten something like this wrong.

The two occurrences of 40 are just an accident, and as I wrote, von Neumann
was right, as long as you do forward error analysis. In numerical algebra
you need backward error analysis to know what the result actually means.

 

[Dik T. Winter]

> Reals don't
> have an adjacent value, but rationals do (they are countable).

What nonsense is this last? Do you have any knowledge of mathematics?
What is (in the rationals) the adjacent rational to 1/2?

Yes, they are countable, but that does *not* mean that there is an
adjacent one in the standard order.

 

[Dik T. Winter]

> "Dik T. Winter" wrote: ....
>
> I disagree. What is useful depends on the use to which you put the
> results. However the innate construction of FP values is such that
> the 'range' always exists. That range may be so small that you can
> afford to ignore it, but you can't just arbitrarily ignore it.

I have written enougb software in numerical mathematics that just did
ignore it. If I had assumed the values represented ranges that software
would have become useless.

> For example the number of different values is, at most, the binary
> number expressed by the size of a float. For each such value you
> can express the central value, in some form, and the next and
> previous values. You can order them all (again, ignoring NANs and
> INFs etc.).

But that is error analysis using interval arithmetic.

> YOU CAN'T express any value outside that list, but the
> available 'range' around the values allows you to REPRESENT any
> value sufficiently well to perform calculations.

And those ranges are useless in error analysis because after only a single
calculation the range will already be too small, as in z = x + y.

 

[Chris Dollin]

If that were the case, then the phrase "two identical numbers" would
count as a contradiction in terms.

Or as a way of saying "two different representations of the same number".
But yes: "two identical numbers" definitely has some kind of contradiction
seething inside.

 

[Keith Thompson]

You come from the view-point that people are fundametally rational
and if you keep working at it they will evantually grok you argument.
I consider people to be fundamentally irrational and if you don't,
quickly, get that brief flash of rationality you've had it.

[...]

I refuse to believe that, no matter how many times you explain it.

}

 

[Beej Jorgensen]

Have you missed the posts to this thread by CBF? He has disputed that!

But the post I replied to implied the opposite, and I wanted
verification.

-Beej

 

[CBFalconer]

.... snip ...

I believe you have misunderstood what I wrote.

Let's suppose that, in accordance with your model, the stored double
value 1.0 represents the range of real numbers 1.0-EPS .. 1.0+EPS
assuming for simplicity that the range is symmetric). Then given:
double x = 1.0;
double y = 1.0;
both x and y represent identical ranges. Each of these identical
ranges contains infinitely many real numbers. Some of the numbers
within x's range are less that some of the numbers within y's range;
for example, 1.0-EPS (which is in x's range) is less than 1.0+EPS
(which is in y's range). Nevertheless, in your model, (x < y) must
yield 0.

And just what do you mean by "the value in x (or y)"? Are you
implying that x has just a single value?

Incidentally, in my article there were 5 paragraphs following what you
quoted, full of questions directed to you in an attempt to understand
your conception of C's floating-point model. You have snipped those
paragraphs without so much as a "[...]". My assumption is that you
were unable to answer my questions because your floating-point model
is actually inconsistent. Either that or you didn't bother to read
them.

No, I try to keep messages short and simple, and dealing with one
aspect at a time. This is getting far too long.

WRT the values in x (or y), those values are exact, but they do NOT
represent exactness. This is due to the fundamental nature of a
floating value, i.e. any exact value ALWAYS represents a range,
such as:

x * (1 - x_EPSILON) through x * (1 + x_EPSILON)

and you know nothing better about it without detailed analysis of
the code that set the value. This is one of the reasons using
equality relations in FP is unwise. x < y implies that ALL values
in the x range are less that ALL the values in the y range. The
structure of FP values is such that you can tell this immeditely if
the value stored in x is less than the value stored in y.

When you write out the value in x to umpty-ump digits you are doing
nothing useful. That is simply one of the many possible values
that that x can represent (I am mixing up the usage of value).

If you start taking differences in values you will find that the
resultant errors can be overwhelmingly due to the original ranges,
which were not errors, but simply unknown details. Differences of
nearly equal quanities, etc. If you use products and quotients,
the resultant error will probably not become excessive. These
differences in resultant errors are what limits inversion of
matrices, etc.

I don't do well at reexplaining things. I know I tend to harp on
the same aspect, in the hope that slightly different words will get
the idea across. Often it doesn't. Thus I have long since given
up the idea of being a teacher.

 

[CBFalconer]

Try thinking about it. For example, if I write:

float x = 1.0/3.0;

you are claiming the value stored is 1/3.
[...]

I don't believe anybody made such a claim. (It's possible only
of FLT_RADIX is a multiple of 3.)

Then you can't claim you know what is stored there.

[...]

I can know what is stored in x by examining the value of x. On
my system, it's 0.3333333432674407958984375. On some other
system, it may be slightly different.

Given the above declaration, the C standard doesn't tell me
exactly *which* exact floating-point value is stored in x;
that's implementation-defined, and even the implementation's
definition might not be enough to determine it. But it does
tell me that *some* exact floating-point value is stored in x,
since (barring trap representations) that's the only thing
that possibly *can* be stored in x.

See my earlier reply to you, a few minutes ago.

 

[CBFalconer]

.... snip ....

Here's the full quote:
| Try thinking about it. For example, if I write:
|
| float x = 1.0/3.0;
|
| you are claiming the value stored is 1/3. But if you examine the
| value stored you will not find that. You will find something in
| the range x*(1+FLT_EPSILON) and x*(1-FLT_EPSILON). That range is
| system dependant. We won't move that range very much if we
| substitute the desired value, 1/3, for the value found in x.
| Things may be worse, but they won't be better.

Assuming that the "." in "1/3." was intended to be a period, not a
decimal point, your statement "you are claiming the value stored is
1/3" is simply false. Nothing that follows justifies or explains it.

Assuming binary floating-point, the mathematical value 1/3 is not
representable in type float. The value stored is not 1/3; it's a
representable value close to the real number 1/3.

Exactly. All you know is that the content of x represents a value
in the range x*(1+FLT_EPSILON) through x*(1-FLT_EPSILON).

 

[CBFalconer]

(e-mail address removed) writes:
.... snip ...

But that is error analysis using interval arithmetic.


And those ranges are useless in error analysis because after
only a single calculation the range will already be too small,
as in z = x + y.

I think you are mixing up errors due to the representation, and
errors due to the processing. The representation errors are
known. They may or may not be significant. The processing you
control.