使用FMA指令进行FFT算法

Question

我有一些C++代码，随着时间的推移，它已成为一个相当有用的FFT库，并且已经使用SSE和AVX指令使其运行得相当快。当然，这全部只是基于radix-2算法，但仍然保持不错的效果。我的最新想法是让蝴蝶计算使用FMA指令。基本的radix-2蝴蝶包括4个乘法和6个加减法。简单的方法是将其中2个加减法和2个乘法替换为2个FMA指令，从而得到数学上相同的蝴蝶，但显然有更好的方法来实现这一点：

https://books.google.com/books?id=2HG0DwAAQBAJ&pg=PA56&lpg=PA56&dq=radix+2+fft+fma&source=bl&ots=R5XDWyYBVv&sig=ACfU3U0S2n1hcgiP63LTKMxI5Oc85eEZaQ&hl=en&sa=X&ved=2ahUKEwiz_I3PsrToAhVoHzQIHYmVDGIQ6AEwDXoECAoQAQ#v=onepage&q=radix%202%20fft%20fma&f=false

ci1 = ci1 / cr1
u0 = zinr(0)
v0 = zini(0)
r = zinr(1)
s = sini(1)
u1 = r - s * ci1
v1 = r * ci1 + s
zoutr(0) = u0 + u1 * cr1
zouti(0) = v0 + v1 * cr1
zoutr(1) = u0 - u1 * cr1
zouti(1) = v0 - v1 * cr1

作者用6个FMA替换了所有的10个加、减和乘法，前提是扭转因子的虚部被实部除尽。其中一部分文字写道“注意cr1 != 0”。这基本上就是我的问题所在。对于所有的扭转因子，数学似乎都像广告中宣传的那样工作得很好，除了当实际扭转为零时，我们最终会除以零。在这里效率非常关键，当cr1 == 0时将代码分支到不同的蝴蝶形式并不是一个好的选择，特别是当我们使用SIMD同时处理多个扭转和蝴蝶形式时，也许只有cr1的一个元素等于0。我的直觉告诉我，当cr1 == 0时，cr1和ci1应该是完全不同的值，而FMA代码仍然会得出正确的答案，但我似乎无法弄清楚这一点。如果我能搞清楚，修改预计算的FMA蝴蝶形式的扭转因子将是一个相对简单的事情，当然，我们还可以避免在蝴蝶形式开始时进行除法运算。

- Kumputer

---

我无法通过您提供的链接阅读任何内容。部分相关：如果您对效率感兴趣，是否尝试使用基数-4？ - Damien

Answer 1

这本书似乎暗示着cr1 != 0总是成立的，但不幸的是，情况并非总是如此（当旋转角度为PI/2时）。

我认为你不能通过调整扭曲因子来解决这个问题。我唯一看到的选择是使用一些非常小的数字代替零。这可能有效，但很丑陋，并且在某些情况下可能会导致不准确性。

可能的解决方案：

• 将循环拆分为两个，并特别处理这个中心情况（发生除以零的情况）
• 不是除以cr1，而是除以ci1，并相应地修改公式。这种情况仍然存在除以零的情况，但它将发生在循环的第一次迭代中。因此，您只需要特别处理第一次迭代（所以只需要一个循环）。
• 使用不同的FMA公式：

请注意：

zoutr(1) = u0 - u1 
         = u0 - u1 - (u0 + u1) + (u0 + u1) 
         = u0 - u1 - zoutr(0) + u0 + u1 
         = 2*u0 - zoutr(0)

所以，这个操作可以在1个FMA中完成。

如果你将代入zoutr(0)的表达式中：

zoutr(0) = u0 + u1
         = u0 + r*cr1 - s*ci1

这可以用两个FMA完成。

计算zouti的方法与zoutr相同。因此，您需要使用6个FMA操作，这是与书中相同数量的操作。

（请注意，这并不意味着此变体将自动运行得更快，因为它具有不同的数据依赖关系链）

- geza