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.

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.

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

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;