使用浮动值控制舍入

Question

我在 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的体验几乎为零

Answer 1

正如Eric提到的，fesetround()原则上是控制float舍入模式的方式，但不幸的是它没有得到很好的支持。

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

使用a_bits = msb(a)获取a中设置为1的最高有效位的索引（有关详细信息，请在[=41中搜索“msb” =] 提升数量）。
用uint32_t a_high = a >> (a_bits - 31)得到高32位。

这一步实际上是一个舍入操作，向下舍入。相反，如果 lsb(a) < a_bits - 31，即如果最低设置位出现在提取的 32 位部分下方，我们可以将 a_high 递增 1。（边缘情况是此增量可能会导致 uint32_t 溢出，我将忽略这一点。）我将使用“a_high_down”来指代向下舍入和“a_high_up”进行四舍五入。

无论哪种方式，我们都有股息的浮点数近似值 a:
```
a ≈ 2^(a_bits − 31) ⋅ a_high.
```
同理，做一个约数的近似表示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

使用浮动值控制舍入

Controlling rounding with floating values

boost

floating-point

bigint

例子