快乐暑假（十）——关于FFT算法中倒置二进制算法的一点补充

本文链接：https://blog.csdn.net/sinat_40872274/article/details/99122715

快速傅里叶变换经过了使用雷德算法进行迭代优化，在效率上已经有了很大的进步。
但是，如果一位一位翻转，尽管位数不是很多，数据个数非常多的时候还是比较慢。所以，我们用一个类似DP的方法来实现这个功能。

~~目前还不知道这个算法的名字~~


 void get_rev(int bit)//bit表示二进制的位数
{  
    for(int i=0;i<(1<<bit);i++)对1~2^bit-1中的所有数做长度为bit的二进制翻转
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
}

~~这个代码一看除了知道使用很多位运算完全不知道在干什么。~~
我们可以把一个二进制看成两部分：

它的前bit-1是一部分
最后一位是一部分

在雷德算法中提到过，二进制翻转就是把它的最后一位当成首位。我们仔细考虑这个问题，其实可以发现这个问题其实有一个最优子结构。

每一个二进制翻转就相当于把他的最后一位当成首位，然后在前面接上它的前bit-1位的进制翻转。

下面的问题就转换为，如何找到一个数的前bit-1位的进制转换是多少。
其实也很好找，“i”的前bit-1位的数值其实就是 ${\lfloor\frac{i}{2}\rfloor}$ （向下取整）的值，也就是i右移一位的的值。

举个例子：
5: 101
它的前2位为：10
5右移一位就是：10
在这里插入图片描述

但是，实际上，i>>1的翻转与“i”的前bit-1位的翻转是有一些不同的，因为我们的二进制翻转始终以bit位为标准，所以i>>1会比“i”的前bit-1位多出一个前导零，而翻转之后就会多出一个“后缀零”，所以“i”的前bit-1位的翻转要去掉那个“后缀零”，也就是“rev[i>>1]>>1”。

还是以5举例，5>>1 = 2 (010),但5的前两位是10，多了一个前置的‘0’，所以我们翻转之后这个前置的‘0’就到了最后，所以我们再次右移把后边的0去掉。

在这里插入图片描述

我们只要把末尾乘上2^(bit-1)变成首位，再加上rev[i>>1]>>1就是我们要的答案了。

原本的代码应该是这样

rev[i]=(rev[i>>1]>>1)+((i&1)<<(bit-1));

但是代码里用了按位或‘ | ’来加快速度。~~反正都用了那么多位运算，就不在乎这一个两个了~~

因为(i&1)^(bit-1)的后bit-1位都是零，0|1=1 0|0=0，0+1=1 0+0=0。此时“按位或”与加法等价。

完整FFT

#include <cstdio>
#include <cmath>
#include <iostream>
#include <algorithm>
#include <cstring>
using namespace std;

const double PI = acos(-1);
const int SZ = 1e4 + 7;
int rev[SZ];
template <typename T> 
struct complex {
	T r, m;
	complex(T _r=0, T _m=0): r(_r), m(_m) {}
	complex operator + (const complex &b) {
		return complex(r+b.r, m+b.m);
	}
	complex operator - (const complex &b) {
		return complex(r-b.r, m-b.m);
	}
	//(a+ib)*(c+id)=ac-bd + i(ad+bc)
	complex operator * (const complex &b) {
		return complex(r*b.r-m*b.m, r*b.m+m*b.r);
	}
};
complex<double> a[SZ], b[SZ], c[SZ<<1];

void get_rev(complex<double> a[],int n)//bit表示二进制的位数
{
	memset(rev,0,sizeof(rev));
	int bit=0;
    while ((1<<bit)<n) bit++;
    for(int i = 0;i < n; i++){//对1~2^bit-1中的所有数做长度为bit的二进制翻转
		rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
		if (i<rev[i]) std::swap(a[i],a[rev[i]]);//不加这条if会交换两次（就是没交换）
	}
}

//op为1是正变换，-1是逆变换 
void fft(complex<double> a[], int n, int op) {
	get_rev(a,n);
	for (int i = 2; i <= n; i <<= 1) {
		complex<double> wn(cos(2*PI*op/i), sin(2*PI*op/i));
		for (int j = 0; j < n; j += i) {
			complex<double> w(1, 0);
			for (int k = j, end = j+i/2; k < end; ++k) {
				complex<double> t1 = a[k], t2 = w * a[k+i/2];
				a[k] = t1 + t2, a[k+i/2] = t1 - t2;
				w = w * wn;
			}
		}
	}
	if (op == -1) for (int i = 0; i < n; ++i) a[i].r /= n;
}
int main()
{
	int n, m;
	//输入多项式的最高幂，以及系数 
	scanf("%d%d", &n, &m);
	for (int i = 0; i <= n; i++) scanf("%lf", &a[i]);
	for (int i = 0; i <= m; i++) scanf("%lf", &b[i]);
	m += n, n = 1; 
	//寻找合适的2的幂，赋值为n 
	while (n < m) n <<= 1;
	//正变换 
	fft(a, n, 1), fft(b, n, 1);
	//0~n-1,n为2的整数幂 
	for (int i = 0; i < n; i++) c[i] = a[i] * b[i];
	//逆变换 
	fft (c, n, -1);
	for (int i = 0; i <= m; i++) printf("%.0f ", c[i].r);
	return 0;
}
/*
3 2
1 1 2 3
1 2 3
1 3 7 10 12 9
*/