ldexp 应该正确舍入

Posted 2023-02-22

【中文标题】ldexp 应该正确舍入

【英文标题】：should ldexp round correctly

【发布时间】：2015-08-21 23:56:41

【问题描述】：

我对 MSVC ldexp 行为感到有点惊讶（它发生在 Visual Studio 2013 中，但也发生在至少 2003 年之前的所有旧版本中......）。

例如：

#include <math.h>
#include <stdio.h>
int main()

    double g=ldexp(2.75,-1074);
    double e=ldexp(3.0,-1074);

    printf("g=%g e=%g \n",g,e);
    return 0;

打印

g=9.88131e-324 e=1.4822e-323

第一个 g 奇怪地圆了... 它是 2.75 * fmin_denormalized，所以我绝对期待第二个结果 e. 如果我评估 2.75*ldexp(1.0,-1074) 我正确地得到与 e 相同的值。

是我的期望太高，还是微软没有遵守某些标准？

【参考方案1】：

虽然问题没有明确说明这一点，但我假设提问者期望的输出是：

g=1.4822e-323 e=1.4822e-323

这是我们对承诺严格遵守 IEEE-754 的 C/C++ 编译器的期望。这个问题同时被标记为C 和C++，我将在这里解决C99，因为这是我手头的标准。

在描述 IEC 60559 浮点运算（其中 IEC 60559 基本上是 IEEE-754 的另一个名称）的附录 F 中，C99 标准规定：

定义__STDC_IEC_559__ 的实现应符合 本附件中的规范。 [...] scalbn 和 scalbln <math.h> 中的函数提供了推荐的 scalb 函数 IEC 60559 附录。

在该附件的下方，F.9.3.6 节规定：

在二进制系统上，ldexp(x, exp) 等价于 scalbn(x, exp)。

C99 标准引用的附录是 1985 版 IEEE-754 的附录，我们发现 scalb 函数定义如下：

Scalab(y, N) 为整数值 N 返回 y × 2^N 而不计算 2^N。

scalb 定义为2的幂次的乘法，乘法必须根据标准根据当前的舍入模式正确舍入。因此，使用符合标准的C99 编译器ldexp() 必须返回正确舍入的结果如果 编译器定义__STDC_IEC_559__。在没有设置舍入模式的库调用时，默认的舍入模式“舍入到最接近或偶数”有效。

我无权访问 MSVC 2013，所以我不知道它是否定义了该符号。这甚至可能取决于编译器标志设置，例如/fp:strict。

在找到我的 C++11 标准副本后，我找不到任何对 __STDC_IEC_559__ 或任何关于 IEEE-754 绑定的语言的引用。根据this question 的回答，这是因为引用 C99 标准包含了该支持。

【讨论】：

这种行为也发生在 /fp:strict 中。 STDC_IEC_559 没有定义，所以它至少不承诺做对......什么标准 MSVC 2013 在 C/C++ 中遵守是另一个问题......

【参考方案2】：

发生这种情况是因为在ldexp 计算过程中，2.75 被截断为 2，这是因为在非规范化数字中，代表“.75”部分的位会从可表示数字的末尾移出并消失.这是一个错误还是设计的行为可以争论。

在计算2.75*ldexp(1.0,-1074) 时发生正常舍入，2.75 变为 3。

编辑：ldexp 应该正确舍入，这是一个错误。

【讨论】：

我认为 OP 是专门要求进行这场辩论的。

@Potatoswatter 我在调试器中花了很多时间来弄清楚 为什么 会发生这种情况，但我忽略了 OP 的要求。我会投票给 bug。

【参考方案3】：

OP 结果不符合 C 规范，因为该规范没有定义计算的精确性。

OP 结果可能未通过 IEEE 754，但这取决于当时使用的舍入模式，未发布。然而 OP 的报告 2.75*ldexp(1.0,-1074) 工作正常，这意味着当时预期的舍入模式已经生效。

使用printf("%la",x) 有助于清楚地看到double 的极限附近发生了什么。

我希望 g 会“四舍五入到最接近 - 平局”，结果是 0x1.8p-1073 - 这确实发生在 Windows 上的 gcc 中。

理想情况下，g 的值应为 0x1.6p-1073

0x0.0p-1073 Zero
0x0.8p-1073 next higher double DBL_TRUE_MIN
0x1.0p-1073 next higher double
0x1.6p-1073 ideal `g` answer, but not available as a double
0x1.8p-1073 next higher double

公平地说，这可能是一个处理器错误 - 它有 happened before。

---

参考

double g=ldexp(2.75,-1074);
printf("%la\n%la\n", 2.75,ldexp(2.75,-1074));
printf("%la\n%la\n", 3.0 ,ldexp(3.0 ,-1074));
double e=ldexp(3.0,-1074);
printf("%la\n%la\n", g,e);
printf("%la\n%la\n", 9.88131e-324, DBL_TRUE_MIN);
printf("g=%g e=%g \n",g,e);

0x1.6p+1
0x1.8p-1073
0x1.8p+1
0x1.8p-1073
0x1.8p-1073
0x1.8p-1073
0x1p-1073
0x1p-1074
g=1.4822e-323 e=1.4822e-323 

【讨论】：

通过调试器中的代码，这是一个算法错误，而不是处理器错误。