Testing for very small doubles with DBL_EPSILON and _isnan()

[James Kuyper]

Expected for that particular value pair. What does it get for
pow(-1e-7, -1.5)? OP may be getting this pair or not.
?

As DIk said "... if a was negative ...". So do you now agree with what
he said? Would you still expect "a NaN and an overflow indication" even
when a is not negative?

 

[Dik T. Winter]

> On Sat, 6 Jun 2009 00:59:35 GMT, "Dik T. Winter" <[email protected]>
> wrote: ....
>
> I took the definition right off of K&R II, page 258.

Yes, it is there indeed, so this needs a place in the errata...

> Your definition is from (or consistent with) n1256 but that document
> gives five possible rounding methods, at least two of which render
> your assertion false.

Eh? I specifically said "with standard rounding", meaning IEEE standard
round to nearest.

 

[JosephKK]

JosephKK wrote:
...

DBL_EPSILON Is not a mathematical abstraction, it is a macro defined in
the C standard library, with a value whose meaning is specified by the C
standard.

When you wrote "Not really.", you were objecting to a true statement
about DBL_EPSILON, not a statement about "the real properties of real
hardware".

When you wrote the sentence which starts "It specifies ...", the subject
of that sentence was "It", which in context had to be a reference to the
use of DBL_EPSILON in the message you were responding to. Again, what
you said about "it" may have been true with regard to some entity
related to "the real properties of real hardware", but you didn't
identify which entity that is; it certainly is not DBL_EPSILON.


I trust DBL_EPSILON to be the difference between 1.0 and the next number
after 1.0 that is representable as a double. Combining that fact with
the definition of C's floating point model given in 5.2.4.2.2p1-2, and
the value of FLT_RADIX, I can calculate the spacing between x and the
next representable number for any value of x that satisfies (DBL_MIN <=
x && x < DBL_MAX) || (-DBL_MAX <= x && x < -DBL_MIN).

Thinking that DBL_EPSILON means anything other than that would be
"overextending" the definition (as CBFalconer so spectacularly
demonstrated in a recent thread), but that's not because of trusting the
programming language standards too much, but because of misunderstanding
what they say.

The C standard says what it says. It defines a few thing in certain
ways. Those definitions are often somewhat at odds with how the
hardware operates and the way that IEEE-754 specifies that they
operate. Thus when the hardware acts differently than what the C
standard specifies, it is not the hardware that is defective.
;
Never mind we are talking past each other.

 

[JosephKK]

It would be better if you read the standard before pronouncing things like
this.

Well fudge. The standard defines it the way that it does. That
definition is in conflict with the way the IEEE-745 compliant hardware
operates. And i am a hardware type.

 

[JosephKK]

As DIk said "... if a was negative ...". So do you now agree with what
he said? Would you still expect "a NaN and an overflow indication" even
when a is not negative?

Some 25 years ago i studied the IEEE-754 specification pretty
thoroughly, it left me with a deep sense of how the hardware would
operate. This came after deep understanding of two other
implementations all the way down to the gate level.

 

[Tim Prince]

Maybe I'm dense. pow(2,-53) is DBL_EPSILON/2. Adding pow(2,-105) results
in something that looks like this:

00111100 10100000 00000000 00000000 00000000 00000000 00000000 00000001
Exp = 970 (-52)
111 11001100
Man = .10000 00000000 00000000 00000000 00000000 00000000 00000001
1.1102230246251568e-16

"standard rounding" would be round to nearest even. 1 + DBL_EPSILON/2
thus would round to 1.0, while 1 + DBL_EPSILON/2*(1+DBL_EPSILON) would
round to 1 + DBL_EPSILON.

 

[James Kuyper]

JosephKK wrote:
....

The C standard says what it says. It defines a few thing in certain
ways. Those definitions are often somewhat at odds with how the
hardware operates and the way that IEEE-754 specifies that they
operate. Thus when the hardware acts differently than what the C
standard specifies, it is not the hardware that is defective.
;
Never mind we are talking past each other.

I'd prefer not to be talking past each other. I'm aware that there are
inconsistencies of the kind you mention in some other parts of the
standard; in particular, I've heard that some of the specifications in
Annex F about how signed zeros and signed infinities are supposed to be
handled don't jibe with the IEEE-754.

