Float comparison

[BartC]

Several years ago, on a dark and stormy Saturday, I came up with this as a
C exercise.

One survey foot = 0.3048006096 meters.
One statute mile = 1.6093472187 kilometers.
One inch = 0.0254 meters (exactly).

If you did this in C, it probably won't be exact. You need to apply CBF's
Theorem.

 

[Phil Carmody]

You need better software.

In what way does better software change the contents of a file?

Phil

 

[Phil Carmody]

If you did this in C, it probably won't be exact. You need to apply
CBF's Theorem.

:-D

Possibly the post with highest utility in recent times.

Phil

 

[Richard Tobin]

It's a jolly good job C doesn't have an exponentiation
operator, or you could have some real fun with 0. ^^ 0. ,
as neither its sign nor its finiteness would be known,
despite the fact that it's clearly 1.

Why clearly?

0 ^^ 0 is clearly 1 where the zeros are cardinals, since there
is exactly one mapping from the empty set to itself. But
writing "0." strongly suggests you are considering some other
variety of numbers.

-- Richard

 

[Ben Bacarisse]

In what way does better software change the contents of a file?

In what way did I say it might?

 

[Dik T. Winter]

> Dik T. Winter wrote: ....
>
> Can you tell the rounding in force? The following is (float)1.0/3.0
>
> 00111110 10101010 10101010 10101011
> Presumably 10101010.1010 was rounded up?

Clearly round to nearest. As the mathematical result is not halfway two
representable numbers, there is no additional problem.

 

[CBFalconer]

So what happened to the ranges of the operands?

In your model, each of x and y represents a range of real values.
For concreteness, let's say the ranges are 0.99 .. 1.01 (that's
obviously oversimplified, but it should suffice to make the point).
Then the possible range of the result of the multiplication should
be 0.9801 .. 1.0201, shouldn't it? This is just over twice the
range of the operands.

But you say that z represents the same range as x and y. Why?

Because that 'range' describes the range of real values that may
have been stored in the fp-object. We know the fp-multiplication
used the fp-exact values, so the result was 1.0. The range says
that is a possible value to have been originally stored in z.

