float128和double-double算法

Question

我已经在Wikipedia中看到，某种程度上实现四精度的方法是使用双精度双精度算术，即使就位而言它的精度不完全相同：https : //en.wikipedia.org/wiki/Quadruple-precision_floating-point_format

在这种情况下，我们使用两个double来存储值。因此，我们进行了两次运算来计算结果，结果的每一倍进行一次运算。

在这种情况下，我们可以在每个双精度数上产生舍入误差，或者是避免这种误差的一种机制？

Answer 1

“在这种情况下，我们使用两个double来存储值。因此，我们需要每次进行两次操作。”

这不是double-double算法的工作原理。您应该期望一个双精度双精度运算可以在6到20个双精度运算中的任意一个位置实现，具体取决于实现的实际运算，融合乘加运算的可用性，一个操作数大于另一个运算数的假设……

例如，这是一种双倍乘法的实现，用于当FMA指令不可用时（取自CRlibm）：

#define Mul22(zh,zl,xh,xl,yh,yl)                      \
{                                                     \
double mh, ml;                                        \
                              \
  const double c = 134217729.;                \
  double up, u1, u2, vp, v1, v2;              \
                              \
  up = (xh)*c;        vp = (yh)*c;            \
  u1 = ((xh)-up)+up;  v1 = ((yh)-vp)+vp;          \
  u2 = (xh)-u1;       v2 = (yh)-v1;                   \
                              \
  mh = (xh)*(yh);                     \
  ml = (((u1*v1-mh)+(u1*v2))+(u2*v1))+(u2*v2);        \
                              \
  ml += (xh)*(yl) + (xl)*(yh);                \
  *zh = mh+ml;                        \
  *zl = mh - (*zh) + ml;                              \
}

仅前8个运算用于将操作数中的每个双精度数精确地划分为两半，以便可以将每侧的一半乘以另一侧的一半，而得到的结果恰好为a double。计算u1*v1，，u1*v2…正是这样做的。

在mh和中获得的值ml可以重叠，因此最后3个操作用于将结果重新归一化为两个浮点数的总和。

在这种情况下，我们可以在每个双精度数上产生舍入误差，或者是避免这种误差的一种机制？

正如评论所说：

/*
 * computes double-double multiplication: zh+zl = (xh+xl) *  (yh+yl)
 * relative error is smaller than 2^-102
 */

您可以在《浮点算术手册》中找到用于获得这些结果的所有机制。