2

给定两个 IEEE-754 双精度浮点数a和b，我想得到精确的商a / b四舍五入到零的整数。
执行此操作的 C99 程序可能如下所示：

#include <fenv.h>
#include <math.h>
#pragma STDC FENV_ACCESS on

double trunc_div(double a, double b) {
  int old_mode = fegetround();
  fesetround(FE_TOWARDZERO);
  double result = a/b;  // rounding occurs here
  fesetround(old_mode);
  return trunc(result);
}

#include <stdio.h>
int main() {
  // should print "6004799503160662" because 18014398509481988 / 3 = 6004799503160662.666...
  printf("%.17g", trunc_div(18014398509481988.0, 3.0));
}

现在假设我只能使用最接近偶数舍入模式：我可能正在使用带有优化的 GCC，为微控制器编译，或者必须使其在 JavaScript 中工作。

我尝试的是使用提供的舍入计算a / b，如果结果的幅度太大，则截断并进行补偿：

double trunc_div(double a, double b) {
  double result = trunc(a/b);
  double prod = result * b;
  
  if (a > 0) {
    if (prod > a || (prod == a && mul_error(result, b) > 0)) {
      result = trunc(nextafter(result, 0.0));
    }
  }
  else {
    if (prod < a || (prod == a && mul_error(result, b) < 0)) {
      result = trunc(nextafter(result, 0.0));
    }
  }

  return result;
}

辅助函数mul_error计算精确的乘法误差（使用 Veltkamp-Dekker 分裂）：

// Return the 26 most significant bits of a.
// Assume fabs(a) < 1e300 so that the multiplication doesn't overflow.
double highbits(double a) {
  double p = 0x8000001L * a;
  double q = a - p;
  return p + q;
}

// Compute the exact error of a * b.
double mul_error(double a, double b) {
  if (!isfinite(a*b)) return -a*b;
  int a_exp, b_exp;
  a = frexp(a, &a_exp);
  b = frexp(b, &b_exp);
  double ah = highbits(a), al = a - ah;
  double bh = highbits(b), bl = b - bh;
  double p = a*b;
  double e = ah*bh - p;  // The following multiplications are exact.
  e += ah*bl;
  e += al*bh;
  e += al*bl;
  return ldexp(e, a_exp + b_exp);
}

某些输入的补偿是否会失败（例如，由于上溢或下溢）？
有更快的方法吗？

编辑：将第一行的mul_errorfrom更改… return a*b为… return -a*b;. 这修复了a = ±∞ 的情况；有限输入没问题。
感谢Eric Postpischil发现错误。

编辑：如果a，b是有限且非零且除法a / b溢出，我想在舍入到零模式下匹配 IEEE-754 除法，它返回最大有限双精度数 ±(2¹⁰²⁴ − 2⁹⁷¹)。

编辑：函数frexp和ldexp只能在必要时调用。这对于具有均匀随机位
的双倍a，b来说是 30% 的加速。

double mul_error(double a, double b) {
  if (!isfinite(a*b)) return -a*b;
  double A = fabs(a), B = fabs(b);
  // bounds from http://proval.lri.fr/gallery/Dekker.en.html
  if (A>0x1p995 || B>0x1p995 || (A*B!=0 && (A*B<0x1p-969 || A*B>0x1p1021))) {
    // ... can overflow/underflow: use frexp, ldexp
  } else {
    // ... no need for frexp, ldexp
  }
}

也许ldexp总是不必要的，因为我们只需要知道 mul_error 与 0 的比较。

编辑：如果您有 128 位整数可用，请按照以下步骤操作。（它比原始版本慢。）

double trunc_div(double a, double b) {
  typedef uint64_t u64;
  typedef unsigned __int128 u128;

  if (!isfinite(a) || !isfinite(b) || a==0 || b==0) return a/b;

  int sign = signbit(a)==signbit(b) ? +1 : -1;
  int ea; u64 ua = frexp(fabs(a), &ea) * 0x20000000000000;
  int eb; u64 ub = frexp(fabs(b), &eb) * 0x20000000000000;
  int scale = ea-53 - eb;
  u64 r = ((u128)ua << 53) / ub;  // integer division truncates
  if (r & 0xFFE0000000000000) { r >>= 1; scale++; }  // normalize
  
  // Scale<0 means that we have fractional bits. Shift them out.
  double d = scale<-63 ? 0 : scale<0 ? r>>-scale : ldexp(r, scale);
  
  // Return the maximum finite double on overflow.
  return sign * (isfinite(d) ? d : 0x1.fffffffffffffp1023); 
}
4

1 回答 1

0

考虑确切的余数r=frem(a,b)。

我们知道，a = b*n + r对于某个整数 n，r 介于 -b/2 和 b/2 之间。

并且介于 -1/2 和 1/2 之间（/ 是这里的精确除法）a/b = n + r/b。r/b

我们可以想象两种情况，float(a/b)四舍五入到上整数部分：

  • 当余数为负时（n 的反号），并且如此之小以至于float(n+r/b)=n
  • 当n自身太大而无法表示为浮点数时

第一种情况的一个例子是

a=ldexp(1.0,53); // 2^53, the successor of 2^53-1
b=nextafter(6361.0,7000.0); // close to exact division because 2^53-1=6361*69431*20394401
r=frem(a,b); // -0.287...

在这种情况下，n=1416003655831四舍五入float(a/b)到 n，余数-r/b小于ulp(n).

请注意，测试a > 0 && fma(result,b,-a) > 0是可以的，但nextafter(result,0.0)在这种情况下不能调整，它会导致非整数结果1416003655830.999755859375。我们宁愿拿result-1when trunc(a/b) < 2^53。

以第二种情况为例：

a=ldexp(1.0,54); // 2^54
b=nextafter(1.0,0.0);
r=frem(a,b); // 2.22...e-16

我们有 n 是 2^54+2，a 和 nextafter(2,2*a) 之间的确切中点，
余数为正数 r，trunc(float(a/b))将四舍五入到 a+4。
并且关于第一种情况下显示的符号的讨论r在这里不起作用，所以不能一概而论......

请注意，通过适当的缩放，第二种情况总是可以减少到第一种情况：

int exp,scale;
double result=a/b;
frexp(result,&exp);
scale=53-exp;
if(scale<0)
    return ldexp( trunc_div(ldexp(a,scale),b) , -scale );

但这并没有实际意义，第一种情况仍然需要调整结果以进行汇总。

因此，正如我们在第一个示例中看到的那样，调整可能无法回答整数，而且这个答案没有显示出更快的方法，可能没有太多收获。

于 2021-12-09T00:27:14.440 回答