多项式乘法与离散傅里叶变换 (3)

$\begin{align*} &A^{[0]}(\omega _4^0) =B^{[0]}(\omega _4^0)+\omega _4^0B^{[1]}(\omega _4^0) =B^{[0]}(\omega _2^0)+\omega _4^0B^{[1]}(\omega _2^0) \\ &A^{[0]}(\omega _4^{-1})=B^{[0]}(\omega _4^{-2})+\omega _4^{-1}B^{[1]}(\omega _4^{-2})=B^{[0]}(\omega _2^{-1})+\omega _4^{-1}B^{[1]}(\omega _2^{-1}) \\ &A^{[0]}(\omega _4^{-2})=B^{[0]}(\omega _4^{-4})+\omega _4^{-2}B^{[1]}(\omega _4^{-4})=B^{[0]}(\omega _2^0)-\omega _4^0B^{[1]}(\omega _2^0) \\ &A^{[0]}(\omega _4^{-3})=B^{[0]}(\omega _4^{-6})+\omega _4^{-3}B^{[1]}(\omega _4^{-6})=B^{[0]}(\omega _2^{-1})-\omega _4^{-1}B^{[1]}(\omega _2^{-1}) \\ \end{align*}$

$\begin{align*} B^{[0]}(x)&=f(\omega _8^0)+f(\omega _8^4)x \\ &=f(\omega _8^0)+xf(\omega _8^4) \end{align*}$

$\begin{align*} &B^{[0]}(\omega _2^0)=a_0+\omega _2^0 a_4 \\ &B^{[0]}(\omega _2^{-1})=a_0+\omega _2^{-1} a_4=a_0-\omega _2^0 a_4 \\ \end{align*}$

相信给出一张图片，你就能明白啦！！！

注意：处理的长度只能是2的次幂，所以对于多项式相乘，应该将长度n变成：大于等于n的最小的2的次幂，再乘以2。以保证不会溢出。 记得除以n。 代码如下： void FFT(complex<double> *a, int n, int inv) //inv=1 FFT inv=-1 IFFT { if (n == 1) return; int mid = n / 2; for (int i = 0; i < mid; i++) { rev[i] = a[i*2]; rev[i + mid] = a[i * 2 + 1]; } for (int i = 0; i < n; i++) { a[i] = rev[i]; } FFT(a, mid, inv); FFT(a + mid, mid, inv); for (int i = 0; i < mid; i++) { complex<double> w(cos(2 * pi * i / n), inv * sin(2 * pi * i / n)); complex<double> l = a[i]; complex<double> r = w * a[i+mid]; a[i] = l+r; a[i + mid] = l-r; } } 优化

上面的代码层层递归，还是很耗时的，如果我们一开始就将各位数放到他们应该去的位置，就可以快很多了。

$$\begin{matrix} a_0 & a_1 & a_2 & a_3 & a_4 & a_5 & a_6 & a_7 \\ \downarrow & \downarrow & \downarrow & \downarrow & \downarrow & \downarrow & \downarrow &\downarrow \\ 000 & 001 & 010 & 011 & 100 & 101 & 110 & 111\\ \downarrow & \downarrow & \downarrow & \downarrow & \downarrow & \downarrow & \downarrow &\downarrow \\ 位反转 & 位反转 & 位反转 & 位反转 & 位反转 & 位反转 & 位反转 & 位反转 & \\ \downarrow & \downarrow & \downarrow & \downarrow & \downarrow & \downarrow & \downarrow &\downarrow \\ 000 & 100 & 010 & 110 & 001 & 101 & 011 & 111 \\ \downarrow & \downarrow & \downarrow & \downarrow & \downarrow & \downarrow & \downarrow &\downarrow \\ a_0 & a_4 & a_2 & a_6 & a_1 & a_5 & a_3 & a_7 \end{matrix}$$

代码如下：

void FFT(complex<double>* a, int n, int inv) { int bit = 0; while ((1 << bit) < n) bit++; for (int i = 0; i < n; i++) { rev[i] = a[(rev32bit(i) >> (32 - (bit)))]; } for (int i = 0; i < n; i++) { a[i] = rev[i]; } for (int mid = 1; mid < n; mid *= 2) { complex<double> w(cos(pi / mid), inv * sin(pi / mid)); for (int i = 0; i < n; i += mid * 2) { complex<double> temp = 1; for (int j = 0; j < mid; j++, temp *= w) { //complex<double> w(cos(pi * j / mid), inv * sin(pi * j / mid)); complex<double> l = a[i + j]; complex<double> r = temp * a[i + j + mid]; a[i + j] = l + r; a[i + j + mid] = l - r; } } } }

32位反转算法代码如下：

unsigned int rev32bit(unsigned int x) { x = ((x & 0x55555555) << 1) | ((x & 0xaaaaaaaa) >> 1); x = ((x & 0x33333333) << 2) | ((x & 0xcccccccc) >> 2); x = ((x & 0x0f0f0f0f) << 4) | ((x & 0xf0f0f0f0) >> 4); x = ((x & 0x00ff00ff) << 8) | ((x & 0xff00ff00) >> 8); x = ((x & 0x0000ffff) << 16) | ((x & 0xffff0000) >> 16); return x; } 注意：上面的代码不能直接用，例如4的位反转为0x20000000。即00000000 00000000 00000000 00000100b反转后为00100000 00000000 00000000 00000000b，我们需要右位移29（32-$\log_2^{长度}$)位，得到0x0001，即为1。 扩展 求解第一个大于等于n的2的次幂