C++ 使用常量整数除数的高效浮点除法

最近的question次会议，是否允许编译器用浮点乘法代替浮点除法，启发了我提出这个问题.

在严格要求代码转换后的结果必须与实际除法运算逐位相同的情况下， 很容易看出，对于二进制IEEE-754算术，这对于2的幂的除数是可能的.只要是互惠的 除数的值是可表示的，乘以除数的倒数可以得到与除法相同的结果.例如，乘以0.5可以替换除以2.0.

假设我们允许任何替换除法的短指令序列，但运行速度要快得多，同时提供位相同的结果，那么人们会想知道这样的替换会对其他除数起什么作用.具体地说，除了普通乘法之外，还允许融合乘法-加法运算. 我在 comments 中指出了以下有关文件:

Nicolas Brisebarre, Jean-Michel Muller, and Saurabh Kumar Raina. Accelerating correctly rounded floating-point division when the divisor is known in advance. IEEE Transactions on Computers, Vol. 53, No. 8, August 2004, pp. 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时首创.相关参考资料如下:

P. W. Markstein. Computation of elementary functions on the IBM RISC System/6000 processor. IBM Journal of Research & Development, Vol. 34, No. 1, January 1990, pp. 111-119个

在马克斯坦的算法中，人们首先计算倒数rc，由此形成初始商q = x * rc.然后，使用FMA AS r = fma (-y, q, x)精确计算除法的剩余部分，并最终将改进的、更精确的商计算为q = fma (r, rc, q).

该算法也存在x是零或无穷大的问题(通过适当的条件执行很容易解决)，但是使用IEEE-754单精度float数据的详尽测试表明，它在许多小整数中的许多除数y的所有可能的被除数x上提供了正确的商.此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);
}

在大多数处理器体系 struct 上，这应该转换为使用谓词、条件移动或 Select 类型指令的无分支指令序列.举一个具体的例子:对于除以3.0f，CUDA 7.5的nvcc编译器为开普勒类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,

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

What set of mathematical conditions can determine a-priori whether the transformation of division into the above code sequence is safe?个答案可以假设所有浮点运算都是在默认舍入模式"舍入到最近或偶数"下执行的.

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