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

Posted 2023-02-16

【中文标题】如何细化浮点除法的结果？

【英文标题】：How to refine the result of a floating point division result?

【发布时间】：2015-02-21 00:25:14

【问题描述】：

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

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

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

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

【问题讨论】：

有一种算法可以基于 FMA（融合乘加）对浮点除法的结果进行舍入，但根据您之前的问题，我认为这对您没有用，就像您的平台一样不提供FMA？ This technical report Alan H. Karp 和 Peter Markstein 于 1993 年指出：“按照最初的公式，Tuckerman 舍入只能用于平方根，不能用于除法。”

您是否尝试过文献检索，例如通过谷歌学术？以下论文可能有用，因为它不假定 FMA 的可用性：Bogdan Pasca, Correctly Rounded Floating-Point Division for DSP-enabled FPGAS。

【参考方案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;

【讨论】：

由于我没有 FMA，是否仍然可以获得正确的结果？在此过程中使用四舍五入或类似的东西。

@starbox：你可能想看看我之前提到的那篇论文，看看这种方法是否适合你。我不确定您的背景是什么：如果您在特定平台的限制下需要此功能用于生产软件或硬件 [您尚未详细解释]，您可能需要考虑聘请顾问。如果你只是想学习各种算法，我建议你拉论文并尝试实现。仔细研究细节是掌握相关技术的可靠方法。

谢谢，是的，我这里有那本书的副本。我正在创建没有 FMA 的自定义硬件，需要完全准确，所以我只是想知道我是否可以从您的代码中删除 FMA 指令并使用舍入到零来获得完全准确。我会尝试实验

@starbox：各种微处理器使用迭代近似进行除法（例如使用 Goldschmidt）并实现了正确的舍入。这需要一些额外的位才能正确舍入。例如，请参阅this paper 关于 Athlon 处理器中的除法。我记得你不赞成额外的位，从我从你过去的问题中回忆起。我不知道针对您当前问题的解决方案，它既不需要 FMA 也不需要额外的位，除了在整数算术中模拟过程。

我用这两个输入试过你的程序，我得到了错误的答案：divide = 0x3F7C15B4; divisor = 0xDE50F303;，当它应该是 0xA09A6CA1 时，答案是 0xA09A6CA2，有什么想法可能是错的吗？