But is there anything that the standard says about DBL_EPSILON
specifically which is "somewhat at odds with the hardware operates and
the way that IEEE-754 specifies that they operate"? If so, how? I'm at a
loss to imagine how it could be. As far as I know, all floating point
representations have a unique next-representable value after 1.0, and
that would seem to be the only requirement for the definition of
DBL_EPSILON to be meaningful.

 

[James Kuyper]

On Mon, 08 Jun 2009 11:05:14 GMT, James Kuyper


Some 25 years ago i studied the IEEE-754 specification pretty
thoroughly, it left me with a deep sense of how the hardware would
operate. This came after deep understanding of two other
implementations all the way down to the gate level.

I asked a question with a yes or no answer; I can't translate what
you've said into either "yes" or "no", nor can I translate it into
anything resembling a claim that the restriction to "yes" and "no" is
inappropriate.

In other words, I don't see how your response can be interpreted as an
answer to my question.

 

[Keith Thompson]

Well fudge. The standard defines it the way that it does. That
definition is in conflict with the way the IEEE-745 compliant hardware
operates. And i am a hardware type.

How so? The definition of DBL_EPSILON isn't about the way the
hardware operates; it's an unambiguous definition of a number. Just
what conflict do you see?

 

[Dik T. Winter]

> Dik T. Winter wrote: ....
>
> Maybe I'm dense. pow(2,-53) is DBL_EPSILON/2. Adding pow(2,-105) results
> in something that looks like this:
Right.

> but adding it to 1.0 it is still too small to make a difference. I don't
> see that rounding gets you where you want to go.

No, adding it to 1.0 results in a number different from 1.0. Adding
pow(2, -53) to 1.0 results in 1.0, double IEEE, round to nearest.

> Standard says 2^-52.

That is the definition of DBL_EPSILON.

 

[Dik T. Winter]

> The C standard says what it says. It defines a few thing in certain
> ways. Those definitions are often somewhat at odds with how the
> hardware operates and the way that IEEE-754 specifies that they
> operate.

In what way is the hardware operation different from the C standard definition
of DBL_EPSILON? Or do you state that in hardware DBL_EPSILON (that is the
difference between 1.0 and the next higher number) is not equal to the
difference between 1.0 and the next higher number? Why not?

 

[Dik T. Winter]

> On Mon, 8 Jun 2009 10:19:57 GMT, "Dik T. Winter" <[email protected]>
> wrote: ....
>
> Well fudge. The standard defines it the way that it does. That
> definition is in conflict with the way the IEEE-745 compliant hardware
> operates. And i am a hardware type.

In what way does it conflict with how the hardware operates? I have no
idea about what you do mean here. DBL_EPSILON is defined in a certain,
precise, way. I do not see how the hardware comes into it.

 

[JosephKK]

In what way does it conflict with how the hardware operates? I have no
idea about what you do mean here. DBL_EPSILON is defined in a certain,
precise, way. I do not see how the hardware comes into it.

I did a study of compliance of 8087 versus the Cyrix and Weitek knock
offs about 25 years ago. The standard has not changed much since. If
you really, really want to know get a copy and study it.

Example: with 53 bits of mantissa, DBL_EPSILON would be closer to
1E-19 rather than 1E-9

 

[JosephKK]

I asked a question with a yes or no answer; I can't translate what
you've said into either "yes" or "no", nor can I translate it into
anything resembling a claim that the restriction to "yes" and "no" is
inappropriate.

In other words, I don't see how your response can be interpreted as an
answer to my question.

Please read it as my background tells me what to expect, based on the
described values. If that conflicts with what you expect, please test
it and tell us all what happened; oops that would require assembler
programming instead of "C". Then again you may be well able to
develop those programs anyway. It is up to you to do what please.

 

[JosephKK]

No, adding it to 1.0 results in a number different from 1.0. Adding
pow(2, -53) to 1.0 results in 1.0, double IEEE, round to nearest.



That is the definition of DBL_EPSILON.

But that is not what the "C" standard says, just for starters. BTW
are quoting from the 1985 IEEE-754 or some other standard?

Plenty of information is readily available for a web search or three.

 

[Flash Gordon]

I did a study of compliance of 8087 versus the Cyrix and Weitek knock
offs about 25 years ago. The standard has not changed much since. If
you really, really want to know get a copy and study it.

Example: with 53 bits of mantissa, DBL_EPSILON would be closer to
1E-19 rather than 1E-9

