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*v1
、u1*v2
、……正是这样做的。
mh
和ml
中得到的值可以重叠,所以最后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 中找到用于实现这些结果的所有机制。
我在维基百科上看到,实现四精度的某种方式是使用双双运算,即使它在位方面的精度不完全相同: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*v1
、u1*v2
、……正是这样做的。
mh
和ml
中得到的值可以重叠,所以最后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 中找到用于实现这些结果的所有机制。