c - 具有恒定整数除数的高效浮点除法

Question

近期 question ，是否允许编译器用浮点乘法代替浮点除法，启发我提出这个问题。

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

然后人们想知道这种替换的其他除数是否有效，假设我们允许任何替换除法但运行速度明显更快的短指令序列，同时提供位相同的结果。除了普通乘法之外，特别允许融合乘加运算。
在评论中，我指出了以下相关论文:

尼古拉斯·布里斯巴尔、让-米歇尔·穆勒和索拉布·库马尔·雷纳。当除数已知时，加速正确舍入的浮点除法。 IEEE 计算机交易，卷。 53，第 8 期，2004 年 8 月，第 1069-1072 页。

论文作者提倡的技术将除数 y 的倒数预先计算为归一化的头尾对 zh:zl 如下:zh = 1/y, zl = fma (-y, zh, 1)/y。之后，除法 q = x/y 然后计算为 q = fma (zh, x, zl * x)。该论文推导出除数 y 必须满足的各种条件才能使该算法工作。正如人们很容易观察到的那样，当头部和尾部的符号不同时，该算法存在无穷大和零的问题。更重要的是，它无法为数量级非常小的红利 x 提供正确的结果，因为商尾 zl * x 的计算存在下溢。

该论文还顺便提及了一种替代的基于 FMA 的除法算法，该算法由 Peter Markstein 在 IBM 时开创。相关引用是:

P.W.马克斯坦。在 IBM RISC System/6000 处理器上计算基本函数。 IBM 研究与开发杂志，卷。 34，第 1 期，1990 年 1 月，第 111-119 页

在 Markstein 的算法中，首先计算倒数 rc，从中形成初始商 q = x * rc。然后，使用 FMA 精确计算除法的余数，如 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);
}

在大多数处理器架构上，这应该转换为无分支指令序列，使用谓词、条件移动或选择类型指令。举个具体的例子:对于除以 3.0f , nvcc CUDA 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 的幂。当然，轶事证据，不是证据。

什么数学条件可以先验地确定除法转换为上述代码序列是否安全？ 答案可以假设所有的浮点运算都是在“舍入到最接近或偶数”的默认舍入模式下执行的。

#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

让我第三次重新开始。我们正在努力加速

    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 方法同样对几乎所有的两个非幂都失败了 y ，但要少得多 x产生不正确的结果(不到所有可能的一半 x ，取决于 y )。
Brisebarre-Muller-Raina 文章表明，naive 方法的最大误差为 ±1.5 ULP。
Markstein 方法为二的幂产生正确的结果 y ，以及奇数 y . (我没有发现 Markstein 方法的失败奇整数除数。)

对于 Markstein 方法，我分析了除数 1 - 19700 ( raw data here )。
绘制失败案例的数量(水平轴中的除数， x 的值数量，其中 Markstein 方法对于所述除数失败)，我们可以看到一个简单的模式发生:

(来源: nominal-animal.net)
请注意，这些图具有水平和垂直轴对数。奇数除数没有点，因为该方法为我测试过的所有奇数除数产生正确结果。
如果我们将 x 轴更改为除数的位反转(相反顺序的二进制数字，即 0b11101101 → 0b10110111, data )，我们有一个非常清晰的模式:

(来源: nominal-animal.net)
如果我们通过点集的中心画一条直线，我们得到曲线 4194304/x . (请记住，该图仅考虑可能的浮点数的一半，因此在考虑所有可能的浮点数时，将其加倍。) 8388608/x和 2097152/x将整个错误模式完全括起来。
因此，如果我们使用 rev(y)计算除数的位反转 y ，然后 8388608/rev(y)是案例数量的一个很好的一阶近似值(在所有可能的浮点数中)，其中 Markstein 方法为偶数、非 2 的幂除数 y 产生不正确的结果. (或者， 16777216/rev(x) 作为上限。)
2016 年 2 月 28 日添加:我使用 Markstein 方法找到了错误案例数量的近似值，给定任何整数(binary32)除数。这是伪代码:

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。 (在每一对中，尾数是相同的，只有两个的幂不同。)
假设我的测试台没有错误，这意味着 Markstein 方法也适用于数量级大于 16777216 的除数(但小于，比如说，1e+30)，如果除数是这样的，当除以 2 的最小幂时产生的商在数量上小于 16777216，并且商是奇数。