舍入误差检测

Question

我有两个整数n和d。这些可以用双 dn(n) 和双 dd(d) 精确表示。 C++ 中是否有可靠的方法来检查

double result = dn/dd

包含舍入错误？如果它只是一个整数除法检查 (n/d) * d==n 是否可行，但使用双精度算法执行此操作可能会隐藏舍入错误。

编辑：发布后不久，令我震惊的是，将舍入模式更改为 round_down 将使 (n/d)*d==n 测试适用于双倍。但如果有更简单的解决方案，我还是想听听。

Answer 1

如果您忽略上溢和下溢（除非代表 d 和 n 的整数类型非常宽，否则您应该能够做到这一点），那么（二进制）浮点除法 dn/dd 是精确的当且仅当 d 是 n 乘以 2 的幂的除数。

检查此问题的算法可能如下所示：

assert(d != 0);
while (d & 1 == 0) d >>= 1; // extract largest odd divisor of d
int exact = n % d == 0;

如果您希望舍入模式在其余时间“最接近”，这比更改 FPU 舍入模式更便宜，并且可能存在可以加快提取最大奇数的位旋转技巧d.

的除数

Answer 2

如果硬件 FMA 可用，那么在大多数情况下（预计 n 不小的情况，如下所示），最快的测试可能是：

#include <cmath>
…
double q = dn/dd;
if (std::fma(-q, dd, dn))
    std::cout << "Quotient was not exact.\n";

如果 nd−q•dd 太小以至于四舍五入，这可能会失败为零，如果其幅度小于最小可表示正值的一半（通常为 2⁻¹⁰⁷⁴），则在舍入到最接近偶数模式中发生。这只有在 dn 本身很小的情况下才会发生。如果需要，我希望我可以计算 dn 的一些界限，并且鉴于 dn = n 和 n 是整数，这不应该发生。

忽略指数边界，测试可整除有效数字的方法是：

#include <cfloat>
#include <cmath>
…
int sink; // Needed for frexp argument but will be ignored.
double fn = std::ldexp(std::frexp(n, &sink), DBL_MANT_DIG);
double fd = std::frexp(d, &sink);
if (std::fmod(fn, fd))
    std::cout << "Quotient will not be exact.\n";

鉴于 n 和 d 是在浮点类型中可以精确表示的整数，我想我们可以显示它们的指数不能使上述测试失败。有些情况下 n 是一个小整数而 d 是一个大整数（从 2¹⁰²³ 到2¹⁰²⁴−2⁹⁷²,包括）我需要考虑。

Answer 3

Is there a reliable way in C++ to check if double result = dn/dd contains a rounding error?

如果您的系统允许访问各种 FP 标志，请在除法后测试 FE_INEXACT。

如果 FP 代码很昂贵，那么至少此代码可用于检查仅整数解决方案。

接下来是 C 解决方案，（我现在无法访问兼容的 C++ 编译器来进行测试）

#include <fenv.h>
// Return 0: no rounding error
// Return 1: rounding error
// Return -1: uncertain
#pragma STDC FENV_ACCESS ON
int Rounding_error_detection(int n, int d) {
  double dn = n;
  double dd = d;

  if (feclearexcept(FE_INEXACT)) return -1;
  volatile double result = dn/dd;
  (void) result;
  int set_excepts = fetestexcept(FE_INEXACT);
  return set_excepts != 0;
}

测试代码

void Rounding_error_detection_Test(int n, int d) {
  printf("Rounding_error_detection(%d, %d) --> %d\n", 
    n, d, Rounding_error_detection(n,d));
}

int main(void) {
  Rounding_error_detection_Test(3, 6);
  Rounding_error_detection_Test(3, 7);
}

输出

Rounding_error_detection(3, 6) --> 0
Rounding_error_detection(3, 7) --> 1

Answer 4

如果商 q=dn/dd 是精确的，它会整除 dn 正好是 dd 次。
因为你有 dd 是整数，你可以用整数除法来测试准确性。
与其用 (dn/dd)*dd==dn 测试商乘以 dd 的结果，舍入误差可以补偿，您应该测试余数。
事实上 std:remainder 总是准确的：

if(std:remainder(dn,dn/dd)!=0)
    std::cout << "Quotient was not exact." << std::endl;

舍入误差检测

Rounding error detection

c++

floating-point

rounding-error