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

Question

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

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

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

Nicolas Brisebarre、Jean-Michel Muller 和 Saurabh Kumar Raina。当除数已知时，加速正确舍入浮点除法。 IEEE 计算机汇刊，卷。 53，第 8 期，2004 年 8 月，第 1069-1072 页。

论文作者提倡的技术将除数 y 的倒数预先计算为归一化的头尾对 z_h:z _l 如下： z_h = 1 / y, z_l = fma (-y, z_{h, 1) / 是的}。稍后，除法 q = x / y 然后计算为 q = fma (z_h, x, z_l * x )。本文推导了除数 y 必须满足的各种条件才能使该算法起作用。正如人们容易观察到的那样，当头尾符号不同时，该算法存在无穷大和零的问题。更重要的是，它无法为幅度非常小的股息 x 提供正确的结果，因为计算商尾 z_l * x em>，遭受下溢。

该论文还顺便提及了一种基于 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，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

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

    q = x / y

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

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

    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 的0.5%，具体取决于y)。

Brisebarre-Muller-Raina 文章表明，朴素方法的最大误差为 ±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 方法产生偶数、非二次除数 y 的结果不正确。 （或者，16777216/rev(x) 表示上限。）

添加于 2016 年 2 月 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 失败案例（但我还没有充分测试大于 8388608 的除数）上产生一个在 ±1 以内的正确误差估计。最后的除法应该不会报告假零，但我不能保证（还）。它没有考虑具有下溢问题的非常大的除数（例如 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）的除数，如果除数除以最小幂时两个产生的商小于 16777216 的数量级，并且商是奇数。

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 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 中。
  • 如果 succ(M) 存在于 U 中，则这对值将 succ(M) - M 贡献给结果的正值。
  • 如果 V 中存在 -succ(M)，则这对值会将负值 M - succ(M) 贡献给结果的负值。
  • 如果 pred(M) 存在于 U 中，则这对值会将负值 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 sup> 的大小，并使用此信息来应用 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 3

我喜欢 @Pascal 的回答，但在优化中，最好有一个简单且易于理解的转换子集，而不是一个完美的解决方案。

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

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

x / 2ⁿ

这与程序中的常量（以及所有可能的以 10 为底的分数）形成对比，它们是以下形式的有理数：

x / (2ⁿ * 5^m)

因此，一种优化将简单地测试 m == 0 的输入和倒数，因为这些数字以 FP 格式精确表示，并且对它们的操作应该产生在格式内准确的数字.

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

.25 .50 .75

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

Answer 4

浮点除法的结果是：

• 标志标志
• 一个有效数字
• 指数
• 一组标志（上溢、下溢、不精确等 - 请参阅 fenv()）

让前 3 件正确（但标志集不正确）是不够的。如果没有进一步的知识（例如，结果的哪些部分实际上很重要，被除数的可能值等），我会假设用乘以常数（和/或复杂的 FMA 混乱）代替常数除法几乎是永远不安全。

另外;对于现代 CPU，我也不认为用 2 个 FMA 替换一个分区总是一种改进。例如，如果瓶颈是指令获取/解码，那么这种“优化”会使性能变差。再举一个例子，如果后续指令不依赖于结果（CPU 可以在等待结果的同时并行执行许多其他指令），FMA 版本可能会引入多个依赖关系停顿并降低性能。对于第三个示例，如果所有寄存器都在使用，那么 FMA 版本（需要额外的“实时”变量）可能会增加“溢出”并降低性能。

请注意，（在许多但并非所有情况下）除以 2 的常数倍数或乘以 2 的常数倍数可以仅通过加法来完成（具体而言，将移位计数添加到指数）。