So? The C standard says that if the correct value of DBL_EPSILON is
1E-19 then that is what the implementation has to define it as. The
specific value mentioned in the standard is a limit on what it is
allowed to be, not the specific value to be used in every implementation.

To be clear, if the difference between 1.0 and the next representable
number in the floating point format used by the hardware is 1E-19, then
the C standard *requires* that DBL_EPSILON be defined as 1E-19 (this
specific exact value would require base 10 floating point rather than
base 2, which is allowed by the standard).

The C++ standard is, I believe, the same in this respect.

 

[user923005]

I did a study of compliance of 8087 versus the Cyrix and Weitek knock
offs about 25 years ago.  The standard has not changed much since.  If
you really, really want to know get a copy and study it.

Telling Dik Winter to study an applied math paper is like telling
George Marsaglia to read you paper on PRNGs or Donald Knuth to read
your article about the art of computer programming.

Example: with 53 bits of mantissa, DBL_EPSILON would be closer to
1E-19 rather than 1E-9

I would say that 1E-19 is off by about a factor of ~ 2220 unless you
know something I don't know.
{My calculation would be: eps is 2.0/pow(2.0,53.0) for a 53 bit
mantissa.}

 

[James Kuyper]

Please read it as my background tells me what to expect, based on the
described values. If that conflicts with what you expect, please test
it and tell us all what happened; oops that would require assembler
programming instead of "C". Then again you may be well able to
develop those programs anyway. It is up to you to do what please.

My question was about what your expectations were; I can't get the
answer to that question by writing a test program; not unless that test
program is hooked up to mind-reading hardware.

Also, this was a question about the behavior of C code, not about
assembler code, so what point would there be to writing assembler code
to determine the answer?

If the floating point hardware implements a near-equivalent of the
pow(a,n) function, that differs from the specification of pow() by
returning a NaN and an overflow indicator when a is a small positive
number, and b is a sufficiently large negative number, then the
implementation of the pow() function is required to contain whatever
additional code is needed to ensure that the behavior of the pow()
function itself meets the requirements of 7.12.1p4:

"If a floating result overflows and default rounding is in effect, or if
the mathematical result is an exact infinity from finite arguments (for
example log(0.0)), then the function returns the value of the macro
HUGE_VAL,HUGE_VALF, or HUGE_VALL according to the return type, with the
same sign as the correct value of the function; if the integer
expression math_errhandling & MATH_ERRNO is nonzero, the integer
expression errno acquires the value ERANGE; if the integer expression
math_errhandling & MATH_ERREXCEPT is nonzero, the ‘‘divide-by-zero’’
floating-point exception is raised if the mathematical result is an
exact infinity and the ‘‘overflow’’ floating-point exception is raised
otherwise."

Note: I'm pretty sure that the parts involving math_errhandling are new
to C99; I'm not sure what the corresponding wording in the C90 standard was.

The required return value in this case, HUGE_VAL, must be positive
(7.12p3), and therefore cannot be a NaN, and must (if __STDC_IEC_559__
is predefined) be positive infinity (F9p1).

There's nothing particularly unusual about a discrepancy between the
behavior of the floating point hardware and the behavior of the
corresponding standard library function. That's precisely the point of
standardizing the library, so people can write code that relies upon the
standard's guarantees, without having to worry about the details of
different hardware implementations.

 

[James Kuyper]

....

But that is not what the "C" standard says, just for starters.

Citation, please? And don't quote the upper limit of 1E-9 again; 2^-52
easily satisfies that limit. If the smallest representable value greater
than 1.0 is 1.0+pow(2,-52), then the C standard says that DBL_EPSILON
must be pow(2,-52).

 

[Dik T. Winter]

> On Tue, 9 Jun 2009 12:44:39 GMT, "Dik T. Winter" <[email protected]>
> wrote: ....
>
> I did a study of compliance of 8087 versus the Cyrix and Weitek knock
> offs about 25 years ago. The standard has not changed much since.

The standard has not changed at all, but the 8087 (and 80187) is not standard
compliant (it was based on a preliminary version of it). As far as I know,
Cyrix and Weitek were standard compliant. Yhe 80287 was the first standard
compliant Intel FP processor.

> If
> you really, really want to know get a copy and study it.

Do you think I do not know that standard?

> Example: with 53 bits of mantissa, DBL_EPSILON would be closer to
> 1E-19 rather than 1E-9

Where is written that DBL_EPSILON is something like 1E-9? With IEEE,
DBL_EPSILON equals 2^(-52).