float128 和 double-double 运算

float128 and double-double arithmetic

我在维基百科上看到,实现四精度的某种方式是使用双双运算,即使它在位方面的精度不完全相同:https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format

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

在这种情况下,我们可以对每个 double 进行舍入误差,或者他们有一种避免这种情况的机制?

“In this case, we use two double to store the value. So we need to make two operations at each time.”

这不是双双算术的工作原理。您应该期望在 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 个操作就是将操作数中的每个 double 精确地分成两半,以便每一侧的一半可以与另一侧的一半相乘,得到的结果正好是 double.计算 u1*v1u1*v2、……正是这样做的。

mhml中得到的值可以重叠,所以最后3次运算是为了将结果重新归一化为两个浮点数的和。

In this case we can have round-off errors on each double or their is a mechanism that avoid this?

正如评论所说:

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

您可以在 Handbook of Floating-Point Arithmetic 中找到用于实现这些结果的所有机制。