在 C++ 中处理 double 除以 2 的大幂

Question

我的问题是，当你必须除以 2 的大幂时，应该如何处理除法，这可能是一件微不足道的事情，但我没有找到任何有用的材料。我基本上是在问两种提议的方法中哪一种（如果有的话）最终更明智和/或是否有更好的方法来实现我的目标。

背景是，在我正在从事的一个项目中，我必须在某一时刻计算

double

类型数字

w

形式的总和，其形式为

w = (c_1/2)*(c_2/2)*...*(c_n/2)

与

c_1,...,c_n

一些其他类型

double

数字。在优化代码的过程中，我首先想到先计算

c_1,...,c_n

的乘积，然后将该乘积除以

2^n

可能是个好主意。由于我无法在这篇文章中真正了解的原因，定义

n

的产品中的

w

可能会变得相当大，甚至可能在 70 或 80 左右。现在我正在有效地计算

w

i 的最终值。 ) 计算上述

c_1,...,c_n

的乘积并将其存储到

w

，ii.) 计算

2^n

的值并将其存储到

long long

类型变量

a

，iii.) 设置为

w = w/a

。

据我所知，

long long

类型变量中可以存储的最大2的幂是63，所以我想知道在for循环中除以

w

是否会更好，比如

for (int k { 0 }; k < n; ++k) {
    w /= 2;
}

避免在上述变量

a

中存储过大的值。我知道可以通过重复除以适当的 2 次幂来将其减少到对数时间，但这一点仍然存在。

或者，我也可以不从乘积中分解出 1/2，而是以旧方式计算

w

。然而，在这种方法中，我不清楚数值稳定性会差多少，或者首先有什么需要担心的。

Answer 1

（从评论中转移）

忽略次正规数和无穷大，我没有看到推迟除以 2 的幂的任何数值稳定性优势：因为

double

值在内部存储为整数尾数乘以整数基数 2 指数，所有除以 2 的沸腾调整双精度数的指数部分；尾数保持不变，并且没有精度损失。

考虑到除以 2 后，您将所有这些数字相乘，这不会影响最终结果的精度，因为在双倍乘法中，指数部分只是相加（再次准确）。

因此，将 c_i/2 个数字相乘，或者将 c_i 相乘并除以

exp2(n)

（这同样会返回

double

的整数值的精确

n

值，因此您无需累加 2 的幂）

long long

变量中的值）最终会产生完全相同的值，因为唯一的区别在于您是在之前还是之后调整指数，这在两种情况下都是“精确”操作。

您可以通过一些模糊测试自行测试：

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

struct XorShift64Star {
    /// PRNG state
    uint64_t y;

    /// Initializes with seed
    explicit XorShift64Star(uint64_t seed = 0) : y(seed) {
        if(y == 0) y = 0x159a55e5075bcd15ULL;
    }

    /// Returns a value in the range [1, 1<<64)
    uint64_t operator()() {
        y ^= (y>>12);
        y ^= (y<<25);
        y ^= (y>>27);
        return y * 0x2545F4914F6CDD1DULL;
    }

    /// 53-bit resolution double in range [0, 2^64)
    double gen_big_double() {
        uint64_t x = (*this)();
        double gexp = x & 0x3f;
        if (x & 0x40) gexp = -gexp;
        return (x >> 11)/9007199254740992.0 * exp2(gexp);
    }
};

int main() {
    XorShift64Star gen;
    double res_half = 0.;
    double res_full = 0.;
    long long iters;
    for (iters = 0; res_half == res_full / exp2(iters); ++iters) {
        double val = gen.gen_big_double();
        res_half *= val / 2.;
        res_full *= val;
    }
    printf("%lld %0.18g %0.18g", iters, res_half, res_full / exp2(iters));
    return 0;
}

我让它运行了分钟一个多小时，但它仍然没有退出；我希望只有当其中一个中间结果达到指数的极限（±2**127）时，这个过程才会真正终止，之后我们进入无穷大（如果超过）或次正规数，其中尾数开始工作另一种方式。