W(nk,N)=e^(-j2pi/Nnk)//蝶形因子
W(nk,N)=W(-nk,N)=W((N-n)k,N)=W(n(N-k),N) 相当于W(-nk,N)+e(-j*2npi)或+e(-j2kpi)对称性
W(nk,N)=W((N+n)k,N)=W(n(N+k),N)周期性
W(nk,N)=W(nmk,mN)=W(n/mk,N/m) = e(-j2pi/N*nk)=e(-j2mpi/Nmnk) 可约性
W(0,N)=1 = e^0 = 1 W(N/2,N)=-1=e-jpi W(k+N/2,N)=-W(k,N) 特殊点
X(k)=Y(k)+W(k,N)Z(k)//因为E(W((2n+1)k,N))=W(k,N)W(20+1)k,N)+W(k,N)W(21+1)k,N)+...=W(k,N)E(W(2nk,N)=W(k,N)E(W(nk,N/2))
X(k+2/N)=Y(k)-W(k,N)Z(k)//这里Y(k)和Z(k)都有N/2的周期又W(k+2/N,N)=-W(k,N)
一直进行奇偶分解到2点DFT，一共log(2,N)级，每级蝶形运算N/2，每次蝶形运算复数乘法1次，复数加法2次
比直接DFT运算量NN降低至N/2*log(2,N)
然后就是进行位交换了，比如第3个为011，位逆序为110交换到第6个处
添加一个bin神链接：http://www.cnblogs.com/kuangbin/archive/2013/07/24/3210389.html
其实这个3层还是没太懂
#include <bits/stdc++.h>
using namespace std;
const int N = 200010;
int a[N];
map<int,int> p;
const double pi = acos(-1.0);
class complex {
public:
double r, i;
complex(){}
complex(double r, double i) {
this->r = r;
this->i = i;
}
complex operator +(const complex &rhs) {
return complex(r + rhs.r, i + rhs.i);
}
complex operator -(const complex &rhs) {
return complex(r - rhs.r, i - rhs.i);
}
complex operator *(const complex &rhs) {
return complex(r*rhs.r - i*rhs.i, r*rhs.i + i*rhs.r);
}
};
//改变输入序列，通过2进制逆序值
void change(complex y[], int len) {
int x = 0;
while ((1 << x) < len)x++;
for (int i = 0;i < len;i++) {
int j = i;
int k = 0;
while (k < x/2) {
if ((j&(1 << k))&&!(j&(1<<(x-k-1)))){
j = j - (1 << k) + (1 << (x - k - 1));
}
else if (!(j&(1 << k)) && (j&(1 << (x - k - 1)))) {
j = j + (1 << k) - (1 << (x - k - 1));
}
k++;
}
if (j > i)
swap(y[i], y[j]);
}
}
void fft(complex y[], int len, int on)
{
change(y, len);
/*complex Y[1000];
for (int i = 0;i < len;i++)
Y[i].r = y[i].r;*/
/*for (int k = 0;k < len;k++)
{
complex sc=complex(0, 0);
for (int j = 0;j < len;j++) {
complex c = complex(Y[j].r * cos(on*2 * pi*k*j / len), -Y[j].r * sin(on*2 * pi*k*j / len));
sc = sc+c;
}
y[k] = sc;
}*/
for (int l = 2;l <= len;l <<= 1) {
complex wn = complex(cos(-on * 2 * pi / l), sin(-on * 2 * pi / l));
for (int j = 0;j < len;j += l) {
complex wnp = complex(1, 0);
for (int k = j;k < j + l / 2;k++) {
complex u = y[k];
complex t = wnp*y[k + l / 2];
y[k] = u + t;
y[k + l / 2] = u - t;
wnp = wnp*wn;
}
}
}
//for (int h = 2; h <= len; h <<= 1)
//{
/*complex wn(cos(-on * 2 * pi / len), sin(-on * 2 * pi / len));
for (int j = 0;j < len;j ++)
{
complex w(1, 0);
for (int k = 0;k < len;k++)
{
complex u = y[k];
complex t = w*y[k + len / 2];
y[k] = u + t;
y[k + len / 2] = u - t;
w = w*wn;
}
}*/
//}
if (on == -1)
for (int i = 0;i < len;i++)
y[i].r /= len;
}
char s1[N], s2[N];
complex a1[N], a2[N];
int ans[N],tmp[N];
int main()
{
while (~scanf("%s%s", s1, s2)) {
int len1 = strlen(s1);
int len2 = strlen(s2);
int tlen = len1 + len2;
int len = 1;
while (len < tlen)len *= 2;//长度>=len1+len2-1的圆周卷积才等于线性卷积，又用的是基2-fft所以长度为2^x
for (int i = 0;i < len;i++) {
if (i < len1)a1[i].r = s1[len1-1-i] - '0';
else a1[i].r = 0;
a1[i].i = 0;
}//补齐0
for (int i = 0;i < len;i++) {
if (i < len2)a2[i].r = s2[len2-1-i] - '0';
else a2[i].r = 0;
a2[i].i = 0;
}//补齐0
//由于输出序列为顺序的所以采取抽取时间的fft方法
fft(a1, len, 1);//fft
fft(a2, len, 1);//fft
for (int i = 0;i < len;i++)
a1[i] = a1[i]*a2[i];
//Y(k)=H(k)*X(k)
//ifft可以用fft子程序实现，先求一次共轭，再fft，序列集体除以len，再求一次共轭，由于虚数部分为0，基本不用求共轭
fft(a1, len, -1);//ifft
for (int i = 0;i < len;i++) {
ans[i] = (int)(a1[i].r + 0.5);
}
for (int i = 1;i <= len;i++) {
ans[i] += ans[i - 1] / 10;
ans[i - 1] %= 10;
}
int x = 0;
len = len1 + len2 - 1;
for (;len>0;len--)
if (ans[len] != 0)break;
for (;len>=0;len--)
printf("%d", ans[len]);
printf("
");
}
return 0;
}