Rounding error detection

Question

I have two integers n and d. These can be exactly represented by double dn(n) and double dd(d). Is there a reliable way in C++ to check if

double result = dn/dd

contains a rounding error? If it was just an integer-division checking if (n/d) * d==n would work but doing that with double precision arithmetic could hide rounding errors.

Edit: Shortly after posting this it struck me that changing the rounding mode to round_down would make the (n/d)*d==n test work for double. But if there is a simpler solution, I'd still like to hear it.

• Geom: @Ron that's the point. I can't represent 1/10 but I can represent 2/4. The goal is to check if it is possible.
• gammatester: Test if d is a power of 2. If not and n is no multiple of d, then n/d is not exactly representable as double.
• Geom: @gammatester: 9/3 is exactly representable as integer and thus also as double.
• gammatester: Already edited.
• Pascal Cuoq: @gammatester d doesn't have to be either a power of two or a divisor of n, it can be the product of the two. Example: d=6, n=3.

Answer 1

If a hardware FMA is available, then, in most cases (cases where n is expected not to be small, per below), the fastest test may be:

#include <cmath>
...
double q = dn/dd;
if (std::fma(-q, dd, dn))
 std::cout << "Quotient was not exact.\n";

This can fail if nd−q•dd is so small it is rounded to zero, which occurs in round-to-nearest-ties-to-even mode if its magnitude is smaller than half the smallest representable positive value (commonly 2⁻¹⁰⁷⁴). That can happen only if dn itself is small. I expect I could calculate some bound on dn for that if desired, and, given that dn = n and n is an integer, that should not occur.

Ignoring the exponent bounds, a way to test the significands for divisibility is:

#include <cfloat>
#include <cmath>
...
int sink; // Needed for frexp argument but will be ignored.
double fn = std::ldexp(std::frexp(n, &sink), DBL_MANT_DIG);
double fd = std::frexp(d, &sink);
if (std::fmod(fn, fd))
 std::cout << "Quotient will not be exact.\n";

Given that n and d are integers that are exactly representable in the floating-point type, I think we could show their exponents cannot be such that the above test would fail. There are cases where n is a small integer and d is large (a value from 2¹⁰²³ to 2¹⁰²⁴−2⁹⁷², inclusive) that I need to think about.

Answer 2

If you ignore overflow and underflow (which you should be able to do unless the integer types representing d and n are very wide), then the (binary) floating-point division dn/dd is exact iff d is a divisor of n times a power of two.

An algorithm to check for this may look like:

assert(d != 0);
while (d & 1 == 0) d >>= 1; // extract largest odd divisor of d
int exact = n % d == 0;

This is cheaper than changing the FPU rounding mode if you want the rounding mode to be "to nearest" the rest of the time, and there probably exist bit-twiddling tricks that can speed up the extraction of the largest odd divisor of d.

• chux: while (d % 2) d /= 2; has the corner advantage of not relying on 2's complement.
• Pascal Cuoq: @chux Last I checked, compilers didn’t replace division by a simple shift when the integer type was signed. Here the division is exact and the transformation valid but an insufficiently smart compiler may generate something more complex for /=2 than for >>=1
• Eric Postpischil: Given two’s complement and non-zero d, the largest odd divisor of. d is d / (d & -d).
• Geom: @Pascal Thank you, this one looks promising. But the while(d & 1.. is a typo, you mean while(d % 1 ==0), right?
• Pascal Cuoq: @Geom d%1 is always zero, so I definitely do not mean that. I don't think there is a typo, but you can use d%2 if you are not confident that it's the same as d&1 (it is).

Answer 3

Is there a reliable way in C++ to check if double result = dn/dd contains a rounding error?

Should your system allow access to the various FP flags, test for FE_INEXACT after the division.

If FP code is expensive, than at least this code can be used to check integer only solutions.

---

A C solution follow, (I do not have access to a compliant C++ compiler to test right now)

#include <fenv.h>
// Return 0: no rounding error
// Return 1: rounding error
// Return -1: uncertain
#pragma STDC FENV_ACCESS ON
int Rounding_error_detection(int n, int d) {
 double dn = n;
 double dd = d;
 if (feclearexcept(FE_INEXACT)) return -1;
 volatile double result = dn/dd;
 (void) result;
 int set_excepts = fetestexcept(FE_INEXACT);
 return set_excepts != 0;
}

Test code

void Rounding_error_detection_Test(int n, int d) {
 printf("Rounding_error_detection(%d, %d) --> %d\n", 
 n, d, Rounding_error_detection(n,d));
}
int main(void) {
 Rounding_error_detection_Test(3, 6);
 Rounding_error_detection_Test(3, 7);
}

Output

Rounding_error_detection(3, 6) --> 0
Rounding_error_detection(3, 7) --> 1

• Pascal Cuoq: For GCC, it's not enough to have the statement (void) result;, unless I am doing something wrong (godbolt.org/g/DrahMk , GCC famously does not support #pragma STDC FENV_ACCESS ON but -frounding-math is advertised as the file-scope workaround. Using a volatile variable result happens to help.
• chux: @PascalCuoq Code amended per your suggestion.

Answer 4

If the quotient q=dn/dd is exact, it will divide dn exactly dd times.
Since you have dd being integer, you could test exactness with integer division.
Instead of testing the quotient multiplied by dd with (dn/dd)*dd==dn where round off errors can compensate, you should rather test the remainder.
Indeed std:remainder is always exact:

if(std:remainder(dn,dn/dd)!=0)
 std::cout << "Quotient was not exact." << std::endl;