# bzoj3771 Triple

### 题目链接

bzoj3771 Triple

“这把斧头，是不是你的？”

“这把斧头，是不是你的？”

“这把斧头，是不是你的？”

“你看看你现在的样子，真是丑陋！”

fft优化一下

### 代码

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>

int x = 0;
char c = getchar();
while(c < '0' || c > '9') c = getchar();
while(c <= '9' && c >= '0')x = x * 10 + c- '0',c = getchar();
return x ;
}
#define pi acos(-1.0)
struct Complex {
double x,y;
Complex(double a = 0,double b = 0) : x(a),y(b)	{};
} ;
Complex operator + (Complex a,Complex b) { return Complex(a.x + b.x,a.y + b.y);}
Complex operator - (Complex a,Complex b) { return Complex(a.x - b.x,a.y - b.y);}
Complex operator * (Complex a,Complex b) { return Complex(a.x * b.x - a.y * b.y,a.x * b.y + a.y * b.x);}

const int maxn = 1000007;
int n;
Complex A[maxn],B[maxn],C[maxn];
void fft(Complex *a,int n,int type) {
for(int i = 0,j = 0;i < n;++ i) {
if(i < j) std::swap(a[i],a[j]);
for(int k = n >> 1;(j ^= k) < k;k >>= 1);
}
for(int m = 2;m <= n;m <<= 1) {
Complex w1 = Complex (cos(2 * pi / m),type * sin(2 * pi / m));
for(int i = 0;i < n;i += m) {
Complex w = Complex(1,0);
for(int k = 0;k < (m >> 1);++ k) {
Complex t = w * a[i + k + (m >> 1)],u = a[i + k];
a[i + k] = u + t;
a[i + k + (m >> 1)] = u - t;
w = w  * w1;
}
}
}
}
Complex ans[maxn];
int main() {
n = read(); int lim = 0;
for(int tmp,i = 1;i <= n;++ i) {
A[tmp] = Complex(1,0);
B[tmp * 2] = Complex(1,0);
C[tmp * 3] = Complex(1,0);
lim = std::max(lim,tmp * 3);
}
for(n = 1;n <= lim;n <<= 1);
fft(A,n,1);
fft(B,n,1);
fft(C,n,1);
for(int i = 0;i < n;++ i)
ans[i] = A[i] + (A[i] * A[i] - B[i]) * Complex(1.0 / 2.0,0) + (A[i] * A[i] * A[i] - A[i] * B[i] * Complex(3.0,0) + Complex(2.0,0) * C[i]) * Complex(1.0 / 6.0,0);
fft(ans,n,-1);
for(int i = 0;i < n;++ i) {
int T = ans[i].x / n + 0.5;
if(T) printf("%d %d\n",i,T);
}
return 0;
}


posted @ 2018-06-13 14:50  zzzzx  阅读(95)  评论(0编辑  收藏  举报