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

Question

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

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

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

Nicolas Brisebarre、Jean-Michel Muller和Saurabh Kumar Raina。当除数事先已知时，加速正确舍入的浮点除法。IEEE Transaction on Compux，Vol.53，No.8,2004年8月，第1069-1072页。

本文作者提倡的技术将除数y的倒数预先计算为一个标准化的头尾对z_{h：zl，如下所示：zh=1/y，zl=fma（-y，zh，1）/y。然后，除数q=x/y被计算为q=fma（zh，x，zl*x）。本文推导了除数y必须满足的各种条件。正如人们容易观察到的那样，当头部和尾部的符号不同时，该算法存在无穷大和零的问题。更重要的是，它将无法为数量非常小的股息x提供正确的结果，因为商尾zl*x的计算受到下溢的影响。}

本文还顺便提到了彼得·马克斯坦（PeterMarkstein）在IBM时首创的另一种基于FMA的除法算法。相关参考是：

P、 W.Markstein。在IBM RISC System/6000处理器上计算基本函数。IBM研究杂志

在Markstein的算法中，首先计算一个倒数rc，由此形成初始商q=x*rc。然后，用FMA精确计算除法的剩余部分，即r=FMA（-y，q，x），并最终计算更精确的改进商，即q=FMA（r，rc，q）。

该算法还存在x为零或无穷大的问题（通过适当的条件执行很容易解决），但使用IEEE-754单精度浮点数据进行的详尽测试表明，它在所有可能的股息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的除法，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的所有可能编码中正确工作。轶事证据，当然不是证明。

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

#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

我喜欢@Pascal的答案，但在优化中，拥有一个简单且易于理解的转换子集通常比完美的解决方案更好。

所有当前和常见的历史浮点格式都有一个共同点：二进制尾数。

因此，所有分数都是形式上的有理数：

x/2n

这与程序中的常量（以及所有可能的10进制分数）形成对比，这些常量是形式上的有理数：

x/（2n*5m）

因此，一个优化只需测试m==0的输入和倒数，因为这些数字以FP格式精确表示，使用它们的操作应生成格式内准确的数字。

因此，例如，在0.01到0.99的（十进制2位数）范围内，将优化除以或乘以以下数字：

.25 .50 .75

其他一切都不会。（我想，先测试一下，哈哈。）

Answer 2

这个问题要求一种方法来识别常量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^的数字-24在幅度上，但不是-0.0f或任何比2^-24更接近零的值。这证明了将常量Y的许多值转换为问题中显示的两种形式之一是合理的。

（*）假设没有它，许多优化是不可能的，并且C编译器已经做出了假设，除非程序明确使用#pragma STDCFENV_ACCESSON

预测上述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中。
  • 如果ucc（M）存在于U中，则这对值对结果的正值有贡献。
  • 如果-suc（M）存在于V中，则这对值将负值M-suc（M）贡献给结果的负值。
  • 如果pred（M）存在于U中，则这对值将负值pred（M）-M贡献给结果的负值。
  • 如果-pred（M）存在于V中，则这对值将值M-pred（M）贡献给结果的正值。

  承认：以上借鉴了Bruno Marre“改进浮点加减约束”的思想

  示例：以下函数的编译：

  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。

  （**）“摆弄符号”是指当x只能是两个有符号零中的一个时，为了使结果正确显示，需要使用-FMA（-C1，x，（-C2）*x）代替FMA（C1，x，C2*x）

Answer 3

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

    q = x / y

其中，y是一个整数常量，q、x和y都是IEEE 754-2008二进制32浮点值。下面，fmaf（a，b，c）表示使用二进制32值的融合乘法加法。

朴素算法是通过预先计算的倒数，

    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中的一半以上，它会产生错误的结果。

布里斯巴雷-穆勒-雷纳方法在几乎所有非幂次的两个y中也同样失败，但产生错误结果的x要少得多（在所有可能的x中，不到0.5%，因y而不同）。

Brisebarre-Muller-Raina的文章表明，朴素方法的最大误差为±1.5 ULPs。

Markstein方法对2的幂产生了正确的结果，也对奇数整数产生了正确的结果。（对于Markstein方法，我没有找到失败的奇数整数除数。）

对于Markstein方法，我分析了除数1-19700（这里的原始数据）。

绘制失效案例的数量（横轴上的除数，Markstein方法对所述除数失效时x的值的数量），我们可以看到一个简单的模式：

请注意，这些图的水平轴和垂直轴都是对数轴。奇数除数没有点，因为这种方法可以为我测试的所有奇数除数生成正确的结果。

如果我们通过点集的中心画一条直线，我们就会得到曲线。（请记住，绘图只考虑了一半可能的浮动，因此在考虑所有可能的浮动时，将其加倍。）<代码>8388608/x和2097152/x完全包含整个错误模式。

因此，如果我们使用除数的位反转来计算除数，那么8388608/rev（y）是一个很好的一阶近似值，可以计算出Markstein方法对两个除数的偶数、非幂次产生错误结果的情况数（在所有可能的浮点数之外）。（或，<代码>16777216/修订版（x） 用于上限。）

新增2016-02-28：我使用Markstein方法，在给定任何整数（二进制32）除数的情况下，找到了错误案例数的近似值。此处为伪代码：

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的除数，我始终看到与除以二的最小幂得出小于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）的除数，如果除数被二的最小幂除时，得到数量级小于16777216的商，并且商是奇数。