具有常数整数除数的高效浮点除法

Question

最近的一个问题是,编译器是否允许用浮点乘法替换浮点除法,这激发了我提出这个问题.

在严格的要求下,代码转换后的结果应与实际的除法运算在位上相同,可以看出,对于二进制IEEE-754算术,这对于2的幂的除数是可能的.只要除数的倒数可以表示,乘以除数的倒数就可以得到与除数相同的结果.例如,乘法0.5可以代替除法2.0.

然后人们想知道其他除数这样的替换是什么工作,假设我们允许任何短指令序列取代除法但运行速度明显更快,同时提供相同的结果.特别是除了普通乘法之外,还允许融合乘法 - 加法运算.在评论中,我指出了以下相关文件:

Nicolas Brisebarre,Jean-Michel Muller和Saurabh Kumar Raina.当事先知道除数时,加速正确舍入的浮点除法.IEEE Transactions on Computers,Vol.53,第8期,2004年8月,第1069-1072页.

本文作者提出的技术预先计算了除数y的倒数作为归一化的头尾对z _h:z _l如下:z _h = 1/y,z _l = fma(-y,z _h,1 )/ y.之后,将除数q = x/y计算为q = fma(z _h,x,z _l*x).本文推导出除数y必须满足的各种条件才能使该算法起作用.正如人们容易观察到的那样,当头尾迹象不同时,该算法存在无穷大和零的问题.更重要的是,它将无法为数量非常小的股息x提供正确的结果,因为商尾的计算z _l*x受到下溢的影响.

本文还提到了另一种基于FMA的划分算法,该算法由Peter Markstein在IBM工作时开创.相关参考是:

PW Markstein.在IBM RISC System/6000处理器上计算基本功能.IBM Journal of Research&Development,Vol.1990年1月34日第1号,第111-119页

在Markstein算法中,首先计算倒数rc,从中形成初始商q = x*rc.然后,用FMA精确地计算除法的余数r = fma(-y,q,x),并且最终计算出改进的,更准确的商,其为q = fma(r,rc,q).

该算法还存在x为零或无穷的问题(可以通过适当的条件执行轻松解决),但使用IEEE-754单精度float数据的详尽测试表明,它为所有可能的红利x提供了正确的商数,对于许多除数y,在这些许多小整数中.这个C代码实现了它:

/* precompute reciprocal */
rc = 1.0f / y;

/* compute quotient q=x/y */
q = x * rc;
if ((x != 0) && (!isinf(x))) {
    r = fmaf (-y, q, x);
    q = fmaf (r, rc, q);
}

在大多数处理器体系结构中,这应该转换为无分支指令序列,使用预测,条件移动或选择类型指令.举一个具体的例子:对于除法3.0f,nvccCUDA 7.5 的编译器为Kepler级GPU生成以下机器代码:

    LDG.E R5, [R2];                        // load x
    FSETP.NEU.AND P0, PT, |R5|, +INF , PT; // pred0 = fabsf(x) != INF
    FMUL32I R2, R5, 0.3333333432674408;    // q = x * (1.0f/3.0f)
    FSETP.NEU.AND P0, PT, R5, RZ, P0;      // pred0 = (x != 0.0f) && (fabsf(x) != INF)
    FMA R5, R2, -3, R5;                    // r = fmaf (q, -3.0f, x);
    MOV R4, R2                             // q
@P0 FFMA R4, R5, c[0x2][0x0], R2;          // if (pred0) q = fmaf (r, (1.0f/3.0f), q)
    ST.E [R6], R4;                         // store q

对于我的实验,我编写了下面显示的微小C测试程序,它按递增顺序逐步执行整数除数,并且每个都对前面的代码序列进行详尽的测试,以防止正确的除法.它打印出通过此详尽测试的除数列表.部分输出如下:

PASS: 1, 2, 3, 4, 5, 7, 8, 9, 11, 13, 15, 16, 17, 19, 21, 23, 25, 27, 29, 31, 32, 33, 35, 37, 39, 41, 43, 45, 47, 49, 51, 53, 55, 57, 59, 61, 63, 64, 65, 67, 69,

为了将替换算法作为优化结合到编译器中,可以安全地应用上述代码转换的除数白名单是不切实际的.到目前为止,程序的输出(以每分钟大约一个结果的速率)表明快速代码x对于y奇数整数或2的幂的除数的所有可能编码都能正常工作.轶事证据,当然不是证明.

什么样的数学条件可以先验地确定划分为上述代码序列的转换是否安全？答案可以假设所有浮点运算都在默认的舍入模式"round to nearest or even"中执行.

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

