Robust unit test for reciprocal approximation

Question

Let $x$ and $y$ be representable floating point numbers. I'm looking for a unit test which can ensure that my user's compiler has not made the reciprocal approximation $\mathrm{fl}(x/y) \approx \mathrm{fl}(x)*\mathrm{fl}(1/y))$ as an optimization. However, I have been unsuccessful in coming up with a compelling example.

It's fairly easy to come up with examples where $\mathrm{fl}(x/y) \ne \mathrm{fl}(x)*\mathrm{fl}(1/y)$, e.g.

\begin{align*} \mathrm{fl}_{64}(1.3)/\mathrm{fl}_{64}(2.9) &= \mathrm{fl}_{64}(0.4482758620689655) \\ \mathrm{fl}_{64}(1.3)*(\mathrm{fl}_{64}(1)/\mathrm{fl}_{64}(2.9)) &= \mathrm{fl}_{64}(0.4482758620689656) \end{align*}

However, for two single numbers, the compiler might evaluate them at compile time and not make the reciprocal approximation, which defeats the point of my unit test to detect this. And even if the numbers are only known at runtime, single divisions are generally not optimized-only array divisions.

How can I design a sequence of double precision representables $x_i$ and a representable $y$ such that $\mathrm{fl}(x_i/y) \ne \mathrm{fl}(x_i)*\mathrm{fl}(1/y)$?

FWIW I'm using C++, though ostensibly that's of little relevance to the question.

Answer 1

Compilers that support turning division into a multiplication with the reciprocal with some optimization flag settings (e.g. -ffast-math or -fp-model=fast) may or may not apply this to any particular division operator. It is therefore impossible to design a unit test that is 100% reliable. It is always possible that the compiler does not apply the transformation to the test case from the unit test, but does apply it to code one cares about.

The following assumes that a system supports the IEEE-754 (2008) floating-point standard, that the double type of C++ is bound to IEEE-754 binary64 type, and that the rounding mode of "to nearest or even" is in effect. These conditions are met by most modern systems, in particular C++ environments running on x86-64 and ARM64. For the proposed technique to work, the standard math function fma() (C++11 or higher) must be available, and must be correctly implemented, preferably by mapping to a corresponding hardware instruction, as spme fma() implementations based on software emulation have been found to be faulty.

Correct rounding of a division q = x / y can be checked by computing the residual x - q * y accurately with the help of a fused multiply-add (FMA): residual = fma (q, -y, x). If we choose to divide by 3.0, residuals of dividends that are consecutive binary64 encodings form a repeating pattern that repeats every three dividends. One third of cases result in a positive residual, one third result in a negative residual, and one third has a zero residual. The code below simply checks these statistics instead of examining the pattern in detail. It also checks whether the residual is bigger than would be expected from division with round-to-nearest-or-even. If division is replaced by multiplication with the reciprocal, both of these checks (residual statistics, maximum residual magnitude) will fail.

For some divisors, it is possible to correctly compute the quotient of a floating-point division by compile-time literal constant without using a division operation, and this includes the case x / 3.0. I am not aware of a compiler that currently does this. To prevent a compiler from discovering the compile-time constant divisor of 3.0, it would therefore be highly advisable to prevent inlining of is_division_correctly_rounded(). Compilers typically offer toolchain-specific function attributes for this, e.g. __attribute__((noinline)), __declspec(noinline).

Note that hexadecimal floating-point literals were added in C++17, so the code below should be compiled accordingly.

#include <cstdio>
#include <cstdlib>
#include <cmath>

#define NBR_PER_CASE   (100)
#define NBR_TEST_CASES (3 * NBR_PER_CASE)

int is_division_correctly_rounded (double dividend, double divisor, double eps)
{
    int count_less = 0;
    int count_more = 0;
    int count_equal = 0;
    int count_too_big = 0;
    for (int i = 0; i < NBR_TEST_CASES; i++) {
        double quotient = dividend / divisor;
        double residual = fma (quotient, -divisor, dividend);
        count_less += (residual == -eps);
        count_more += (residual == +eps);
        count_equal += (residual == 0.0);
        count_too_big += (fabs (residual) > eps);
        dividend = nextafter (dividend, INFINITY);
    };
    return ((count_less == NBR_PER_CASE) && (count_more == NBR_PER_CASE) &&
            (count_equal == NBR_PER_CASE) && (count_too_big == 0)) ? 1 : 0;
}

int main (void)
{
    if (is_division_correctly_rounded (2.5, 3.0, 0x1.0p-53)) {
        printf ("division appears to be correctly rounded to nearest-even\n");
    } else {
        printf ("division NOT correctly rounded to nearest-even\n");
    }
    return EXIT_SUCCESS;
}
```