Testing for very small doubles with DBL_EPSILON and _isnan()

[JosephKK]

No, that would result in an infinity, not in a NAN. The only way that could
yield a NAN would be if a was negative (and with a well-designed library,
b non-integral). In all other cases a number or an infinity should be
returned (with the now common IEEE standard arithmetic).

So I suspect your a has become negative.

As much as i don't like to argue with you i expect a NaN and an
overflow indication (which may be going unhandled). But then i am
coming from the hardware side.

 

[Martin Eisenberg]

On Thu, 04 Jun 2009 19:30:04 -0700, Keith Thompson


Not really. It specifies the guaranteed accuracy of
mathematical operations, not the resolution.

Yes it does. C++98 final draft, 18.2.1.2/20 says:

| static T epsilon() throw();
|
| Machine epsilon: the difference between 1 and the least value
| greater than 1 that is representable. 8)
|
| 8) Equivalent to FLT_EPSILON, DBL_EPSILON, LDBL_EPSILON.

The C++ standard does not define accuracy, but IEEE754 mandates that
the relative error of basic operations not exceed 1 ulp which roughly
means the interval [-2*epsilon, epsilon].


Martin

 

[JosephKK]

I'd suggest printing the values of all these variables (r2, m_alpha, and
all the deltas) before the computation of `output'. That will probably
bring enlightenment.


As mentioned earlier, my guess would be that r2 comes out negative, so
you could test for that. It's also conceivable that one of the deltas
got set to a NaN somehow, depending on what initializeDeltaXYZ does.

I do not see any way for r2 to ever come out negative, deltaX is a
different matter.

 

[JosephKK]

You misunderstand DBL_EPSILON. It is the difference between 1.0 and the
next larger number. With standard rounding:
1.0 + 3 * DBL_EPSILON / 4 > 1.0

No. See my previous explanation. It represents the accuracy ratio.

 

[Martin Eisenberg]

The C++ standard does not define accuracy, but IEEE754 mandates
that the relative error of basic operations not exceed 1 ulp
which roughly means the interval [-2*epsilon, epsilon].

Correction: [-1,1]*epsilon.


Martin

 

[Nate Eldredge]

As much as i don't like to argue with you i expect a NaN and an
overflow indication (which may be going unhandled). But then i am
coming from the hardware side.

On my system (amd64, gcc, FreeBSD), print("%f\n", pow(0.0, -1.0)) prints
"inf".

 

[Nate Eldredge]

I do not see any way for r2 to ever come out negative, deltaX is a
different matter.

Oops, I didn't read it carefully. r2 is a sum of squares.

 

[James Kuyper]

....

Not really. ...

Yes, really. The C standard defines it as: "the difference between 1 and
the least value greater than 1 that is representable in the given
floating point type ..." (5.2.4.2.2p11). Since this is cross-posted to
C++, I'll also mention that the C++ standard incorporates the C standard
library by reference (1.2p2) and specifically mentions both <float.h>
(D.5p1) and DBL_EPSILON (C.2p4) as features of that library which are so
incorporated, and specifies that for <cfloat>　â€The contents are the

> ... It specifies the guaranteed accuracy of mathematical
operations, not the resolution The difference can be pretty subtle,
but consider: if you are dealing with largish quantities, say
Avogadro's number, and you want to know the number of atoms delta
resolution that can be represented it would still be many trillions.
But the accuracy ratio scales directly with the size of the numbers
represented, and this is but one of the contributing components to a
property called numerical instability.

The definition of DBL_EPSILON is cited above. It's quite simple, and
makes no mention of any of those issues. Those issues are certainly
related to it, but not in any way that justifies saying "Not really.".

 

[James Kuyper]

incorporated, and specifies that for <cfloat>　â€ÂThe contents are the
same as the Standard C library header <float.h>â€Â． (18.2.2p4).

Sorry about that - but I've finally figured out what's going wrong. My
computer has software installed that lets me easily insert characters in
encodings and fonts designed for showing both Chinese and English text.d
It switches modes automatically when I hit Shift-Space, or Ctrl-Space,
key combinations that I've been hitting by accident far too frequently
lately.

That was supposed to say:

 

[Richard Bos]

I'd suggest printing the values of all these variables (r2, m_alpha, and
all the deltas) before the computation of `output'. That will probably
bring enlightenment.

