BZOJ 3456 城市规划 ——NTT
搞出递推式。
发现可以变成三个函数的乘积。
移项之后就可以求逆+NTT做了。
miskoo博客中有讲
#include <map> #include <cmath> #include <queue> #include <cstdio> #include <cstring> #include <iostream> #include <algorithm> using namespace std; #define F(i,j,k) for (int i=j;i<=k;++i) #define D(i,j,k) for (int i=j;i>=k;--i) #define ll long long #define mp make_pair #define md 1004535809 #define g 3 #define maxn 500005 int rev[maxn],n; int ksm(int a,int b) { int ret=1; for (;b;b>>=1,a=(ll)a*a%md) if (b&1) ret=(ll)ret*a%md; return ret; } void NTT(int *x,int n,int flag) { F(i,0,n-1) if (rev[i]>i) swap(x[rev[i]],x[i]); for (int m=2;m<=n;m<<=1) { int wn=ksm(g,((md-1)/m*flag+md-1)%(md-1)); for (int i=0;i<n;i+=m) { int w=1; for (int j=0;j<(m>>1);++j) { int u=x[i+j],v=(ll)x[i+j+(m>>1)]*w%md; x[i+j]=(u+v)%md; x[i+j+(m>>1)]=(u-v+md)%md; w=(ll)w*wn%md; } } } if (flag==-1) { int inv=ksm(n,md-2); F(i,0,n-1) x[i]=(ll)x[i]*inv%md; } } int fac[maxn],fac_inv[maxn],C[maxn],G[maxn],F[maxn],N,Inv_G[maxn]; void Get_Inv(int *a,int *b,int n) { static int tmp[maxn];if (n==1){b[0]=ksm(a[0],md-2);return;} Get_Inv(a,b,n>>1);F(i,0,n-1)tmp[i]=a[i],tmp[i+n]=0; int L=0;while(!(n>>L&1))L++; F(i,0,(n<<1)-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<L); NTT(tmp,n<<1,1);NTT(b,n<<1,1); F(i,0,(n<<1)-1) tmp[i]=(ll)b[i]*(2LL-(ll)tmp[i]*b[i]%md+md)%md; NTT(tmp,n<<1,-1);F(i,0,n-1) b[i]=tmp[i],b[n+i]=0; } int main() { scanf("%d",&n); fac[0]=1;F(i,1,maxn-1) fac[i]=(ll)fac[i-1]*i%md; fac_inv[maxn-1]=ksm(fac[maxn-1],md-2); D(i,maxn-2,0) fac_inv[i]=(ll)fac_inv[i+1]*(i+1)%md; for (N=1;N<=n;N<<=1); F(i,0,n) C[i]=(ll)ksm(2,(ll)i*(i-1)/2%(md-1))*fac_inv[i-1]%md; F(i,0,n) G[i]=(ll)ksm(2,(ll)i*(i-1)/2%(md-1))*fac_inv[i]%md; Get_Inv(G,Inv_G,N); NTT(C,N<<1,1);NTT(Inv_G,N<<1,1); F(i,0,(N<<1)-1) F[i]=(ll)C[i]*Inv_G[i]%md; NTT(F,N<<1,-1); printf("%d\n",(ll)F[n]*fac[n-1]%md); }
Po姐讲了另外一种方法。
哈哈哈,完全不会,抄抄抄
#include <map> #include <cmath> #include <queue> #include <cstdio> #include <cstring> #include <iostream> #include <algorithm> using namespace std; #define F(i,j,k) for (int i=j;i<=k;++i) #define D(i,j,k) for (int i=j;i>=k;--i) #define ll long long #define mp make_pair #define maxn 500005 #define md 1004535809 #define g 3 int rev[maxn]; int ksm(int a,int b) { int ret=1; for (;b;b>>=1,a=(ll)a*a%md) if (b&1) ret=(ll)ret*a%md; return ret; } void NTT(int *x,int n,int flag) { F(i,0,n-1) if (rev[i]>i) swap(x[rev[i]],x[i]); for (int m=2;m<=n;m<<=1) { int wn=ksm(g,((md-1)/m*flag+md-1)%(md-1)); for (int i=0;i<n;i+=m) { int w=1; for (int j=0;j<(m>>1);++j) { int u=x[i+j],v=(ll)x[i+j+(m>>1)]*w%md; x[i+j]=(u+v)%md; x[i+j+(m>>1)]=(u-v+md)%md; w=(ll)w*wn%md; } } } if (flag==-1) { int inv=ksm(n,md-2); F(i,0,n-1) x[i]=(ll)x[i]*inv%md; } } int n,G[maxn],F[maxn],Inv_G[maxn],N,fac[maxn],fac_inv[maxn],Der_G[maxn]; void Get_Inv(int *a,int *b,int n) { static int tmp[maxn];if (n==1){b[0]=ksm(a[0],md-2);return;} Get_Inv(a,b,n>>1);F(i,0,n-1)tmp[i]=a[i],tmp[n+i]=0; int L=0;while(!(n>>L&1)) L++; F(i,0,(n<<1)-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<L); NTT(tmp,n<<1,1);NTT(b,n<<1,1); F(i,0,(n<<1)-1) tmp[i]=(ll)b[i]*(2-(ll)tmp[i]*b[i]%md+md)%md; NTT(tmp,n<<1,-1); F(i,0,n-1) b[i]=tmp[i],b[i+n]=0; } int main() { scanf("%d",&n); fac[0]=1;F(i,1,maxn-1) fac[i]=(ll)fac[i-1]*i%md; fac_inv[maxn-1]=ksm(fac[maxn-1],md-2); for (N=1;N<=n;N<<=1); D(i,maxn-2,0)fac_inv[i]=(ll)fac_inv[i+1]*(i+1)%md; F(i,0,n)G[i]=(ll)ksm(2,(ll)i*(i-1)/2%(md-1))*fac_inv[i]%md; Get_Inv(G,Inv_G,N); F(i,1,N-1) Der_G[i-1]=(ll)G[i]*i%md; int L=0;while(!(N>>L&1)) L++; F(i,0,(N<<1)-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<L); NTT(Der_G,N<<1,1);NTT(Inv_G,N<<1,1); F(i,0,(N<<1)-1)F[i]=(ll)Der_G[i]*Inv_G[i]%md; NTT(F,N<<1,-1); printf("%d\n",(ll)F[n-1]*ksm(n,md-2)%md*fac[n]%md); }