将FMA指令用于FFT算法

Question

我有一小段C++代码，随着时间的推移，它已经变成了一个有点有用的快速傅立叶变换库，并且使用SSE和AVX指令使它运行得相当快。诚然，这一切都只是基于基数-2算法，但它仍然适用。我最新的渴望是让蝶形计算与FMA指令一起工作。基本的基2蝶形运算由4个乘法和6个加法或减法组成。一种简单的方法是用2个FMA指令替换2个加法和减法以及2个乘法，从而得到一个数学上相同的蝶形，但显然有更好的方法来做到这一点：

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蝶形的预计算旋转因子将是一件相对简单的事情，当然，我们也可以避免在蝶形开始时的除法操作。

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)的表达式替换成u1

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

这可以用2个FMA来完成。

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

(注意，这并不意味着这个变体会自动运行得更快，因为它有不同的数据依赖链)