int main (void)
{
    float r, q, x, y, rc;
    volatile union {
        float f;
        unsigned int i;
    } arg, res, ref;
    int err;

    y = 1.0f;
    printf ("PASS: ");
    while (1) {
        /* precompute reciprocal */
        rc = 1.0f / y;

        arg.i = 0x80000000;
        err = 0;
        do {
            /* do the division, fast */
            x = arg.f;
            q = x * rc;
            if ((x != 0) && (!isinf(x))) {
                r = fmaf (-y, q, x);
                q = fmaf (r, rc, q);
            }
            res.f = q;
            /* compute the reference, slowly */
            ref.f = x / y;

            if (res.i != ref.i) {
                err = 1;
                break;
            }
            arg.i--;
        } while (arg.i != 0x80000000);

        if (!err) printf ("%g, ", y);
        y += 1.0f;
    }
    return EXIT_SUCCESS;
}

Answer 1

这个问题要求一种方法来识别常量的值,Y这样可以安全地x / Y使用FMA 转换为更便宜的计算,用于所有可能的值x.另一种方法是使用静态分析来确定值x可以采用的过度近似,从而可以应用通常不健全的变换,因为知道变换的代码与原始分区不同的值不会发生.

使用很好地适应浮点计算问题的浮点值集合的表示,即使从函数开头开始的向前分析也可以产生有用的信息.例如:

float f(float z) {
  float x = 1.0f + z;
  float r = x / Y;
  return r;
}

假设默认的舍入到最近模式(*),在上面的函数x中只能是NaN(如果输入是NaN),+ 0.0f,或者数量大于2到²⁴的数字,但不是-0.0f或任何接近零的东西都不是2到²⁴.这证明了对于常数的许多值,转换为问题中所示的两种形式之一Y.

(*)假设没有这种假设,许多优化是不可能的,并且C编译器已经制作,除非程序明确使用 #pragma STDC FENV_ACCESS ON

预测上述信息的转发静态分析x可以基于表达式可以作为以下元组的浮点值集合的表示:

• 可能的NaN值集合的表示(由于NaN的行为未被指定,选择是仅使用布尔值,true意味着可以存在一些NaN,并且false表示不存在NaN.),
• 四个布尔标志分别表示+ inf,-inf,+0.0,-0.0的存在,
• 负有限浮点值的包含间隔,和
• 正有限浮点值的包含间隔.

为了遵循这种方法,静态分析器必须理解C程序中可能发生的所有浮点运算.为了说明,用于+在分析的代码中处理的值集U和V之间的相加可以实现为:

• 如果NaN存在于其中一个操作数中,或者如果操作数可以是相反符号的无穷大,则结果中存在NaN.
• 如果0不能是添加U值和V值的结果,则使用标准区间运算.对于U中最大值和V中最大值的舍入到最近相加,得到结果的上界,因此这些边界应该用舍入到最近计算.
• 如果0可以是添加U的正值和V的负值的结果,那么令M是U中的最小正值,使得-M存在于V.
  • 如果u中存在succ(M),那么这对值将succ(M)-M贡献给结果的正值.
  • 如果-succ(M)存在于V中,那么这对值将负值M-succ(M)贡献给结果的负值.
  • 如果在U中存在pred(M),那么这对值将负值pred(M)-M贡献给结果的负值.
  • 如果-pred(M)存在于V中,那么这对值将值M-pred(M)贡献给结果的正值.
• 如果0可以是添加负值U和正值V的结果,则执行相同的工作.

致谢:上述借鉴了"改善浮点加法和减法约束",Bruno Marre和Claude Michel

示例:编译以下函数f:

float f(float z, float t) {
  float x = 1.0f + z;
  if (x + t == 0.0f) {
    float r = x / 6.0f;
    return r;
  }
  return 0.0f;
}

问题中的方法拒绝将函数中的除法f转换为替代形式,因为6不是可以无条件转换除法的值之一.相反,我建议的是从函数的开头应用一个简单的值分析,在这种情况下,确定它x是一个有限的浮点数+0.0f或至少2到^24个数量级,并使用此信息来应用Brisebarre等人的转型,对x * C2不会下流的知识充满信心.

为了明确,我建议使用如下所示的算法来决定是否将分割转换为更简单的分类:

1. 是否Y可以根据他们的算法使用Brisebarre等人的方法转换其中一个值？
2. 他们的方法中的C1和C2是否具有相同的符号,或者是否可以排除红利无限的可能性？
3. 他们的方法中的C1和C2是否具有相同的符号,或者只能x采用0的两个表示中的一个？如果在C1和C2具有不同符号并且x只能是零的一个表示的情况下,记住用基于FMA的计算的符号来调整(**)以使其在零时产生正确的x零.
4. 可以保证股息的大小足以排除下跌的可能性x * C2吗？

如果对四个问题的答案为"是",则可以在正在编译的函数的上下文中将除法转换为乘法和FMA.上述静态分析用于回答问题2,3和4.

(**)"摆弄标志"是指使用-FMA(-C1,x,( - C2)*x)代替FMA(C1,x,C2*x),这样才能使结果出来当x只能是两个带符号的零之一时正确

Answer 2

让我第三次重启.我们正在努力加速

    q = x / y

