[BZOJ3456]城市规划

题链

sol：  Fn=Wn−(n−1)!×∑i=1nFi(i−1)!∗Wn−i(n−i)!Fn=Wn−(n−1)!×∑i=1nFi(i−1)!∗Wn−i(n−i)

 我们可以分治FFT

#pragma GCC optimize("-O3")
#include<bits/stdc++.h>
#define mo 1004535809
#define LL long long
#define N 270301
using namespace std;
int D[N],L,m,n;
LL gg,aa[2007],ni[2007],g[N],nig[N],f[N],x[N],y[N],z[N],wn,w,X,Y,NI,b[N],a[N],nia[N];
inline LL qsm(LL x,LL y) {
    static LL anw;
    for (anw=1;y;y>>=1,x=x*x%mo) if (y&1) anw=anw*x%mo; 
    return anw;
}
inline void Mo(LL &x) {
    if (x<0) x=x%mo+mo; else if (x>=mo) x%=mo; 
}
inline void FNT(LL* a,int n,int y){
    for (int i=0;i<n;i++) {
        D[i]=(D[i>>1]>>1)|((i&1)<<(L-1));
        if (D[i]>i) swap(a[i],a[D[i]]);
    }
    for (int i=1;i<n;i<<=1) {
        wn=y?g[i]:nig[i];
        for (int j=0,w=1;j<n;j+=i<<1,w=1) 
            for (int k=j;k<j+i;k++,w=w*wn%mo) {
                X=a[k]; Y=w*a[k+i]%mo;
                a[k]=(X+Y)%mo;a[k+i]=(X-Y+mo)%mo;
        }
    }
  NI=ni[L]; if (!y) for (int i=0;i<n;i++) (a[i]*=NI)%=mo;
}
void FFT(int n){
    for(m=1,L=1;m<=n;m<<=1,L++); m<<=1;
    FNT(x,m,1); FNT(y,m,1);
    for (int i=0;i<m;i++) z[i]=(x[i]*y[i])%mo;
    FNT(z,m,0);
}
#define Mid (l+r>>1)
#define max(a,b) ((a)>(b)?(a):(b))
void cqh(int l,int r){
    if (l==r) {f[l]=(b[l]+mo-(a[l-1]*f[l])); Mo(f[l]); return;}
    cqh(l,Mid);
    int cnt=0;
    for (int i=0;i<r-l;i++) x[cnt++]=(b[i+1]*nia[i+1])%mo;
    for (int i=Mid-l;~i;i--) y[i]=(f[l+i]*nia[l+i-1])%mo;
    FFT(max(r-Mid,Mid-l+1));
    for (int i=Mid+1;i<=r;i++) f[i]+=z[i-l-1],Mo(f[i]);
    memset(x,0,gg*(m+1)),memset(y,0,gg*(m+1)),memset(z,0,gg*(m+1));
    cqh(Mid+1,r);
}
signed main () {
    gg=sizeof gg;
  freopen("bzoj_3456.in","r",stdin);
  freopen("bzoj_3456.out","w",stdout);
    cin>>n; nia[0]=a[0]=ni[0]=aa[0]=1;
    for (int i=1;i<=n;i++) a[i]=a[i-1]*i%mo; 
    nia[n]=qsm(a[n],mo-2);for (int i=n;i;i--) nia[i-1]=nia[i]*i%mo;
    for (int i=1;i<=n;i++) g[i]=qsm(3,(mo-1)/i/2),nig[i]=qsm(g[i],mo-2);
    for (int i=0;i<=n;i++) b[i]=qsm(2,1ll*i*(i-1)/2);
    for (int i=1;i<=1000;i++) aa[i]=aa[i-1]*2%mo;
    ni[1000]=qsm(aa[1000],mo-2);
    for (int i=1000;i;i--) ni[i-1]=(ni[i]<<1)%mo;
    cqh(1,n);
    cout<<f[n]; return 0;
}

 我们还可以 多项式求逆。

#pragma GCC optimize("-O2")
#include<bits/stdc++.h>
#define LL long long
#define mo 1004535809
#define N 300005
//#define int LL
using namespace std;
LL jc[N],ny[N],b[N],c[N],ni_b[N],tmp[N];
int rev[N],n,m,L;
inline LL qsm(LL x,LL y=mo-2){
    static LL anw;
    for (anw=1;y;y>>=1,x=x*x%mo) if (y&1) anw=anw*x%mo;
    return anw;
}
void NTT(LL *a,int x,int n,int L){
    for (int i=0;i<n;i++) {
        rev[i]=(rev[i>>1]>>1)|((i&1)<<L-1);
        if (i<rev[i]) swap(a[i],a[rev[i]]);
    }
    LL X,Y;
    for (int i=1;i<n;i<<=1) {
        LL wn=qsm(3,(mo-1)/i/2);
        if (!(~x)) wn=qsm(wn);
        for (int j=0;j<n;j+=i<<1)  {
            LL w=1;
            for (int k=0;k<i;k++,w=w*wn%mo) {
                X=a[j+k]; Y=w*a[j+k+i]%mo;
                a[j+k]=(X+Y)%mo; a[j+k+i]=((X-Y)%mo+mo)%mo;
            } 
        }
    }
    if (!(~x)) {
        LL wc=qsm(n);
        for (int i=0;i<n;i++) a[i]=a[i]*wc%mo; }
}
void get(LL* a,LL* b,int n,int L) {
    if (n==1) {
        b[0]=qsm(a[0]); return;
    }
    get(a,b,n>>1,L-1);
    memcpy(tmp,a,sizeof (a[0])*n);
    memset(tmp+n,0,sizeof (a[0])*n);
    NTT(tmp,1,n<<1,L+1); NTT(b,1,n<<1,L+1);
    for (int i=(n<<1)-1;~i;--i)
      b[i]=(2+mo-tmp[i]*b[i]%mo)%mo*b[i]%mo;
    NTT(b,-1,n<<1,L+1);
    memset(b+n,0,sizeof (b[0])*n);
}
signed main () {
//    freopen("a.in","r",stdin);
    scanf("%d",&n);
    jc[0]=ny[0]=1;
    for (int i=1;i<=n;i++) jc[i]=jc[i-1]*i%mo,ny[i]=qsm(jc[i]);
    for (int i=0;i<=n;i++) b[i]=qsm(2,(LL)i*(i-1)/2)*ny[i]%mo;
    for (int i=1;i<=n;i++) c[i]=qsm(2,(LL)i*(i-1)/2)*ny[i-1]%mo;
    for (m=1;m<=n;m<<=1) L++;
    get(b,ni_b,m,L);
    NTT(c,1,m<<1,L+1); NTT(ni_b,1,m<<1,L+1);
    for (int i=(m<<1)-1;~i;i--) c[i]=c[i]*ni_b[i]%mo;
    NTT(c,-1,m<<1,L+1);
    printf("%lld
",(LL)c[n]*jc[n-1]%mo);
    return 0;
}