I want to know whether the program defined below can return 1 assuming:
max/x
nor in f*x
)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
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 x
and ix
as exact fractions, we can write the integer division: 1 = ix * x + r
where r
is 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...