标签：round end AC complex rac command limits opera max
真是个愉悦的故事>_<
$\newcommand{\pretty}[1]{\begin{align*}#1\end{align*}}$设$\pretty{A(x)=\sum\limits_{i=1}^nx^{a_i},B(x)=\sum\limits_{i=1}^nx^{2a_i},C(x)=\sum\limits_{i=1}^nx^{3a_i}}$，那么只拿一个斧头，价值为$k$的的方案数是$\left[x^k\right]A(x)$，拿两个斧头，价值为$k$的方案数是$\left[x^k\right]\dfrac{A^2(x)-B(x)}2$（把“一把斧头拿两次”的情况容斥掉再消序），拿三把斧头，价值为$k$的方案数是$\left[x^k\right]\dfrac{A^3(x)-3A(x)B(x)+2C(x)}6$（容斥掉“一个斧头选两次，另一个斧头选一次”$\text{abb,bab,bba!}$和“一个斧头选三次”的情况，最后消序）
#include<stdio.h>
#include<math.h>
typedef double du;
typedef long long ll;
struct complex{
du x,y;
complex(du a=0,du 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);}
int rev[262144],N,iN;
void pre(int n){
int i,j,k;
for(N=1,k=0;N<n;N<<=1)k++;
for(i=0;i<N;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));
}
void swap(complex&a,complex&b){
complex c=a;
a=b;
b=c;
}
void fft(complex*a,int on){
int i,j,k;
complex t;
for(i=0;i<N;i++){
if(i<rev[i])swap(a[i],a[rev[i]]);
}
for(i=2;i<=N;i<<=1){
for(j=0;j<N;j+=i){
for(k=0;k<i>>1;k++){
t=complex(cos(k*M_PI/(i/2)),on*sin(k*M_PI/(i/2)))*a[i/2+j+k];
a[i/2+j+k]=a[j+k]-t;
a[j+k]=a[j+k]+t;
}
}
}
if(on==-1){
for(i=0;i<N;i++)a[i].x/=N;
}
}
int max(int a,int b){return a>b?a:b;}
complex a[262144],b[262144],c[262144];
int main(){
int n,m,i,x;
ll s;
scanf("%d",&n);
m=0;
for(i=0;i<n;i++){
scanf("%d",&x);
a[x].x+=1;
b[x*2].x+=1;
c[x*3].x+=1 ;
m=max(m,x*3);
}
pre(m<<1|1);
fft(a,1);
fft(b,1);
fft(c,1);
for(i=0;i<N;i++)a[i]=(a[i]*a[i]*a[i]-3*a[i]*b[i]+2*c[i])*(1/6.)+(a[i]*a[i]-b[i])*.5+a[i];
fft(a,-1);
for(i=1;i<=m;i++){
s=llround(a[i].x);
if(s)printf("%d %lld\n",i,s);
}
}
标签：round end AC complex rac command limits opera max
原文地址：https://www.cnblogs.com/jefflyy/p/8893410.html
