使用浮点值控制舍入

编程问答 2022-05-27

问题描述

我正在使用带有 boost 的大整数（128 位或更多）。我想通过使用浮点来避免一些非常缓慢的除法。在许多情况下，我可以只使用商的高位，我想通过使用浮点来使用它。 说我有：

ccp_int a,b; // big integers +ve
double fa,fb; // a float that would be big enough to hold the exponent.
fa/fb >= a/b  so a floor of a divide
(fa - 1) / fb + 1 <= (a - 1) / b + 1  a ceiling of a divide

所以大概这意味着能够分配具有正确舍入的 fa = a say 并且还可以在不同方向上进行除法舍入。 这是你可以用浮点控制的东西吗？我对 fp 的体验几乎为零。

解决方法

就像 Eric 提到的，fesetround() 原则上是控制浮点舍入模式的方式，但不幸的是它没有很好的支持。

这里有一个解决方法的想法。手动提取高 32 位，然后使用 POD 算法计算除法。所以本质上，模拟浮点转换和除法，以便完全控制舍入的方式。为了简化这个大纲，我假设 a 和 b 都是正数，至少为 2^32，希望能清楚如何为一般实现填充这些细节：

  1. 使用 a_bits = msb(a) 获取 a 中设置为 1 的最高有效位的索引（有关详细信息，请在 boost 上的 this page 中搜索“msb”号）。

  2. 使用 uint32_t a_high = a >> (a_bits - 31) 获取高 32 位。

    这一步实际上是一个舍入操作，向下舍入。为了取而代之，我们可以将 a_high 加 1 if lsb(a) < a_bits - 31，也就是说，如果最低设置位出现在提取的 32 位部分下方。 （一个极端的情况是这个增量可能会导致 uint32_t 溢出，我将忽略它。）我将使用“a_high_down”来指代四舍五入和“a_high_up”来指代四舍五入。

    无论哪种方式，我们都有股息a的浮动近似值：

    a ≈ 2^(a_bits − 31) ⋅ a_high.

  3. 同理，对除数b做一个近似表示：

    b ≈ 2^(b_bits − 31) ⋅ b_high.

那么，除法 a/b 大约是

a/b ≈ 2^(a_bits − b_bits − 32) ⋅ (2^32 ⋅ a_high) // b_high,

其中 // 表示整数除法向下舍入。分子 2^32 ⋅ a_high 应通过将 a_high 转换为 uint64_t 并向上移动来形成。

一些注意事项：

  • 这种调整很重要，以便商保持一些精度（因为 a_high 和 b_high 具有相似的幅度，因此简单地将 a_high // b_high 作为整数除法只会给出 1 位精度）。
  • 我们可以将 a_high 提取为 a 的高64 位，而不是这种 32 位移位，以获得更高的精度。

要从上到下限制商，计算

a/b ≤ 2^(a_bits − b_bits − 32) ⋅ ((2^32 ⋅ a_high_up − 1) // b_high_down + 1),a/b ≥ 2^(a_bits − b_bits − 32) ⋅ (2^32 ⋅ a_high_down) // b_high_up,

其中 // 表示整数除法向下舍入。 rhs 上的除法可以在本机 uint64_t 算术中完成。最后，将这些除法结果转换为 cpp_int 并左移 (a_bits − b_bits − 32)。

示例

Ex 1.对于值

a = 2^100 ⋅ 123456 = 156499072501776288991176990922899456,b = 2^70 ⋅ 123455 = 145749938535668012464209920,

确切的商大约是 a/b = 1073750521.434887。上面的过程准确地近似了这一点，前 10 位数字是正确的：

a_bits = floor(log2(a)) = 116
a_high = a >> (116 − 31) = 4045406208

b_bits = floor(log2(b)) = 86
b_high = b >> (86 − 31) = 4045373440

a/b ≈ 2^(a_bits − b_bits − 32) ⋅ (2^32 ⋅ a_high) // b_high
    = 2^(−2) ⋅ (2^32 ⋅ 4045406208) // 4045373440
    = 0.25 ⋅ 4295002085
    = 1073750521.25.

Ex 2.对于值

a = 10^50 = 100000000000000000000000000000000000000000000000000,b = 9^50 = 515377520732011331036461129765621272702107522001,

精确商约为 a/b = 194.03252175，并且

a_bits = floor(log2(a)) = 166
a_high_down = a >> (166 − 31) = 2295887403
a_high_up = a_high_down + 1 = 2295887404

b_bits = floor(log2(b)) = 158
b_high_down = b >> (158 − 31) = 3029116820
b_high_up = b_high_down + 1 = 3029116821

a/b ≤ 2^(a_bits − b_bits − 32) ⋅ ((2^32 ⋅ a_high_up − 1) // b_high_down + 1)
    = 2^(−24) ⋅ (9860761315478339583 // 3029116820 + 1)
    = 2^(−24) ⋅ 3255325530
    = 194.03252184

a/b ≥ 2^(a_bits − b_bits − 32) ⋅ (2^32 ⋅ a_high_down) // b_high_up
    = 2^(−24) ⋅ (9860761311183372288 // 3029116821)
    = 2^(−24) ⋅ 3255325526
    = 194.03252161
bigintboostboostfloating-point