日常学习——FFT
FFT相关详解：
http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform
https://blog.csdn.net/ggn_2015/article/details/68922404
练习题：
模板题
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define inf 0x7f7f7f7f
using namespace std;
const int maxn=2e5;
const double pi=acos(-1);
int n,m,l;
int sum[maxn*2+10],rev[maxn+10];
int read()
{
char ch;int x=0,f=1;
for (;ch<'0'||ch>'9';ch=getchar()) if (ch=='-') f=-1;
for (;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0';
return x*f;
}
struct complex
{
double x,y;
complex(){};
complex(double _x,double _y){x=_x,y=_y;}
complex operator +(const complex &a){return complex(x+a.x,y+a.y);}
complex operator -(const complex &a){return complex(x-a.x,y-a.y);}
complex operator *(const complex &a){return complex(x*a.x-y*a.y,x*a.y+y*a.x);}
}A[maxn+10],B[maxn+10],Ans[maxn*2+10];
int get()
{
int ch=getchar();
for (;ch<'0'||ch>'9';ch=getchar());
return ch-'0';
}
void fft(complex *a,int len,int flag)
{
for (int i=0;i<len;i++) if (rev[i]<i) swap(a[i],a[rev[i]]);
for (int i=2;i<=len;i<<=1)
{
complex wn(cos(2*pi/i),sin(2*pi/i*flag));
for (int j=0;j<len;j+=i)
{
complex nw(1,0);
for (int k=0;k<(i>>1);k++,nw=nw*wn)
{
complex x=a[j+k],y=nw*a[j+k+(i>>1)];
a[j+k]=x+y,a[j+k+(i>>1)]=x-y;
}
}
}
if (flag==-1) for (int i=0;i<len;i++) a[i].x/=len;
}
int main()
{
n=read();
for (int i=n-1;~i;i--) A[i].x=get();
for (int i=n-1;~i;i--) B[i].x=get();
for (n*=2,m=1;m<=n;m<<=1)l++;
for (int i=0;i<m;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
fft(A,m,1),fft(B,m,1);
for (int i=0;i<m;i++) Ans[i]=A[i]*B[i];
fft(Ans,m,-1);
for (int i=0;i<m;i++) sum[i]=(int)(Ans[i].x+0.5);
for (int i=0;i<m;i++)
{
sum[i+1]+=sum[i]/10;
sum[i]%=10;
}
while((!sum[n-1])&&n) n--;
while(sum[n]) n++;
for (int i=n-1;~i;i--) putchar(sum[i]+'0');
printf("\n");
return 0;
}
模板题+1。计算\(C_k=\sum_{i=k}^{n-1}{a_i \times b_{i-k}}\),将\(b[]\)翻转，则\(C_k=\sum_{i=k}^{n-1}{a_i \times b_{n+1+k-i}}\)，即\(C_k=\sum\limits_{i+j=n+1+k,0<=i,j<=2 \times n}{a_i \times b_j}\)。直接FFT即可。
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define inf 0x7f7f7f7f
using namespace std;
const int maxn=4e5;
const double pi=acos(-1);
int n,m,l;
int rev[(maxn<<2)+5];
int read()
{
char ch;int x=0,f=1;
for (;ch<'0'||ch>'9';ch=getchar()) if (ch=='-') f=-1;
for (;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0';
return x*f;
}
struct complex
{
double x,y;
complex(){};
complex(double _x,double _y){x=_x,y=_y;}
complex operator +(const complex &a){return complex(x+a.x,y+a.y);};
complex operator -(const complex &a){return complex(x-a.x,y-a.y);};
complex operator *(const complex &a){return complex(x*a.x-y*a.y,x*a.y+y*a.x);};
}a[maxn+10],b[maxn+10],ans[(maxn<<2)+5];
void fft(complex *a,int len,int flag)
{
for (int i=0;i<len;i++) if (rev[i]<i) swap(a[i],a[rev[i]]);
for (int i=2;i<=len;i<<=1)
{
complex wn(cos(2*pi/i),sin(2*pi/i*flag));
for (int j=0;j<len;j+=i)
{
complex nw(1,0);
for (int k=0;k<(i>>1);k++,nw=nw*wn)
{
complex x=a[j+k],y=nw*a[j+k+(i>>1)];
a[j+k]=x+y,a[j+k+(i>>1)]=x-y;
}
}
}
if (flag==-1) for (int i=0;i<len;i++) a[i].x/=len;
}
int main()
{
n=read();
for (int i=0;i<n;i++) a[i].x=read(),b[n+1-i].x=read();
for (m=1;m<=n*2;m<<=1)l++;
for (int i=0;i<m;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
fft(a,m,1),fft(b,m,1);
for (int i=0;i<m;i++) ans[i]=a[i]*b[i];
fft(ans,m,-1);
for (int i=n+1;i<=n*2;i++) printf("%d\n",(int)(ans[i].x+0.5));
return 0;
}
假设移动a位，整体增加c，则答案为\(\sum_{i=1}^{n}{(x_{((i+a-1) \% n+1)}-y_{i}+c)^{2}}\)。把此式展开\(\sum _{i=1} ^{n} {x_{( (i+a-1) \% n+1 )} ^{2} +y_{i} ^{2} +c ^{2} +2 \times c \times x_{( (i+a-1) \% n+1 )} -2 \times y_{i} \times c -2 \times x_{( (i+a-1) \% n+1)} \times y_{i}}\)。
\(\sum _{i=1} ^{n} {x_{( (i+a-1) \% n+1 )} \times y_{i}}\)可以用FFT求，之后枚举a和c就行了。
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define ll long long
#define inf 0x7f7f7f7f
using namespace std;
const int maxn=5e4;
const double pi=acos(-1);
int n,m,k,l,s,ans=1e9;
int a[maxn+10],b[maxn+10],rev[maxn*4+10];
int read()
{
char ch;int x=0,f=1;
for (;ch<'0'||ch>'9';ch=getchar()) if (ch=='-') f=-1;
for (;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0';
return x*f;
}
struct complex
{
double x,y;
complex(){};
complex (double _x,double _y){x=_x,y=_y;};
complex operator +(const complex &a){return complex(x+a.x,y+a.y);}
complex operator -(const complex &a){return complex(x-a.x,y-a.y);}
complex operator *(const complex &a){return complex(x*a.x-y*a.y,x*a.y+y*a.x);}
}x1[maxn*4+10],x2[maxn*4+10];
void fft(complex *a,int len,int flag)
{
for (int i=0;i<len;i++) if (rev[i]<i) swap(a[i],a[rev[i]]);
for (int i=2;i<=len;i<<=1)
{
complex wn(cos(2*pi/i),sin(2*pi/i*flag));
for (int j=0;j<len;j+=i)
{
complex nw(1,0);
for (int k=0;k<(i>>1);k++,nw=nw*wn)
{
complex x=a[j+k],y=nw*a[j+k+(i>>1)];
a[j+k]=x+y,a[j+k+(i>>1)]=x-y;
}
}
}
if (flag==-1) for (int i=0;i<len;i++) a[i].x/=len;
}
int main()
{
n=read(),k=read();
for (int i=0;i<n;i++) a[i]=read(),x1[i].x=a[i];
for (int i=0;i<n;i++) b[i]=read(),x2[n-i-1].x=b[i],s+=b[i];
for (m=1;m<=n*2;m<<=1)l++;
for (int i=0;i<m;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
fft(x1,m,1),fft(x2,m,1);
for (int i=0;i<m;i++) x1[i]=x1[i]*x2[i];
fft(x1,m,-1);
for (int c=0;c<=k;c++)
{
int sum=0;
for (int i=0;i<n;i++) sum+=(a[i]+c)*(a[i]+c)+b[i]*b[i];
ans=min(ans,sum-2*((int)(x1[n-1].x+0.5)+s*c));
for (int i=1;i<n;i++)
{
int tmp=(int)(x1[n+i-1].x+0.5)+(int)(x1[i-1].x+0.5)+s*c;
ans=min(ans,sum-2*tmp);
}
}
printf("%d\n",ans);
return 0;
}
bzoj4555 [Tjoi2016&Heoi2016]求和
求\(f_n=\sum_{i=0}^{n}\sum_{j=0}^{i}S_{i,j} \times 2^j \times j!\),对答案\(\%98244353\)。其中\(S_{i,j}\)为第二类斯特林数，意义为将n个不同的元素拆分成m个集合的方案数,所以当m>n时\(S_{n,m}\)为0。而\(S_{n,m} \times m!\)就是将n个不同的元素拆分成m个不同的集合的方案数。用容斥可得到\(S_{n,m} \times m!=\sum_{i=0}^{m}{(-1)^i \times \binom{m}{i} \times (m-i)^n}\)。
好，那我们开始推式子。
\(f_n=\sum_{i=0}^{n}\sum_{j=0}^{i}S_{i,j} \times 2^j \times j!\)
\(f_n=\sum_{i=0}^{n}\sum_{j=0}^{n}{2^j \times\sum_{k=0}^{j}}(-1)^k \times \binom{j}{k} \times (j-k)^i\)
\(f_n=\sum_{j=0}^{n}2^j \times \sum_{k=0}^{j}(-1)^k \times \frac{j!}{k!(j-k)!} \times \sum_{i=0}^{n} (j-k)^i\)
\(f_n=\sum_{j=0}^{n} 2^j \times j! \times \sum_{k=0}^{j} \frac{(-1)^k}{k!} \times \frac{\sum_{i=0}^{n}(j-k)^i}{(j-k)!}\)
后面的\(\sum_{k=0}^{j} \frac{(-1)^k}{k!} \times \frac{\sum_{i=0}^{n}(j-k)^i}{(j-k)!}\)就是一个卷积的形式，因为要膜，所以直接上NTT即可。
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define ll long long
#define inf 0x7f7f7f7f
using namespace std;
const int mod=998244353;
const int maxn=1e5+5;
int ans,n,m,l;
int rev[maxn<<2],inv[maxn<<2],x[maxn<<2],y[maxn<<2];
int read()
{
char ch;int x=0,f=1;
for (;ch<'0'||ch>'9';ch=getchar()) if (ch=='-') f=-1;
for (;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0';
return x*f;
}
int power(int a,int k)
{
int sum=1;
for (;k;k>>=1,a=1ll*a*a%mod)
if (k&1)
sum=1ll*sum*a%mod;
return sum;
}
void ntt(int *a,int len,int flag)
{
for (int i=0;i<len;i++) if (rev[i]<i) swap(a[rev[i]],a[i]);
for (int i=2;i<=len;i<<=1)
{
int gn=power(3,((mod-1)+flag*(mod-1)/i)%(mod-1));
for (int j=0;j<len;j+=i)
{
int ng=1;
for (int k=0;k<(i>>1);k++,ng=1ll*ng*gn%mod)
{
int x=a[j+k],y=1ll*ng*a[j+k+(i>>1)]%mod;
a[j+k]=(x+y)%mod,a[j+k+(i>>1)]=1ll*((x-y)%mod+mod)%mod;
}
}
}
if (flag==-1)
{
int invlen=power(len,mod-2);
for (int i=0;i<len;i++) x[i]=1ll*x[i]*invlen%mod;
}
}
int main()
{
n=read();
for (m=1;m<n*2;m<<=1)l++;
for (int i=0;i<m;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
inv[0]=1;
for (int i=1;i<=n;i++) inv[i]=1ll*inv[i-1]*i%mod;
inv[n]=power(inv[n],mod-2);
for (int i=n-1;i;i--) inv[i]=1ll*inv[i+1]*(i+1)%mod;
for (int i=0;i<=n;i++) x[i]=(inv[i]*(i&1?-1:1)+mod)%mod;
y[0]=1,y[1]=n+1;
for (int i=2;i<=n;i++) y[i]=1ll*(power(i,n+1)-1)*power(i-1,mod-2)%mod*inv[i]%mod;
ntt(x,m,1),ntt(y,m,1);
for (int i=0;i<m;i++) x[i]=1ll*x[i]*y[i]%mod;
ntt(x,m,-1);
int two_power=1,fact=1;ans=x[0];
for (int i=1;i<=n;i++)
{
two_power=1ll*two_power*2%mod,fact=1ll*fact*i%mod;
ans=(ans+1ll*two_power*fact%mod*x[i]%mod)%mod;
}
printf("%d\n",ans);
return 0;
}
