0

In my program I calculate two values (doubles) x and y and then find the number 1/(x^2 + y^2). My problem arises when x or y are so close to zero that the fraction gives a NaN. In reality, the fraction should give zero when both variables approach zero. I tried to use a comparison x==0 || y==0, but that doesn't work because they are doubles.

Is there a computationally effective method to take this into account?

#include <iostream>

using namespace std;

val(double x, double y)
{
    return 1/(x*x + y*y);
}

int main()
{
    //determine x and y....
    val(x, y)
    return 0;
}
doctorlove
  • 18,872
  • 2
  • 46
  • 62
BillyJean
  • 1,537
  • 1
  • 22
  • 39
  • @Rob013 doesn't work, because sometimes `x=...e-312`, that is what I get when I use `cout << x << endl;`. So `x==0.0` outputs `0` – BillyJean Jul 13 '13 at 09:47
  • I don't get why "in reality, the fraction should give zero when both variables approach zero." When a constant divides an infinitesimally small number, it should approach infinity, not 0. – cheeyos Jul 13 '13 at 09:50
  • @cheeyos because that is how my function is constructed, it is a piecewise function – BillyJean Jul 13 '13 at 09:52
  • Why do you want it zero in this case? Mathematically this fraction does not go to zero when x and y go towards 0. Consider `x=y=.1`, then `frac=50`, `x=y=.01` then `frac=5000` and so on, so it does clearly *not* go to zero. – KillianDS Jul 13 '13 at 09:58
  • @KillianDS I suppose BillyJean is using a "volcano function" with a `f(x, y) = 0 if x=y=0` condition. It should be similar to this: http://www.wolframalpha.com/input/?i=1%2F%28x^2%2By^2%29 . Btw, I suppose you should limit the height of your volcano and use nio's solution, even if it's not valid for "every small double". – Cob013 Jul 13 '13 at 10:09
  • This question is defective, and a correct answer cannot be given until it is repaired. When e is tiny, 1/e is huge or, if overflow occurs, infinity. It is not NaN. So the question states a false premise. If it is corrected by explaining that some other function is being computed, e.g., `x/(x^2+y^2)`, then it could return a NaN for small x and y, and suggestions to avoid that could be provided. If it is corrected by explaining that NaN is not being returned, or that some other bug led to this behavior, then there is no need to alter the evaluation of this expression. – Eric Postpischil Jul 13 '13 at 13:28

3 Answers3

1

If you want to compare a double value with zero,you should do just like this:

bool IsZero(double dest)
{
if(dest > -0.0000001 && dest < 0.0000001)
return true;
else
return false;
}
minicaptain
  • 1,196
  • 9
  • 16
  • 1
    you know you can reduce this to the much easier oneliner `return dest > -0.0000001 && dest < 0.0000001` right? – KillianDS Jul 13 '13 at 10:05
  • This answer proposes that an error be corrected by compounding it with more error. That is, we are confronted with a situation in which `dest` is near zero but, due to rounding or other errors, we are unsure whether it is zero. This answer proposes that we treat it as zero, thus possibly increasing previous errors. The correct behavior depends on knowing more about the circumstances. Since the statements in the question are inconsistent (asserting that NaN is returned when IEEE 754 calls for infinity to be returned), we cannot answer the question with the current information. – Eric Postpischil Jul 13 '13 at 13:31
0

I don't know why you're getting NaN, unless x or y is already NaN, or the calculation you're doing is not the one you've shown, or you're using some non-IEEE floating-point implementation.

With IEEE floating-point, a small floating-point value squared is always representable, either as a smaller non-zero value or as 0. In any event, don't get distracted by the computations in the denominator. The underlying problem is that at some point, 1/d, even when d is not 0, becomes too large to represent as a double and you get an infinite result. If you want to turn infinities into zeros the simplest approach is to check for infinity and replace it with 0:

#include <math.h>
double res = // whatever
if (isinf(res))
    res = 0;

This assumes C99 or C++11.

Pete Becker
  • 74,985
  • 8
  • 76
  • 165
-1
#include <float.h>

// could also be your custom consant
// see also http://stackoverflow.com/a/77735/213682
// also this http://stackoverflow.com/a/2729750/213682
const double my_epsilon = DBL_MIN;

double val(double x, double y)
{
    double denom = x*x + y*y;
    if (denom < my_epsilon) {
        return 0.0; // or some other value indicating an error
    }
    return 1.0 / denom;
}

On my machine, 1.0 / DBL_MIN gives a mighty number:

44942328371557897693232629769725618340449424473
55766431835752028943316895137524078317711933060
18840052800284699678483394146974422036041556232
11857659868531094441973356216371319075554900311
52352986327073802125144220953767058561572036847
82776352068092908376276711465745599868114846199
29076208839082406056034304.000000
Alexei Sholik
  • 7,287
  • 2
  • 31
  • 41