# bzoj 3771: Triple 快速傅里叶变换 FFT

#### 题解:

$A(x) = x + x^2 + x^3 + ...$

$A(x)*A(x)$

$B(x) = x^2 + x^4 + x^6 + ...$

$C(x) = x^3 + x^6 + x^9 + ...$

$\frac{A(x)}{1!} + \frac{A(x)^2-B(x)}{2!} + \frac{A(x)^3 - 3A(x)B(x) + 2C(x)}{3!}$

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long ll;
x=0;char ch;bool flag = false;
while(ch=getchar(),ch<'!');if(ch == '-') ch=getchar(),flag = true;
while(x=10*x+ch-'0',ch=getchar(),ch>'!');if(flag) x=-x;
}
inline int cat_max(const int &a,const int &b){return a>b ? a:b;}
inline int cat_min(const int &a,const int &b){return a<b ? a:b;}
const int maxn = 208000;
const double pi = acos(-1);
struct complex{
double x,y;
complex(){}
complex(double a,double b){x=a;y=b;}
complex operator + (const complex &r){return complex(x+r.x,y+r.y);}
complex operator - (const complex &r){return complex(x-r.x,y-r.y);}
complex operator * (const complex &r){return complex(x*r.x-y*r.y,x*r.y+y*r.x);}
complex operator / (const double &r){return complex(x/r,y/r);}
};
void FFT(complex *x,int n,int p){
for(int i=0,t=0;i<n;++i){
if(i > t) swap(x[i],x[t]);
for(int j=n>>1;(t^=j) < j;j>>=1);
}
for(int m=2;m<=n;m<<=1){
complex wn(cos(p*2*pi/m),sin(p*2*pi/m));
for(int i=0;i<n;i+=m){
complex w(1,0),u;
int k = m>>1;
for(int j=0;j<k;++j,w=w*wn){
u = x[i+j+k]*w;
x[i+j+k] = x[i+j] - u;
x[i+j] = x[i+j] + u;
}
}
}
if(p == -1) for(int i=0;i<n;++i) x[i] = x[i]/n;
}
complex a[maxn],b[maxn],c[maxn],d[maxn];
int main(){
int nn = 0;
for(int i=0,x;i<n;++i){
a[x].x += 1.0;
b[x<<1].x += 1.0;
c[x*3].x += 1.0;
nn = cat_max(nn,x*3);
}
int len = 1;
for(int i=1;(i>>1)<nn;i<<=1) len = i;
FFT(a,len,1);FFT(b,len,1);
for(int i=0;i<len;++i) d[i] = a[i]*a[i]*a[i]/6.0 + (a[i]*a[i] - a[i]*b[i] - b[i])/2.0 + a[i];
FFT(d,len,-1);
for(int i=0;i<len;++i) d[i] = d[i] + c[i]/3.0;
for(int i=0;i<len;++i){
int x = (int)(d[i].x + 0.5);
if(x == 0) continue;
printf("%d %d\n",i,x);
}
getchar();getchar();
return 0;
}



posted @ 2017-01-29 09:38  Sky_miner  阅读(380)  评论(0编辑  收藏  举报