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

问题描述 投票:0回答:1

我的问题是,当你必须除以 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
。然而,在这种方法中,我不清楚数值稳定性会差多少,或者首先有什么需要担心的。

c++ precision division numerical-methods numerical-stability
1个回答
2
投票

(从评论中转移)

忽略次正规数和无穷大,我没有看到推迟除以 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)时,这个过程才会真正终止,之后我们进入无穷大(如果超过)或次正规数,其中尾数开始工作另一种方式。

© www.soinside.com 2019 - 2024. All rights reserved.