Exact double division

Question

Consider the following function:

auto f(double a, double b) -> int
{
  return std::floor(a/b);
}

So I want to compute the largest integer k such that k * b <= a in a mathematical sense. As there could be rounding errors, I am unsure whether the above function really computes this k. I do not worry about the case that k could be out of range. What is the proper way to determine this k for sure?

Answer 1

It depends how strict you are. Take a double b and an integer n, and calculate bn. Then a will be rounded. If a is rounded down, then it is less than the mathematical value of nb, and a/b is mathematically less than n. You will get a result if n instead of n-1.

On the other hand, a == b*n will be true. So the “correct” result could be surprising.

Your condition was that “kb <= a”. If we interpret this as “the result of multiplying kb using double precision is <= a”, then you’re fine. If we interpret it as “the mathematically exact product of k and b is <= a”, then you need to calculate k*b - a using the fma function and check the result. This will tell you the truth, but might return a result of 4 if a was calculated as 5.0 * b and was rounded down.

Answer 2

The problem is that float division is not exact.

a/b can give 1.9999 instead of 2, and std::floor can then give 1.

One simple solution is to add a small value prior calling std::floor:

std::floor (a/b + 1.0e-10);

Result:

result = 10 while 11 was expected
With eps added, result = 11

Test code:

#include <iostream>
#include <cmath>

int main () {
    double b = atan (1.0);
    int x = 11;
    double a = x * b;
    int y = std::floor (a/b);
    std::cout << "result = " << y << " while " << x << " was expected\n";
    
    double eps = 1.0e-10;
    int z = std::floor (a/b + eps);
    std::cout << "With eps added, result = " << z << "\n";
    
    return 0;
}