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

4 回答

• 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 的百分之几，取决于 y ） .

  Brisebarre-Muller-Raina文章显示，天真方法的最大误差为±1.5 ULPs .

  Markstein方法得到两个 y 的幂的正确结果，也对于奇数整数 y . （我没有找到Markstein方法的失败奇数整数除数 . ）

  ---

  对于Markstein方法，我分析了1 - 19700（raw data here）的除数 .

  绘制失败案例的数量（水平轴上的除数，Markstein逼近所述除数的 x 的值的数量），我们可以看到一个简单的模式：

  Markstein failure cases http://www.nominal-animal.net/answers/markstein.png

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

  如果我们将x轴更改为反转位（反向的二进制数字，即0b11101101→0b10110111，data）除数，我们有一个非常清晰的模式：Markstein failure cases, bit reverse divisor http://www.nominal-animal.net/answers/markstein-failures.png

  如果我们在点集的中心绘制一条直线，我们得到曲线 4194304/x . （请记住，该图只考虑了一半可能的浮点数，因此在考虑所有可能的浮点数时，加倍它 . ） 8388608/x 和 2097152/x 完全包含整个错误模式 .

  因此，如果我们使用 rev(y) 计算除数 y 的位反转，那么 8388608/rev(y) 是一个良好的一阶近似的情况数（在所有可能的浮点数中），其中Markstein方法产生一个偶数，非幂的不正确结果两个除数 y . （或者， 16777216/rev(x) 为上限 . ）

  添加2016-02-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失效案例中，这会产生一个正确的误差估计值±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.（在每对中，尾数相同，只有2的幂变化 . ）

  假设我的测试平台没有错误，这意味着Markstein方法在大小上也会使用大于16777216的除数（但是比1e 30小），如果除数是这样的，除以得到的最小二次幂商数小于16777216，商为奇数 .

  回复于 2024-04-27T01:00:56+08:00
• 7

  这个问题要求一种方法来识别常量 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.

  • 如果su中存在succ（M），那么这对值将succ（M）-M贡献给结果的正值 .

  • 如果-succ（M）存在于V中，则该对值将负值M-succ（M）贡献给结果的负值 .

  • 如果在U中存在pred（M），那么这对值将负值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 只能是零的一个表示的情况下，记住用基于FMA的计算的符号来调整（**）以使其在 x 为零时产生正确的零 .

  • 可以保证股息的大小足以排除 x * C2 下溢的可能性吗？

  如果对四个问题的答案为“是”，则可以在编译的函数的上下文中将除法转换为乘法和FMA . 上述静态分析用于回答问题2,3和4 .

  （**）“摆弄标志”是指使用-FMA（-C1，x，（ - C2）* x）代替FMA（C1，x，C2 * x），这样才能使结果出来当x只能是两个带符号的零之一时正确

  回复于 2024-04-27T01:00:56+08:00
• 6

  浮点除法的结果是：

  • 标志旗

  • 有意义的

  • 一个指数

  • 一组标志（溢出，下溢，不精确等等 - 见 fenv() ）

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

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

  注意（在许多但不是所有情况下）除以2的常数倍可以单独加法（具体地说，向指数添加移位计数） .

  回复于 2024-04-27T01:00:56+08:00
• 1

  我喜欢@Pascal 's answer but in optimization it'通常更好地拥有一个简单且易于理解的变换子集，而不是一个完美的解决方案 .

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

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

  x / 2n

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

  x /（2n * 5m）

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

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

  .25 .50 .75
  

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

  回复于 2024-04-27T01:00:56+08:00