ieee754 浮点数 1/x * x > 1.0

Question

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

IEEE754 浮点运算
没有溢出（max/x 和 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;
}

在我看来，它应该在某个时候 return 1，因为 roundToNearest(max / x) 通常可以大于 max/x。我能够找到相反情况的数字，其中 f * x < max，但我没有显示 f * x > max 的输入示例，而且我不知道如何找到一个。有人可以帮忙吗？

编辑：我知道 x 的值如果在 10^(-6) 和 10^6 之间（仍然留下很多（太多可能的双精度值）），但我知道我不必处理溢出、下溢或子- 正常数字！另外，我刚刚意识到因为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 确实超过精确结果 1/x.
的一半 ulp 这不是一个正确的答案，因为确切的结果不是浮点数，也没有 ulp，但你明白了......
如果我有时间，我可能会带回正确的证明，但我敢打赌你会发现它已经完成了，可能在 SO

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

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

Answer 2

在没有下溢或上溢的情况下，指数无关；如果 M/x*x > M，则 (M/p) / (x/q) * (x/q) > (M/p) 对于两个 p 和 q 的任意幂。所以让我们考虑 2⁵² ≤ x < 2⁵³ 和 M = 2¹⁰⁵。我们可以消除 x = 2⁵² 因为这会产生精确的浮点运算，所以 2⁵² < x < 2⁵³.

2¹⁰⁵ 除以 x 得到整数商 q 和整数余数 r，其中 2⁵² q < 2⁵³、0 < r < x 和 2¹⁰⁵ = q•x + r.

为了M/x*x超过M，除法和乘法都必须四舍五入。由于除法向上取整，x/2 ≤ r.

四舍五入后，2¹⁰⁵ 除以 x 的浮点数除法结果为 q+1。然后精确（未四舍五入）乘法产生 (q+1)•x = q•x + x = q•x + x + r - r = q•x + r + x − r = 2¹⁰⁵ + x − r。由于 x/2 < r，x − r ≤ x/2，因此将这个精确结果四舍五入，得到 2¹⁰⁵。（“<”的情况总是向下舍入，“=”的情况向下舍入，因为 2¹⁰⁵ 具有低偶数位。）

因此，对于 2 的幂 M 和指数范围内的所有算术，M/x*x > M 永远不会出现四舍五入到偶数的情况。

Answer 3

canfail(1, pow(2, 1023) * (2 - pow(2, -51))) 将 return 1.

ieee754 浮点数 1/x * x > 1.0

ieee754 floating point 1/x * x > 1.0

c

floating-point

ieee-754