如何细化浮点除法结果? [英] How to refine the result of a floating point division result?

问题描述

我有一个使用牛顿-拉普森算法来计算浮点数平方根除法的算法.我的结果并不完全准确，有时会偏离1 ulp.

I have an an algorithm for calculating the floating point square root divide using the newton-raphson algorith. My results are not fully accurate and sometimes off by 1 ulp.

我想知道是否有一种用于浮点除法的精炼算法来获得最终的精度.我将tuckerman检验用于平方根，但是有类似的除法算法吗?还是可以将塔克曼测验进行分组划分?

I was wondering if there is a refinement algorithm for floating point division to get the final bits of accuracy. I use the tuckerman test for square root, but is there a similar algorithm for division? Or can the tuckerman test be adapted for division?

我也尝试使用此算法，但未获得完全准确的结果:

I tried using this algorithm too but didn't get full accuracy results:

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

推荐答案

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

One practical way of correctly rounding the result of iterative division is to produce a preliminary quotient to within one ulp of the mathematical result, then use the exactly-computed residual to compute the final result.

精确计算残差的首选工具是融合乘加(FMA)操作.这种方法的许多基础工作(无论是在数学上还是在实际实现方面)都归功于彼得·马克斯坦(Peter Markstein)，后来又由其他研究人员进行了完善.Markstein的结果很好地总结在他的书中:

The tool of choice for the exact computation of residuals is the fused-multiply add (FMA) operation. Much of the foundational work of this approach (both in terms of the mathematics and of practical implementations) is due to Peter Markstein and was later refined by other researchers. Markstein's results are nicely summarized in his book:

Peter Markstein，IA-64和基本功能:速度和精度.普伦蒂斯厅2000.

使用马克斯坦(Markstein)方法进行正确舍入的直接方法是首先计算正确舍入的倒数，然后将其乘以乘积来计算正确舍入的商.股利，然后是基于残差的最终四舍五入步骤.

A straightforward approach to correctly-rounded division using Markstein's approach is to first compute a correctly-rounded reciprocal, then compute the correctly-rounded quotient by multiplying it with the dividend, followed by the final residual-based rounding step.

残差可用于直接计算最终的舍入结果，如下面代码中商数舍入所示(我注意到此代码序列导致一次舍入不正确的舍入结果)10 ¹¹ 除法，然后用比较选择成语的另一个实例替换它，这是Markstein使用的技术.另外，它也可以用作类似于Tuckerman舍入法的双向比较和选择过程的一部分，在下面的代码中为倒数舍入而显示.

The residual can be used to compute the final rounded result directly, as is shown for the quotient rounding in the code below (I noticed that this code sequence resulted in an incorrectly rounded result in one out of 10¹¹ divisions, and replaced it with another instance of the comparison-and-select idiom) which is the technique used by Markstein. Alternatively it may be used as part of a two-sided comparison-and-select process somewhat akin to Tuckerman rounding, which is shown for the reciprocal rounding in the code below.

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

There is one caveat with regard to the reciprocal computation. Many commonly used iterative approaches (including the one I used below), when combined with Markstein's rounding technique, deliver an incorrect result if the mantissa of the divisor consists entirely of 1-bits.

解决这个问题的一种方法是特别对待这种情况.在下面的代码中，我选择了一种双向比较和选择方法，该方法还允许在舍入前使误差略大于一个ulp，从而消除了在倒数迭代本身中使用FMA的需要.

One way of getting around this is to treat this case specially. In the code below I instead opted for a two-sided comparison-and-select approach, which also allows errors slightly larger than one ulp prior to rounding and thus eliminates the need to use FMA in the reciprocal iteration itself.

请注意，为了使代码简洁明了且易于理解，我在下面的C代码中省略了对次正规结果的处理.我将自己限制在标准C库函数中，以执行诸如提取浮点操作数的部分，组合浮点数以及应用一进制增量和减量之类的任务.大多数平台将为这些机器提供特定于机器的选项，并具有更高的性能.

Please note that I omitted the handling of sub-normal results in the C code below to keep the code concise and easy to follow. I have limited myself to standard C library functions for tasks like extracting parts of floating-point operands, assembling floating-point numbers, and applying one-ulp increments and decrements. Most platforms will offer machine-specific options with higher performance for these.

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;
}

这篇关于如何细化浮点除法结果?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！