如何细化浮点除法的结果？

Question

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

我想知道是否有用于浮点除法的改进算法来获得最终的精度。我对平方根使用 tuckerman 检验，但是否有类似的除法算法？或者 tuckerman 测试可以适用于除法吗？

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

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

Answer 1

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

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

Peter Markstein，IA-64 和基本函数：速度和精度。 Prentice-Hall 2000。

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

残差可以直接用来计算最终的舍入结果，如下面代码中的商舍入所示（我注意到这个代码序列导致了一个不正确的舍入结果10¹¹ 的划分，并将其替换为 Markstein 使用的技术的另一个比较和选择成语的实例。或者，它可以用作双边比较和选择过程的一部分，有点类似于 Tuckerman 舍入，下面的代码中显示了倒数舍入。

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

解决此问题的一种方法是对这种情况进行特殊处理。在下面的代码中，我选择了一种双边比较和选择方法，它还允许在舍入之前略大于 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;
}