仅使用就近偶舍入模式计算IEEE-754双精度trunc(a/b)的方案咨询

原始需求与标准实现 #

我们需要对两个IEEE-754双精度浮点数a和b的精确商a/b做向零取整，得到整数结果。原本可以通过C99标准的fesetround函数切换舍入模式为FE_TOWARDZERO实现，示例代码如下：

#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;  // 舍入操作在此处执行
  fesetround(old_mode);
  return trunc(result);
}

#include <stdio.h>
int main() {
  // 应当输出"6004799503160662"，因为18014398509481988 / 3 = 6004799503160662.666...
  printf("%.17g", trunc_div(18014398509481988.0, 3.0));
}

场景限制 #

但部分场景下仅能使用nearest-ties-to-even（就近偶舍入）模式，比如开启优化的GCC编译环境、微控制器编译场景、JavaScript运行环境都存在这类限制。

现有实现方案 #

误差补偿方案 #

目前已尝试的方案为：按就近偶舍入模式计算a/b后取整，再通过Veltkamp-Dekker拆分实现的mul_error函数计算乘法误差，对结果幅值偏大的场景做补偿，实现代码如下：

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用于计算精确乘法误差，实现代码如下：

// 返回a的最高26个有效位
// 假设fabs(a) < 1e300，避免乘法溢出
double highbits(double a) {
  double p = 0x8000001L * a;
  double q = a - p;
  return p + q;
}

// 计算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;  // 后续乘法均为精确运算
  e += ah*bl;
  e += al*bh;
  e += al*bl;
  return ldexp(e, a_exp + b_exp);
}

优化后的mul_error实现 #

现有mul_error函数已可通过边界判断减少frexp、ldexp的调用，在随机输入下可获得30%的性能提升，对应优化代码框架如下：

double mul_error(double a, double b) {
  if (!isfinite(a*b)) return -a*b;
  double A = fabs(a), B = fabs(b);
  // 边界值来自Dekker乘法的验证结论
  if (A>0x1p995 || B>0x1p995 || (A*B!=0 && (A*B<0x1p-969 || A*B>0x1p1021))) {
    // 可能出现溢出/下溢：使用frexp、ldexp处理
  } else {
    // 无需调用frexp、ldexp
  }
}

128位整数参考实现 #

另有基于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;  // 整数除法自动向零取整
  if (r & 0xFFE0000000000000) { r >>= 1; scale++; }  // 归一化处理
  
  // scale<0表示存在小数位，直接移出小数位即可
  double d = scale<-63 ? 0 : scale<0 ? r>>-scale : ldexp(r, scale);
  
  // 溢出时返回双精度最大有限值
  return sign * (isfinite(d) ? d : 0x1.fffffffffffffp1023); 
}

待咨询问题 #

• 上述补偿方案是否会在部分输入场景下失效，例如溢出、下溢场景？
• 是否存在运算效率更高的实现方式？

补充说明 #

1. 已修正mul_error函数边界处理逻辑，修复了a为±∞时的异常问题；若a、b为有限非零值且a/b运算溢出，需匹配IEEE-754向零舍入模式的除法行为，返回双精度最大有限值±(2¹⁰²⁴ − 2⁹⁷¹)。

内容的提问来源于stack exchange，提问作者Řrřola