使用常数整数除法器进行高效的浮点除法

最近的一个问题 ,是不是允许编译器用浮点乘法来取代浮点运算,这就激励我提出这个问题。

在严格的要求下,代码转换后的结果应该与实际的除法运算按位相同,但对于二进制IEEE-754algorithm来说,对于二的幂因子是可能的。 只要除数的倒数是可以代表的,乘以除数的倒数即可得出与该除数相同的结果。 例如,乘以0.5可以代替2.0

然后人们想知道这样的replace是如何工作的,假设我们允许任何简短的指令序列来代替除法,但运行速度明显更快,同时提供相同的结果。 除了普通的乘法之外,还特别允许融合的乘加运算。 我在评论中指出了以下相关文章:

Nicolas Brisebarre,Jean-Michel Muller和Saurabh Kumar Raina。 提前知道除数时加速正确舍入的浮点除法。 IEEE Transactions on Computers,Vol。 53,第8号,2004年8月,第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必须满足的各种条件才能使该algorithm正常工作。 正如人们容易观察到的那样,当头部和尾部的符号不同时,该algorithm存在无穷大和零的问题。 更重要的是,由于商尾zl * x的计算遭受下溢,所以它将不能提供非常小的股息x的正确结果。

本文还提供了另一种基于FMA的分割algorithm,由Peter Markstein在IBM任职时开创。 相关的参考资料是:

PW Markstein。 计算IBM RISC System / 6000处理器上的基本function。 IBM研究与发展杂志,Vol。 34,No.1,1990年1月,第111-119页

在Marksteinalgorithm中,首先计算一个倒数rc ,由此形成初始商q = x * rc 。 然后,用FMA精确计算除法的余数为r = fma(-y,q,x) ,最终计算出一个改进的更精确的商,即q = fma(r,rc,q)

这个algorithm对x也是有问题的,它们是零或者无穷(可以很容易的用适当的条件执行来解决),但是使用IEEE-754单精度float数据的详尽的testing表明它为所有可能的分红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); } 

在大多数处理器体系结构中,这应该转换成无分支指令序列,使用预测,条件移动或selecttypes指令。 举一个具体的例子:对于除以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 

对于我的实验,我写了如下所示的微小的Ctesting程序,按递增顺序通过整数除数,并且对每个程序都进行了彻底的testing,以对照正确的划分对上述代码序列进行testing。 它打印一个通过这个穷举testing的除数列表。 部分输出如下所示:

 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, 

为了将replacealgorithm结合到编译器中作为优化,可以安全地应用上述代码变换的除数的白名单是不切实际的。 到目前为止,程序的输出(以每分钟大约一个结果的速率)表明,对于那些奇数整数或二的幂的因子,快速代码在x所有可能的编码中正确地工作。 轶事证据,当然不是certificate。

什么样的math条件可以事先确定划分到上述代码序列的转换是否安全? 答案可以假定所有的浮点运算都是以“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; } 

让我第三次重新启动。 我们正在努力加速

  q = x / y 

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

天真的algorithm是通过预先计算的倒数,

  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) 

Marksteinalgorithm将天真的方法与两个融合的乘法相结合,如果初始方法在最不重要的地方产生一个单位内的结果,则通过预先计算得出正确的结果

  C1 = 1.0f / y C2 = -y 

这样可以近似使用divison

  t1 = x * C1 t2 = fmaf(C1, t1, x) q = fmaf(C2, t2, t1) 

天真的方法适用于所有的两个权力,但否则它是非常糟糕的。 例如,对于除数7,14,15,28和30,对于所有可能x一半以上,它会产生不正确的结果。

Brisebarre-Muller-Raina方法对于两个y几乎所有非功率也类似地失败,但是less得多的x产生不正确的结果(小于所有可能的x的百分之五,取决于y )。

Brisebarre-Muller-Raina文章显示,天真方法中的最大误差是±1.5 ULP。

Markstein方法对于两个y幂,对于奇数y也产生正确的结果。 (我还没有findMarkstein方法的失败奇数整数除数。)


对于Markstein方法,我已经分析了除数1 – 19700( 原始数据在这里 )。

绘制失败案例的数量(水平轴上的除数,对于所述除数,Markstein方法失效的x的数值),我们可以看到一个简单的模式:

马克斯坦失败案例answers/markstein.png

请注意,这些图表的水平轴和垂直轴都是对数的。 奇数除数是没有点的,因为这个方法对于我testing的所有奇数除数都能得到正确的结果。

如果我们把x轴改成反转位(反转的二进制数字,即0b11101101→0b10110111, 数据 )的除数,我们有一个非常清晰的模式: Markstein故障情况,位反向除数http://www.nominal- animal.net/answers/markstein-failures.png

如果我们通过点集的中心绘制一条直线,我们得到曲线4194304/x 。 (记住,情节只考虑一半可能的花车,所以当考虑所有可能的花车时,将其加倍。) 8388608/x2097152/x完全包含整个错误模式。

因此,如果我们使用rev(y)计算除数y的位反转,那么8388608/rev(y)是Markstein方法产生不正确情况的情况数(在所有可能的float中)的一个好的一阶近似值导致一个平均的,非二次除数的除数y 。 (或者,上限为16777216/rev(x) 。)

增加了2016-02-28:给出了任何整数(二进制32)除数,我发现使用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 

在我testing的Markstein故障情况下(但我还没有充分testing大于8388608的除数),这产生了一个正确的误差估计,误差在±1以内。 最后的划分应该是这样的,它报告没有虚假的零,但我不能保证(还)。 它没有考虑到非常大的因子(比如说0x1p100,或者1e + 30,以及更大的数量级),这些因素都有下溢的问题 – 无论如何,我肯定会排除这样的因数。

