# How to divide two doubles accurately

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.

-
When you say that -1 <= a/b <= 1, do you mean in theory or in practice? Might you already have roundoff errors in a and/or b before you do the division? –  Steve314 Oct 8 '09 at 0:55
Please provide the actual numbers you're actually using to get the actual value of 13 or 14. –  S.Lott Oct 8 '09 at 1:31
Are you dealing with subnormal numbers? (en.wikipedia.org/wiki/Subnormal_number) –  Jason S Oct 8 '09 at 2:55
If be can be arbitrarily small then a/b will be arbitrarily large. –  Bill Oct 8 '09 at 16:29

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.

-
absolutely correct. –  Jason S Oct 8 '09 at 2:51
I think there may be a problem with a and b before I even get to this step. I'm going to have to get back to you guys. –  ldog Oct 8 '09 at 3:36
Ok so if anyone cares, the actual case of my error happened before I got to the division step and the condition abs(a)<=abs(b) was being broken before division. –  ldog Oct 8 '09 at 21:57

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;
}
}
``````
-
If -1 <= a/b <= +1 is theoretically satisfied but wrong in practice due to rounding errors, the OP might want to return +1 or -1 (depending on the signs of a and b) instead of 0. –  sellibitze Oct 8 '09 at 8:59

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.

-

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)

-

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
``````
-

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?

-
That's not strictly true. The IEEE standard specifies a few different rounding modes, but your point is valid: division hardware (and software emulations) are accurate to within the precision available. The problem here is lack of precision of the underlying representation, not a bad algorithm. –  Andy Ross Oct 8 '09 at 2:36
True. I'm assuming that a user who is (apparently) not particularly savvy about floating-point has probably not changed the rounding mode from the default. The rounding error in a non-default rounding mode can be as large as 1ulp, which is still fantastically smaller error than the error that the questioner is reporting. –  Stephen Canon Oct 8 '09 at 2:48

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.

-
No, this is true of only the particular two numbers a and b. I'm getting inaccuracies because I'm dividing by a number that is much too small, smaller than machine precision. I need a way to check when I can guarantee the result us atmost a certain value from the true result. –  ldog Oct 8 '09 at 0:53
The only number smaller than machine precision is 0. If you are not dividing by 0, then you are not dividing by a number smaller than machine precision. –  Sam Harwell Oct 8 '09 at 1:04
I wonder if he's talking about subnormal numbers. –  Jason S Oct 8 '09 at 2:54
@280z28: yes thats right, the only number smaller is 0. –  ldog Oct 8 '09 at 16:39
{"found_math": false, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 0, "mathjax_display_tex": 0, "mathjax_asciimath": 0, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0, "katex": 0, "math-container": 0, "wp-katex-eq": 0, "align": 0, "equation": 0, "x-ck12": 0, "texerror": 0}