如何细化浮点除法结果？

Question

我有一个使用 newton-raphson 算法计算浮点平方根除法的算法。我的结果并不完全准确，有时会相差 1 ulp。

我想知道是否有用于浮点除法的优化算法来获得最终的精度。我对平方根使用塔克曼检验，但是否有类似的除法算法？还是 tuckerman 测试可以适用于除法？

我也试过使用这个算法，但没有得到完全准确的结果：

z= divisor
r_temp = divisor*q
 r = dividend - r_temp
result_temp = r*z
q + result_temp

Answer 1

正确舍入迭代除法结果的一种实用方法是生成一个初步商，使其在数学结果的一个 ulp 以内，然后使用精确计算的残差来计算最终结果。

精确残差计算的首选工具是融合乘加 (FMA) 运算。这种方法的大部分基础工作（无论是在数学方面还是在实际实现方面）都归功于 Peter Markstein，后来由其他研究人员进行了改进。马克斯坦的结果在他的书中得到了很好的总结：

Peter Markstein，IA-64 和初等函数：速度和精度。 Prentice-Hall 2000.

使用 Markstein 方法进行正确舍入除法的一种直接方法是首先计算正确舍入的倒数，然后通过将其与股息，然后是最终的基于残差的舍入步骤。

残差可用于直接计算最终的舍入结果，~~如下面代码中的商舍入所示~~ （我注意到此代码序列导致错误四舍五入的结果是 10¹¹ 个除法中的一个，并将其替换为 comparison-and-select 习语的另一个实例），这是 Markstein 使用的技术。或者，它可以用作有点类似于 Tuckerman 舍入的双边比较和 select 过程的一部分，这在下面的代码中显示为倒数舍入。

关于倒数计算有一个警告。许多常用的迭代方法（包括我在下面使用的方法）与 Markstein 的舍入技术结合使用时，如果除数的尾数完全由 1 位组成，则会产生不正确的结果。

解决这个问题的一种方法是特殊处理这种情况。在下面的代码中，我选择了双边比较和 select 方法，这也允许在舍入之前误差略大于 1 ulp，从而消除了在倒数迭代本身中使用 FMA 的需要。

请注意，我在下面的 C 代码中省略了对次正常结果的处理，以保持代码简洁易懂。我将自己限制在标准 C 库函数上，以完成诸如提取部分浮点操作数、组装浮点数以及应用 1 ulp 递增和递减等任务。大多数平台将提供具有更高性能的特定于机器的选项。

float my_divf (float a, float b)
{
    float q, r, ma, mb, e, s, t;
    int ia, ib;

    if (!isnanf (a+b) && !isinff (a) && !isinff (b) && (b != 0.0f)) {
        /* normal cases: remove sign, split args into exponent and mantissa */
        ma = frexpf (fabsf (a), &ia);
        mb = frexpf (fabsf (b), &ib);
        /* minimax polynomial approximation to 1/mb for mb in [0.5,1) */
        r =        - 3.54939341e+0f;
        r = r * mb + 1.06481802e+1f;
        r = r * mb - 1.17573657e+1f;
        r = r * mb + 5.65684575e+0f;
        /* apply one iteration with cubic convergence */
        e = 1.0f - mb * r;
        e = e * e + e;
        r = e * r + r;
        /* round reciprocal to nearest-or-even */
        e = fmaf (-mb, r, 1.0f); // residual of 1st candidate
        s = nextafterf (r, copysignf (2.0f, e)); // bump or dent 
        t = fmaf (-mb, s, 1.0f); // residual of 2nd candidate
        r = (fabsf (e) < fabsf (t)) ? r : s; // candidate with smaller residual
        /* compute preliminary quotient from correctly-rounded reciprocal */
        q = ma * r;
        /* round quotient to nearest-or-even */
        e = fmaf (-mb, q, ma); // residual of 1st candidate
        s = nextafterf (q, copysignf (2.0f, e)); // bump or dent 
        t = fmaf (-mb, s, ma); // residual of 2nd candidate
        q = (fabsf (e) < fabsf (t)) ? q : s; // candidate with smaller residual
        /* scale back into result range */
        r = ldexpf (q, ia - ib);
        if (r < 1.17549435e-38f) {
            /* sub-normal result, left as an exercise for the reader */
        }
        /* merge in sign of quotient */
        r = copysignf (r, a * b);
    } else {
        /* handle special cases */
        if (isnanf (a) || isnanf (b)) {
            r = a + b;
        } else if (b == 0.0f) {
            r = (a == 0.0f) ? (0.0f / 0.0f) : copysignf (1.0f / 0.0f, a * b);
        } else if (isinff (b)) {
            r = (isinff (a)) ? (0.0f / 0.0f) : copysignf (0.0f, a * b);
        } else {
            r = a * b;
        }
    }
    return r;
}

如何细化浮点除法结果？

How to refine the result of a floating point division result?

math

floating-point

division

ieee-754