ieee754 floating point 1/x * x > 1.0

Question

I want to know whether the program defined below can return 1 assuming:

• IEEE754 floating point arithmetics
• no overflow (neither in max/x nor in f*x)
• no nan or inf (obviously)
• 0 < x and 0 < n < 32
• no unsafe math optimization

int canfail(int n, double x) {
    double max = 1ULL << n; // 2^n
    double f = max / x;
    return f * x > max;
}

In my opinion, it should sometime return 1, as roundToNearest(max / x) can in general be greater than max/x. I'm able to find numbers for the opposite case, where f * x < max, but I have no examples of input that show f * x > max and I have no idea of how to find one. Can somebody help ?

EDIT: I know the value of x if in a range between 10^(-6) and 10^6 (that still leaves a lot (too much possible double values), but I know I will not have to deal with overflow, underflow or sub-normal numbers ! In addition, I just realized that because max is a power of two and we don't deal with overflow, the solution will be the same by fixing max=1 as it is exactly the same computation, but shifted.

Therefore, the problem correspond to finding a positive, normal double value x such that `(1/x) * x > 1.0 !!

I made a little program to try to find a solution:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdint.h>
#include <omp.h>

int main( void ) {
    #pragma omp parallel
    {
        unsigned short int xsubi[3] = {
            omp_get_thread_num(),
            omp_get_thread_num(),
            omp_get_thread_num()
        };

        #pragma omp for
        for(int64_t i=0; i<INT64_MAX; i++) {
            double x = fmod(nrand48(xsubi), 1048576.0);
            if(x<0.000001)
                continue;

            double f = 1.0 / x;
            if(f * x > 1.0) {
                printf("found !!! x=%.30f\n", x);
                fflush(stdout);
            }
        }
    }
    return 1;
}

If you change the sign of the comparison, you will find some value quickly. However, it seems to run forever with f * x > 1.0

Answer 1

In the absence of underflow or overflow, the exponents are irrelevant; if M/x*x > M, then (M/p) / (x/q) * (x/q) > (M/p) for any powers of two p and q. So let’s consider 2⁵² ≤ x < 2⁵³ and M = 2¹⁰⁵. We can eliminate x = 2⁵² since this yields exact floating-point arithmetic, so 2⁵² < x < 2⁵³.

Division of 2¹⁰⁵ by x yields integer quotient q and integer remainder r, with 2⁵²q < 2⁵³, 0 < r < x, and 2¹⁰⁵ = q•x + r.

In order for M/x*x to exceed M, both the division and the multiplication must round up. Since the division rounds up, x/2 ≤ r.

With rounding up, the result of floating-point division of 2¹⁰⁵ by x yields q+1. Then the exact (not rounded) multiplication yields (q+1)•x = q•x + x = q•x + x + r - r = q•x + r + x − r = 2¹⁰⁵ + x − r. Since x/2 < r, x − r ≤ x/2, so rounding this exact result rounds down, yielding 2¹⁰⁵. (The “<” case always rounds down, and the “=” case rounds down because 2¹⁰⁵ has the low even bit.)

Therefore, for powers of two M and all arithmetic within exponent bounds, M/x*x > M never occurs with round-to-nearest-ties-to-even.

Answer 2

Multiplication by a power of two is just a scaling of exponent, it does not change the problem: so it's the same as finding x such that (1/x) * x > 1.

One solution is brute force search.
For same reasons, we can limit the search of such x in the interval (1.0,2.0(

A better approach is to analyze error bounds without brute force.

Let's note ix the nearest floating point to 1/x.

Considering xand ixas exact fractions, we can write the integer division: 1 = ix * x + r where ris the remainder
(these are all fractions with denominators being powers of 2, so we have to multiply the whole equation by appropriate power of 2 to really have integer division).
In other words, ix = 1/x - r/x, where -r/x is the rounding error of inversion.
When we multiply the inverse approximation by x, the exact value is ix*x = 1 - r.
We know that the floating point result will be rounded to the nearest float to that exact value.
So, assumming default rounding mode to nearest, tie to even, the question asked is whether -r can exceed 0.5 ulp.

The short answer is never!
Suppose |r| > 0.5 ulp, then the rounding error -r/x does exceed half ulp of exact result 1/x.
This is not a proper answer, because the exact result is not a floating point and does not have an ulp, but you get the idea...
I might come back with a correct proof if i have time, but my bet is that you can find it already done, possibly on SO

EDIT Why can you find (1/x) * x < 1?

Simply because 1.0 is at a binade limit, so below 1, we have to prove that r<0.25 ulp, what we cannot...