ieee754 浮点数 1/x * x > 1.0

Question

我想知道下面定义的程序是否可以返回 1 假设：

• IEEE754 浮点运算
• 没有溢出（既不是 inmax/x也不是 in f*x）
• 没有 nan 或 inf（显然）
• 0 < x 和 0 < n < 32
• 没有不安全的数学优化

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。我能够找到相反情况下的数字 where f * x < max，但我没有显示的输入示例，f * x > max我不知道如何找到一个。有人可以帮忙吗？

编辑：如果在 10^(-6) 和 10^6 之间的范围内，我知道 x 的值（仍然有很多（可能的双值），但我知道我不必处理溢出、下溢或次正规数！另外，我刚刚意识到，因为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;
}

如果你改变比较的符号，你会很快找到一些值。然而，它似乎永远运行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超过精确结果的一半 ulp 1/x。
这不是一个正确的答案，因为确切的结果不是浮点并且没有 ulp，但你明白了......
如果我有时间，我可能会带着正确的证明回来，但我打赌你可以发现它已经完成了，可能在SO上

编辑 为什么你能找到(1/x) * x < 1？

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