If we eliminate consideration of the programming that set x and y
(the x = 1.0; y = 1.0 now we still know that the multiplication
used the exact fp-values of x and y, but we don't know exactly what
that represented. Now the 'range' of the result is unchanged - it
depends only on the stored value. But the possible error has to
consider the possible errors in the input to the multiplication,
i.e. the 'range' of those input values. The 'range' is not an
error measure - it is a limit to the precision of storage.

 

[CBFalconer]

what I call CBF's theorum:
"the range represented by z = x * y consists of the extreme range
values of the values calculated by multiplying the extremes of the
ranges x and y"

No. The range of z x*y depends on the fp-object value of z
alone. This is not an error measure. It is a minimum error
measure.

 

[CBFalconer]

.... snip ...

It's a jolly good job C doesn't have an exponentiation
operator, or you could have some real fun with 0. ^^ 0. ,
as neither its sign nor its finiteness would be known,
despite the fact that it's clearly 1.

No, it's unknown. Any power of 0 is 0, except the 0th power. The
0th power of anything (except 0) is 1. If the fp-systems has NaNs
it should return one for 0 raised to the 0th power.

 

[Dik T. Winter]

> Several years ago, on a dark and stormy Saturday, I came up with this as
> a C exercise.
>
> One survey foot = 0.3048006096 meters.
> One inch = 0.0254 meters (exactly).

From 1866 in the US 1 meter = 39.37 inches exactly. The above value for the
icnh became standard in the US in 1959 and in the UK in 1963. The older
standard gives your value of the survey foot.

> One US gallon = 3.785411784 liters (exactly).

Since 1707 a US gallon was standardized at 231 cubic inches. Depending on the
value of the inch it gives you the exact value in liters. In 1995 the UK
defined the gallon as 4.54609 liters exactly.

> One pound = 453.59237 grams (exactly).

From 1866 in the US, 1000 grams = 2.2046 pounds exactly, changed in 1893 to
2.20462 pounds and in 1894 to 2.20462234 pounds. Since 1959 the value above
in the US, and the same value in the UK from 1963.

For the US see: <http://physics.nist.gov/Pubs/SP447/contents.html>.

So we see that from 1866 the US system was based on the metric system. It was
only in 1995 that the UK based all of their system on it.

 

[CBFalconer]

.... snip ...

Let me clarify something.

As a floating constant of type double, appearing in a C program,
there's no effective difference between
0.333333333333333314829616256247390992939472198486328125
and
0.3333333333333333
Both evaluate (on my system) to exactly the same double value.

As real numbers, though, they are quite different. They're
close, but they differ by a known amount.

And thus, when we are evaluating the fp-object-value created, there
is no difference whatsoever. If we find a difference, we are
evaluating things not in that fp-object, i.e. the programming
behind its existance.

There is nothing wrong with saying 'Those are real values. Here is
what the fp-system does to them.'. We can evaluate the things that
are identical to the fp-system quickly by evaluating the 'range' of
the fp-object-value. We can easily evaluate that 'range' by
computing xmax and xmin for x=1/3, remembering that those values
are real and computed for their effect on the fp-system only.

What I have read from
<http://docs.sun.com/source/806-3568/ncg_goldberg.html> so far does
not disagree with me in the slightest. Goldbergs paper is more
thorough than I have been so so far. My writings here have been
gradually refined. His paper did the refining before publication.

 

[CBFalconer]

.... snip ...


I thought we had established that you were using EPSILON as a
general term for FLT_EPSILON, DBL_EPSILON, and LDB_EPSILON.

Now it appears that "epsilon" (lower case), "EPSILON" (upper
case), and "DBL_EPSILON" (the constant defined by the C
standard) are *three* different things.

I have directly asked you about this, and you have yet to give a
coherent answer; it just becomes more complicated every time you
try to define it.

I have been using lower case to define real arithmetic values.
EPSILON defines a defined value, as in the C standard. I think I
have been consistent in this action.

 

[CBFalconer]

.... snip ...

So in the context of C, does the word EPSILON refer to any of
FLT_EPSILON, DBL_EPSILON, or LDBL_EPSILON? Yes or no?

It refers to whichever, if any, is applicable to the fp-system
under discussion.

Upthread, you wrote:

Think about WHY EPSILON changes at the power of two. It has to
do with the rounding in computing the real value x*(1+EPSILON)
when that expression is handed to the fp-system.

Would you care to explain that remarkable statement?

What's remarkable? When you compute and use that expression in the
fp-system, the results are exactly equal because that EPSILON is
the smallest value that can get noticed. When x doubles it size,
so does the resultant epsilon, and is suddently is big enough to be
noticed.

Compute that x*(1+EPSILON), using the appropriate float.h value of
EPSILON, for various x's from 1.0 to something less that 2.0. Note
that you get the same result from x+EPSILON in that range. (Here
range means from x=1.0 to less that 2.0).

That EPSILON should correspond to one half the weight of the least
significant bit in the significand of 1.0.

 

[CBFalconer]

.... snip ...


Why complicate it like that? The "least value greater than 1"
*is* the successor value.

When you talk about inserting real values into the FP system,
you need to consider rounding errors. Rounding errors are not
relevant to the definition of the *_EPSILON constants.

Not errors, just rounding rules. Yes they are. That is why I
interpreted it as I did.

Think of this - you have a floating value x, and you add something
that is identical (d), except all significand bits are zero except
the one one position to the right of the significand. When that is
normalized it will change the exponent, and have a single most
significant significand bit, which will then be replaced by the
sign bit. I.e. the significand has become a zero after
normalization.

When this is added to x you expect it to get rounded, and x will
become one larger in its least significant bit. Otherwise your
system is truncating, and this leads to large errors.

 

[CBFalconer]

.... snip ...


How can b1-p (badly formatted but presumably you know what it
means) be interpreted to mean anything but one thing? For you
example system it is 1/8.

That is the representable form of the sum, not the input.

 

[CBFalconer]

.... snip ...

But that you are a designer of floating-point hardware does not
make you knowledgable in the handling of floating-point in
numerical analysis. Moreover, it does not show that you did make
the best possible. One of the most atrocious examples of badly
designed floating-point hardware can be found on the Cray-1 (it
could get the multiplication wrong by 4 "ulp"s, and it did not
do full division).

Hardware AND software. Only the first was purely hardware.
Believe it or not, I did learn something in 45 years of building
arithmetic systems. I didn't end up with the Crays problems.

 

[CBFalconer]

Apparently you know better than an excellent numerical
mathematician as Kahan, who basicly designed the IEEE system and
was involved in the design of the 8087 (which implemented a
preliminary form of the standard).

Why not? I presented the reasons. You cited somebody as God.

 

[Kenny McCormack]

a bunch of interesting, but alas (and more's the pity) off-topic stuff:
....

From 1866 in the US, 1000 grams = 2.2046 pounds exactly, changed in 1893 to
2.20462 pounds and in 1894 to 2.20462234 pounds. Since 1959 the value above
in the US, and the same value in the UK from 1963.

None of which is mentioned in the C standards documents.

Where is Kiki when you need him?

 

[CBFalconer]

CBFalconer wrote:
.... snip ...


Hmm. Se you think you know better than the IEEE committees
responsible for writing their floating point standards? I
somehow think that they might know just a little bit more than
you about floating point arithmetic.

I may well. I am at least willing to defend things.

Do you at least accept that your rounding mode is NOT what is
usually done in this day and age?

I don't know. You are probably correct. That doesn't mean the
usual is right. If it was, there wouldn't be any choice of
methods, would there?

 

[CBFalconer]

So you clearly do not follow the rules as set down in the IEEE
standard from 1985. That the last rounds down rather than up
gives better statistics.

As I said in the part you quoted. However I later saw a reason to
suspect the accuracy of that, and published it. Your response was
a reference to some god or other.