其中y为整数常数,并且q,x和y都是IEEE 754-2008 binary32浮点值.下面,fmaf(a,b,c)表示a * b + c使用binary32值的融合乘法 - 加法.

天真的算法是通过预先计算的倒数,

    C = 1.0f / y

这样在运行时一个(快得多)乘法就足够了:

    q = x * C

Brisebarre-Muller-Raina加速度使用两个预先计算的常数,

    zh = 1.0f / y
    zl = -fmaf(zh, y, -1.0f) / y

这样在运行时,一个乘法和一个融合乘法 - 加法就足够了:

    q = fmaf(x, zh, x * zl)

Markstein算法将朴素方法与两个融合乘法相加结合起来,如果天真方法在最不重要的位置产生1个单位内的结果,通过预先计算得到正确的结果

    C1 = 1.0f / y
    C2 = -y

因此,可以使用近似来区分

    t1 = x * C1
    t2 = fmaf(C1, t1, x)
    q  = fmaf(C2, t2, t1)

天真的方法适用于所有两种力量y,但除此之外它非常糟糕.例如,对于除数7,14,15,28和30,它会产生超过一半可能的不正确结果x.

Brisebarre-Muller-Raina方法同样对几乎所有非2的幂都失败了y,但x产生的错误结果却少得多(不到所有可能的半数x,取决于y).

Brisebarre-Muller-Raina文章显示,天真方法的最大误差为±1.5 ULPs.

Markstein方法得到2的幂的正确结果y,也对于奇数整数y.(我没有找到Markstein方法的失败奇整数除数.)

对于Markstein方法,我已经分析了除数1 - 19700(原始数据).

绘制失败案例的数量(横轴的除数x,Markstein逼近除数的值的数量),我们可以看到一个简单的模式:

Markstein失败案例http://www.nominal-animal.net/answers/markstein.png

请注意,这些图具有水平轴和垂直轴对数.奇数除数没有点,因为该方法可以为我测试的所有奇数除数产生正确的结果.

如果我们将x轴更改为除数的位反转(反向二进制数字,即0b11101101→0b10110111,数据),我们有一个非常清晰的模式: Markstein失败案例,位反转除数http://www.nominal- animal.net/answers/markstein-failures.png

如果我们在点集的中心绘制一条直线,我们得到曲线4194304/x.(请记住,该图只考虑了一半可能的浮点数,因此在考虑所有可能的浮点数时,加倍它.) 8388608/x并2097152/x完全包含整个错误模式.

因此,如果我们rev(y)用来计算除数的位反转y,那么 8388608/rev(y)就是一个良好的一阶近似的案例数(在所有可能的浮点数中),其中Markstein方法产生一个偶数,非幂的不正确的结果两个除数y.(或者,16777216/rev(x)上限.)

添加2016-02-28:在给定任何整数(binary32)除数的情况下,我找到了使用Markstein方法的错误情况数的近似值.这是伪代码:

function markstein_failure_estimate(divisor):
    if (divisor is zero)
        return no estimate
    if (divisor is not an integer)
        return no estimate

    if (divisor is negative)
        negate divisor

    # Consider, for avoiding underflow cases,
    if (divisor is very large, say 1e+30 or larger)
        return no estimate - do as division

    while (divisor > 16777216)
        divisor = divisor / 2

    if (divisor is a power of two)
        return 0

    if (divisor is odd)
        return 0

    while (divisor is not odd)
        divisor = divisor / 2

    # Use return (1 + 83833608 / divisor) / 2
    # if only nonnegative finite float divisors are counted!
    return 1 + 8388608 / divisor

在我测试的Markstein失效案例中,这产生了一个正确的误差估计值±1(但我还没有充分测试大于8388608的除数).最终的划分应该是它没有报告错误的零,但我不能保证它(还).它没有考虑具有下溢问题的非常大的除数(比如0x1p100,或1e + 30,并且幅度更大) - 无论如何我绝对会将这些除数从加速中排除.

在初步测试中,估计似乎非常准确.我没有绘制比较估计值和除数1到20000的实际误差的图,因为这些点在图中完全重合.(在此范围内,估计值是精确的,或者太大.)基本上,估计值会准确地再现此答案中的第一个图.

Markstein方法的失败模式是规则的,非常有趣.该方法适用于两个除数的所有幂和所有奇数整数除数.

对于大于16777216的除数,我一直看到与除数相同的误差,除数除以2的最小幂,得到小于16777216的值.例如,0x1.3cdfa4p + 23和0x1.3cdfa4p + 41,0x1. d8874p + 23和0x1.d8874p + 32,0x1.cf84f8p + 23和0x1.cf84f8p + 34,0x1.e4a7fp + 23和0x1.e4a7fp + 37.(在每对中,尾数是相同的,只有2的幂变化.)

假设我的测试平台没有错误,这意味着Markstein方法在大小上也会使用大于16777216的除数(但是比1e + 30小),如果除数是这样的除数,除以2的最小幂.产生的数量小于16777216,并且商是奇数.