在初步testing中,这个估计似乎是非常准确的。 我没有画一个比较估计值和因子1到20000的实际误差的情节,因为所有的点都完全相同。 (在这个范围内,估计是准确的,或者是一个太大。)本质上,估计重现了这个答案中的第一个图。


Markstein方法的失败模式是规则的,非常有趣。 该方法适用于两个因子的所有权力,以及所有的奇数整除因子。

对于大于16777216的除数,我始终如一地看到与除数除以2的最小次方产生小于16777216的值相同的错误。例如,0x1.3cdfa4p + 23和0x1.3cdfa4p + 41,0×1。 d8874p + 23和0x1.d8874p + 32,0×1.cf84f8p + 23和0x1.cf84f8p + 34,0×1.e4a7fp + 23和0x1.e4a7fp + 37。 (在每一对中,尾数是一样的,只有两个的幂是变化的。)

假设我的testing平台没有错误,这意味着Markstein方法也可以使用大于16777216的除数(数值小于1e + 30),如果除数除以2的最小次幂产量小于16777216的商,商数很奇怪。

这个问题要求提供一种方法来确定常数Y的值,这使得将x / Y转换为更便宜的使用FMA计算所有可能的x值的计算是安全的。 另一种方法是使用静态分析来确定x可以采取的值的过度近似,从而可以应用通常不稳定的转换,即知道转换后的代码与原始分割不同的值不会发生。

使用适合于浮点计算问题的浮点值集合的表示,即使是从函数开始的前向分析,也可以产生有用的信息。 例如:

 float f(float z) { float x = 1.0f + z; float r = x / Y; return r; } 

假设默认的四舍五入模式(*),在上面的函数中, x只能是NaN(如果input是NaN),+ 0.0f,或者一个大于2 -24的数字,而不是-0.0f或比2-24更接近零的任何东西。 这certificate了对于常数Y许多值,转换成问题所示的两种forms之一。

(*)的假设,没有它,许多优化是不可能的,并且C编译器已经创build,除非程序明确使用#pragma STDC FENV_ACCESS ON


转发静态分析,预测上面的x的信息可以基于一组expression式可以作为一个元组的浮点值的表示:

  • 一组可能的NaN值的表示(由于NaN的行为是低标准的,所以select只使用一个布尔值, 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。
    • 如果-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分解转换成另一种forms,因为6不是可以无条件地转换分解的价值之一。 相反,我所build议的是从函数的开始开始应用一个简单的值分析,在这种情况下,确定x+0.0f或至less2 -24的有限浮点数,并且使用这个信息来应用Brisebarre等人的转变,对x * C2不下溢的知识充满信心。

为了明确起见,我build议使用如下的algorithm来决定是否将分割转换为简单的分割:

  1. Y是根据他们的algorithm使用Brisebarre等人的方法可以转化的价值之一吗?
  2. 从他们的方法中,C1和C2有相同的符号,还是有可能排除红利是无限的可能性?
  3. 他们的方法C1和C2是否有相同的符号,或者x只能取0的两个表示之一? 如果在C1和C2有不同符号的情况下, x只能是0的一个表示forms,请记住用FMA计算的符号来摆弄(**),使x在零时产生正确的零点。
  4. 能否保证股利的大小足以排除x * C2下溢的可能性?

如果对这四个问题的答案是“是”,那么在编译函数的情况下,可以将该部分转换成乘法和FMA。 上述静态分析用于回答问题2,3和4。

(**)“摆弄符号”是指在必要时使用-FMA(-C1,x,(-C2)* x)代替FMA(C1,x,C2 * x)正确的时候,x只能是两个有符号的零之一

浮点除法的结果是:

  • 一个标志标志
  • 一个有效数字
  • 一个指数
  • 一组标志(溢出,下溢,不精确等 – 见fenv()

获得前3个部分是正确的(但标志集不正确)是不够的。 没有进一步的知识(例如结果的哪一部分哪一部分实际上是重要的,股息的可能值等),我会假设用常数(和/或复杂的FMA混乱)乘以常数代替除法是几乎从来不安全。

此外; 对于现代的CPU,我也不会认为用2个FMAreplace一个分区总是一个改进。 例如,如果瓶颈是取指令/解码,那么这个“优化”会使性能变差。 又如,如果后续指令不依赖于结果(在等待结果的同时CPU可以并行执行许多其他指令),则FMA版本可能会引入多个依赖性延迟,从而使性能变差。 对于第三个例子,如果所有寄存器都被使用,那么FMA版本(需要额外的“实时”variables)可能会增加“溢出”并使性能变差。

注意(在许多情况下,但不是全部情况下)除以2的常数倍数的分割或乘法可以通过单独添加(具体地,向指数添加移位计数)来完成。

我喜欢@Pascal的答案,但是在优化中,最好是有一个简单的,被很好理解的转换子集,而不是一个完美的解决scheme。

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

因此,所有的分数都是有理数的forms:

x / 2 n

这与程序中的常量(以及所有可能的基数为10的分数)形成鲜明对比,这些常数是以下forms的有理数:

x /(2 n * 5 m

因此,一个优化可以简单地testinginput和m == 0的倒数,因为这些数字完全以FP格式表示,并且与它们的操作应该产生在格式内精确的数字。

因此,例如,在.010.99的(十进制2位)范围内,除以下数字将被优化:

 .25 .50 .75 

而其他一切都不会。 (我想,先试一下,哈哈。)