I agree.

As mentioned earlier, my guess would be that r2 comes out negative,

That seems unlikely; floating point inaccuracies are one thing, but
changing the sign of a sum of (non-complex) squares to negative is
another game.

so you could test for that. It's also conceivable that one of the deltas
got set to a NaN somehow, depending on what initializeDeltaXYZ does.

Sounds more likely to me. Or that m_alpha is.

Richard

 

[Keith Thompson]

Isn't it true for *all* negative numbers?

Of course (and for 0.0), but the intended point of the comparison was
to determine whether x is "really small".

 

[JosephKK]

On Thu, 04 Jun 2009 19:30:04 -0700, Keith Thompson


Not really. It specifies the guaranteed accuracy of
mathematical operations, not the resolution.

Yes it does. C++98 final draft, 18.2.1.2/20 says:

| static T epsilon() throw();
|
| Machine epsilon: the difference between 1 and the least value
| greater than 1 that is representable. 8)
|
| 8) Equivalent to FLT_EPSILON, DBL_EPSILON, LDBL_EPSILON.

The C++ standard does not define accuracy, but IEEE754 mandates that
the relative error of basic operations not exceed 1 ulp which roughly
means the interval [-2*epsilon, epsilon].


Martin

I was talking in clc, was not mindful of crossposting. Besides being
more of a hardware type i care more about IEEE-754 terminology.
j

 

[JosephKK]

Sorry about that - but I've finally figured out what's going wrong. My
computer has software installed that lets me easily insert characters in
encodings and fonts designed for showing both Chinese and English text.d
It switches modes automatically when I hit Shift-Space, or Ctrl-Space,
key combinations that I've been hitting by accident far too frequently
lately.

That was supposed to say:

Well shit. The standards say what the standards say. I am corrected
on what they say. It does not alter the real properties of the real
hardware, nor the mathematical properties thereof. Overextend your
trust in the programming language standards at your own risk.

 

[JosephKK]

Oops, I didn't read it carefully. r2 is a sum of squares.

Just the same i can't figure out the issue or the underlying problem
well enough to help OP much. I wonder if OP is trying to "white room"
a well known mathematical algorithm.
?

 

[JosephKK]

On my system (amd64, gcc, FreeBSD), print("%f\n", pow(0.0, -1.0)) prints
"inf".

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

 

[JosephKK]

No. See my previous explanation. It represents the accuracy ratio.

Yeah, i see what the standards say. I am talking about how the
hardware and the math works.

 

[Dik T. Winter]

> On Thu, 04 Jun 2009 19:30:04 -0700, Keith Thompson <[email protected]>
> wrote: ....
>
> Not really. It specifies the guaranteed accuracy of mathematical
> operations, not the resolution.

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

 

[Dik T. Winter]

> On Sat, 6 Jun 2009 00:59:35 GMT, "Dik T. Winter" <[email protected]>
> wrote:
>

>
> No. See my previous explanation. It represents the accuracy ratio.

Have you looked in the standard how it is defined?

 

[Dik T. Winter]

> On Sat, 06 Jun 2009 15:39:34 -0700,

>
> Yeah, i see what the standards say. I am talking about how the
> hardware and the math works.

When you use DBL_EPSILON you should be aware what the definition is
according to the standard. The hardware does not define that macro.

On the other hand, if you are looking for the smallest positive value that,
when added to 1.0 yields a value different from 1.0, for double IEEE that
would depend on the rounding mode. With standard rounding it is:
pow(2, -53) + pow(2, -105)

 

[James Kuyper]

JosephKK wrote:
....

Well shit. The standards say what the standards say. I am corrected
on what they say. It does not alter the real properties of the real
hardware, nor the mathematical properties thereof. ...

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.

... Overextend your
trust in the programming language standards at your own risk.

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.