具有常数整数除数的高效浮点除法
Efficient floating-point division with constant integer divisors
最近的一个，是否允许编译器用浮点乘法代替浮点除法，启发了我问这个问题。
在严格要求下，代码转换后的结果必须与实际的除法运算按位完全相同，
很容易看出，对于二进制 IEEE-754 算法，这对于 2 的幂的除数是可能的。只要互惠
除数的可表示性，乘以除数的倒数得到与除法相同的结果。例如，乘以 0.5 可以代替除以 2.0.
然后人们想知道还有哪些其他除数可以替代除法，假设我们允许任何短的指令序列来替代除法但运行速度明显更快，同时提供位相同的结果。除了普通乘法之外，特别允许融合乘加运算。
在评论中，我指出了以下相关论文：
Nicolas Brisebarre、Jean-Michel Muller 和 Saurabh Kumar Raina。当除数事先已知时，加速正确舍入的浮点除法。 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);
}
在大多数处理器架构上，这应该转化为无分支的指令序列，使用预测、条件移动或 select 类型的指令。举个具体的例子：除以3.0f，CUDA 7.5的nvcc编译器为Kepler-class 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 正确工作两个。轶事证据，当然不是证据。
什么数学条件可以先验判断除法变换成上面的代码序列是否安全？答案可以假设所有的浮点运算都进行在 "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;
}
这个问题要求找到一种方法来识别常量 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 中最小的正值，使得 V 中存在 -M。
- 如果 U 中存在 succ(M)，则这对值将 succ(M) - M 贡献给结果的正值。
- 如果-succ(M)存在于V中，则这对值将负值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 大小，并使用此信息应用 Brisebarre 等人的变换，确信 x * C2 不会下溢。
明确地说，我建议使用如下算法来决定是否将除法转换为更简单的算法：
Y 是可以根据他们的算法使用 Brisebarre 等人的方法转换的值之一吗？- 他们方法中的C1和C2是同号吗，还是可以排除被除数是无穷大的可能？
- 他们的方法中的 C1 和 C2 是否具有相同的符号，或者
x 可以仅采用 0 的两种表示形式之一？如果在C1和C2符号不同，x只能是0的一种表示的情况下，记得fiddle(**)加上FMA计算的符号，使其产生当 x 为零时更正零。
- 能否保证分红的幅度足够大，排除
x * C2下溢的可能？
如果四个问题的答案都是“是”，那么除法可以在编译函数的上下文中转化为乘法和FMA。上述静态分析用于回答问题 2.、3. 和 4.
(**) “摆弄符号”意味着在需要时使用 -FMA(-C1, x, (-C2)*x) 代替 FMA(C1, x, C2*x)当 x 只能是两个带符号的零之一时，使结果正确出来
浮点除法结果为：
- 标志旗
- 有效数
- 指数
- 一组标志（上溢、下溢、不精确等 - 参见
fenv()）
仅正确地获得前 3 个部分（但标志集不正确）是不够的。在没有进一步知识的情况下（例如，结果的哪些部分实际上很重要，股息的可能值等）我会假设用常数乘法代替常数除法（and/or 一个令人费解的 FMA 混乱）几乎从来都不安全。
此外；对于现代 CPUs，我也不认为用 2 个 FMA 替换一个分区总是一种改进。例如，如果瓶颈是指令 fetch/decode，那么这个 "optimisation" 会使性能变差。再例如，如果后续指令不依赖于结果（CPU 可以在等待结果的同时并行执行许多其他指令），FMA 版本可能会引入多个依赖性停顿并使性能变差。对于第三个示例，如果所有寄存器都被使用，则 FMA 版本（需要额外的 "live" 变量）可能会增加 "spilling" 并使性能变差。
请注意（在许多但不是所有情况下）除以或乘以 2 的常数倍数可以单独通过加法来完成（具体来说，将移位计数添加到指数）。
我喜欢 @Pascal 的回答，但在优化中，拥有简单且易于理解的转换子集通常比拥有完美的解决方案更好。
所有当前和常见的历史浮点格式都有一个共同点：二进制尾数。
因此，所有分数都是以下形式的有理数：
x / 2n
这与程序中的常量（以及所有可能的以 10 为底的分数）形成对比，这些常量是以下形式的有理数：
x / (2n * 5m)
因此，一种优化将简单地测试 m == 0 的输入和倒数，因为这些数字完全以 FP 格式表示，并且对它们的操作应该产生以下数字在格式内是准确的。
因此，例如，在 .01 到 0.99 的（十进制 2 位数）范围内除以或乘以以下数字将得到优化：
.25 .50 .75
而其他一切都不会。 （我想，先测试一下，哈哈。）
让我第三次重启。我们正在努力加速
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 算法通过预先计算
C1 = 1.0f / y
C2 = -y
这样可以使用
来近似除法
t1 = x * C1
t2 = fmaf(C1, t1, x)
q = fmaf(C2, t2, t1)
朴素的方法适用于 2 y 的所有幂，但除此之外就很糟糕了。例如，对于除数 7、14、15、28 和 30，超过一半的可能 x.
会产生不正确的结果
Brisebarre-Muller-Raina 方法对于几乎所有 2 的非幂 y 都同样失败，但更少的 x 产生不正确的结果（不到所有可能 y 的一半 x，因 y) 而异。
Brisebarre-Muller-Raina 文章表明，朴素方法的最大误差为 ±1.5 ULP。
Markstein 方法对 2 的幂 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) 是案例数的一个很好的一阶近似值（在所有可能的情况下float)，其中 Markstein 方法对偶数、非二次方除数 y 产生错误结果。 （或者，16777216/rev(x) 为上限。）
已添加 2016-02-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）的除数，如果除数是这样的，当除以商小于 16777216 的最小二乘方，商为奇数。
最近的一个
在严格要求下，代码转换后的结果必须与实际的除法运算按位完全相同，
很容易看出，对于二进制 IEEE-754 算法，这对于 2 的幂的除数是可能的。只要互惠
除数的可表示性，乘以除数的倒数得到与除法相同的结果。例如，乘以 0.5 可以代替除以 2.0.
然后人们想知道还有哪些其他除数可以替代除法，假设我们允许任何短的指令序列来替代除法但运行速度明显更快，同时提供位相同的结果。除了普通乘法之外，特别允许融合乘加运算。 在评论中，我指出了以下相关论文：
Nicolas Brisebarre、Jean-Michel Muller 和 Saurabh Kumar Raina。当除数事先已知时，加速正确舍入的浮点除法。 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);
}
在大多数处理器架构上，这应该转化为无分支的指令序列，使用预测、条件移动或 select 类型的指令。举个具体的例子：除以3.0f，CUDA 7.5的nvcc编译器为Kepler-class 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 正确工作两个。轶事证据，当然不是证据。
什么数学条件可以先验判断除法变换成上面的代码序列是否安全？答案可以假设所有的浮点运算都进行在 "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;
}
这个问题要求找到一种方法来识别常量 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 中最小的正值，使得 V 中存在 -M。
- 如果 U 中存在 succ(M)，则这对值将 succ(M) - M 贡献给结果的正值。
- 如果-succ(M)存在于V中，则这对值将负值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 大小，并使用此信息应用 Brisebarre 等人的变换，确信 x * C2 不会下溢。
明确地说，我建议使用如下算法来决定是否将除法转换为更简单的算法：
Y是可以根据他们的算法使用 Brisebarre 等人的方法转换的值之一吗？- 他们方法中的C1和C2是同号吗，还是可以排除被除数是无穷大的可能？
- 他们的方法中的 C1 和 C2 是否具有相同的符号，或者
x可以仅采用 0 的两种表示形式之一？如果在C1和C2符号不同，x只能是0的一种表示的情况下，记得fiddle(**)加上FMA计算的符号，使其产生当x为零时更正零。 - 能否保证分红的幅度足够大，排除
x * C2下溢的可能？
如果四个问题的答案都是“是”，那么除法可以在编译函数的上下文中转化为乘法和FMA。上述静态分析用于回答问题 2.、3. 和 4.
(**) “摆弄符号”意味着在需要时使用 -FMA(-C1, x, (-C2)*x) 代替 FMA(C1, x, C2*x)当 x 只能是两个带符号的零之一时，使结果正确出来
浮点除法结果为：
- 标志旗
- 有效数
- 指数
- 一组标志（上溢、下溢、不精确等 - 参见
fenv()）
仅正确地获得前 3 个部分（但标志集不正确）是不够的。在没有进一步知识的情况下（例如，结果的哪些部分实际上很重要，股息的可能值等）我会假设用常数乘法代替常数除法（and/or 一个令人费解的 FMA 混乱）几乎从来都不安全。
此外；对于现代 CPUs，我也不认为用 2 个 FMA 替换一个分区总是一种改进。例如，如果瓶颈是指令 fetch/decode，那么这个 "optimisation" 会使性能变差。再例如，如果后续指令不依赖于结果（CPU 可以在等待结果的同时并行执行许多其他指令），FMA 版本可能会引入多个依赖性停顿并使性能变差。对于第三个示例，如果所有寄存器都被使用，则 FMA 版本（需要额外的 "live" 变量）可能会增加 "spilling" 并使性能变差。
请注意（在许多但不是所有情况下）除以或乘以 2 的常数倍数可以单独通过加法来完成（具体来说，将移位计数添加到指数）。
我喜欢 @Pascal 的回答，但在优化中，拥有简单且易于理解的转换子集通常比拥有完美的解决方案更好。
所有当前和常见的历史浮点格式都有一个共同点：二进制尾数。
因此，所有分数都是以下形式的有理数：
x / 2n
这与程序中的常量（以及所有可能的以 10 为底的分数）形成对比，这些常量是以下形式的有理数：
x / (2n * 5m)
因此，一种优化将简单地测试 m == 0 的输入和倒数，因为这些数字完全以 FP 格式表示，并且对它们的操作应该产生以下数字在格式内是准确的。
因此，例如，在 .01 到 0.99 的（十进制 2 位数）范围内除以或乘以以下数字将得到优化：
.25 .50 .75
而其他一切都不会。 （我想，先测试一下，哈哈。）
让我第三次重启。我们正在努力加速
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 算法通过预先计算
C1 = 1.0f / y
C2 = -y
这样可以使用
来近似除法 t1 = x * C1
t2 = fmaf(C1, t1, x)
q = fmaf(C2, t2, t1)
朴素的方法适用于 2 y 的所有幂，但除此之外就很糟糕了。例如，对于除数 7、14、15、28 和 30，超过一半的可能 x.
Brisebarre-Muller-Raina 方法对于几乎所有 2 的非幂 y 都同样失败，但更少的 x 产生不正确的结果（不到所有可能 y 的一半 x，因 y) 而异。
Brisebarre-Muller-Raina 文章表明，朴素方法的最大误差为 ±1.5 ULP。
Markstein 方法对 2 的幂 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) 是案例数的一个很好的一阶近似值（在所有可能的情况下float)，其中 Markstein 方法对偶数、非二次方除数 y 产生错误结果。 （或者，16777216/rev(x) 为上限。）
已添加 2016-02-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）的除数，如果除数是这样的，当除以商小于 16777216 的最小二乘方，商为奇数。