c - ieee754 float 1/x * x> 1.0

Question

我想知道以下定义的程序是否可以返回1:

IEEE754浮点算法

没有溢出(max/x和f*x都没有)

no nan或inf(显然)

0

没有不安全的数学优化

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

我认为它应该有时返回1，因为 roundToNearest(max / x)通常可以大于 max/x。
我能够找到相反情况下的数字 f * x < max，但是我没有显示 f * x > max的输入示例，也不知道如何找到一个。有人可以帮忙吗？

编辑:
我知道x的值是否在10 ^(-6)和10 ^ 6之间(仍然留下很多(太多的 double 值))，但是我知道我不必处理上溢，下溢或子运算-正常数字!
另外，我刚刚意识到，因为 max是2的幂，并且我们不处理溢出，所以 的解决方案将是相同的，因为它固定了max=1 ，因为它是完全相同的计算，但是发生了变化。

因此，问题对应于找到一个正的，正常的双值 x，使得((1/x)* x> 1.0 !!!

我编写了一个小程序来尝试找到解决方案:

#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;
}

如果更改比较的符号，您将很快找到一些值(value)。但是，它似乎与 f * x > 1.0一起永远运行

Answer 1

乘以2的幂只是指数的缩放，它不会改变问题:因此，它与找到x这样(1/x) * x > 1相同。

一种解决方案是蛮力搜索。
出于相同的原因，我们可以将搜索x的时间间隔限制为(1.0,2.0(

更好的方法是在没有蛮力的情况下分析误差范围。

让我们注意ix最接近1/x的浮点。

将x和ix视为精确分数，我们可以写整数除法:1 = ix * x + r其中r是余数
(这些都是分母为2的幂的分数，因此我们必须将整个方程乘以2的适当幂才能真正具有整数除法)。
换句话说，ix = 1/x - r/x，其中-r/x是求反的舍入误差。
当我们将反近似值乘以x时，精确值为ix*x = 1 - r。
我们知道浮点结果将四舍五入到最接近该精确值的浮点数。
因此，假设默认的舍入模式最接近，甚至为偶数，问的问题是-r是否可以超过0.5 ulp。

简短的答案永远不会!
假设|r| > 0.5 ulp，那么舍入误差-r/x确实超过了精确结果1/x的一半ulp。
这不是一个正确的答案，因为确切的结果不是浮点数，没有浮点数，但是您知道了...
如果我有时间，我可能会带回正确的证明，但我敢打赌，您可能已经发现它已经完成，可能是这样

编辑
为什么可以找到(1/x) * x < 1？

仅仅因为1.0处于binade极限，所以低于1，我们必须证明r<0.25 ulp，我们不能...