How to divide two doubles accurately

Question

I have two

double a, b;

I know that the following is true

-1 <= a/b <= 1

however b can be arbitrarily small. When I do this naively and just compute the value

 a/b

the condition specified above does not hold in some cases and I get values like much greater than 1 in absolute value (like 13 or 14.)

How can I ensure that when I divide by b, I get a value such that the condition mentioned above can be enforced. In the case where I can not guarantee this I am happy to set the computed value a/b to 0.

Answer 1

I am a little confused by your question, but if you are suggesting that -1 <= a/b <= 1 holds true for all real value a and b, this is definitely not the case. Consider:

1 / 0.5 = 2

If you want to check whether the division will be within [-1, 1], why not just do the division and then act however you want when it falls between -1 and 1.

Answer 2

What you need to enforce is abs(a)≤abs(b). If that condition holds, then -1≤a/b≤1, regardless of floating-point precision used. Your logic error is occurring before the division, since at the division point abs(a)>abs(b) which violates your a-priori requirement.

Answer 3

Division on an IEEE-754 system is a correctly rounded operation, which means that if no overflow or underflow occurs, the result will always be within 0.5 "ulp" of the mathematical "infinitely precise" result. In non-FP-nerd speak, this means that the result will always be within a factor of about 2^-53 of the exact answer. Since you know that the infinitely precise result is between -1 and 1, overflow cannot occur; underflow can, but that would result in numbers very, very near zero, not on the order of 13.

Either your condition does not actually hold, or you are on a system that does not have IEEE-754 arithmetic, or there is a bug in your code. Can you post the values of a and b that are generating this result, and the code that you are using to do the division and print the result?

Answer 4

If you really want to test out your theory that this is a precision issue (though the other answers make it seem like you can't possibly be getting the results you're getting) you should be using GMP. You already said you're "dividing by a number that is much too small, smaller than machine precision." GMP's C++ bindings have arithmetic operators for everything and are actually kind of pleasant to use compared to the C bindings, give them a try!

      mpf_class f(1.5);        // default precision
      mpf_class f(1.5, 500);   // 500 bits of precision (at least)
      double out = f.get_d();  // get a double out of an mpf_class

Answer 5

It's very unlikely you're actually triggering data loss due to underflow. While it is possible doubles have an incredible range and you're not likely to hit it.

I would think the problem lies somewhere before this. You've either got a logic bug or you are simply eating up the available precision with a bunch of operations. Be especially wary of additions and subtractions. 1E20 + 1 = 1E20.

If it's due to eating up the precision then you'll have to redesign your routine or resort to an arbitrary-precision library for your math. (Beware--SLOW)

Answer 6

The only situations I can think of where this would be true, are if a/b should be between -1 and 1, but is not, because both a and b are derived from calculations that make them very sensitive to floating-point errors.

For instance, if you try to calculate b = 2-2*cos(theta) and a = theta*theta for theta in radians, near zero, you will run into issues because the calculation of b involves the subtraction of numbers very close to 1, even though mathematically the answer is very close to 1.0 for small angles. (If I run the calculation of a/b for theta = 0.4e-7 using the Spidermonkey Javascript engine, I get 1.029 even though it really should be darned close to 1.0)

280Z28's answer is where you want to be, however without knowing more about how a and b are derived, and under what conditions you get undesirable results, it's hard to say what the best way is to ensure the calculations give meaningful numbers.

Answer 7

You surely have some problem before you reach this point. Regardless, you can elimate that problem here if you are prepared to accept a 0 result.

float Calculate(float a, float b)
{
    if (abs(a) > abs(b))
    {
        return 0;
    }
    else
    {
        return a/b;
    }
}