Codechef : CLOWAY(特征多项式+二项式反演)

传送门

题解:
先预处理出 Gi,k G i , k 表示第 i i 个点走k步回到自己的个数(可以停留)。

询问 (l,r,k) ( l , r , k ) 时,我们把 l,r l , r Gi,k G i , k 乘起来就得到了恰好 k k 步走回自己的概率。
不过注意这样走可能存在某一步一张图都没有走,这时候是不合法的。
不过我们发现一个只走了j步的方案会被统计 (kj) ( k j ) 次,这是一个二项式反演。用任意模数NTT即可。

现在考虑如何求 Gi,k G i , k ,显然我们可以直接矩阵快速幂,时间复杂度 O(n3k) O ( n 3 k )
暴力 O(n4) O ( n 4 ) 插值出特征多项式后就可以 O(nk) O ( n k ) 了。 时间复杂度为 O(T(n4+nk)+T2klogk) O ( T ∗ ( n 4 + n k ) + T 2 k log ⁡ k ) ,卡卡常就过了。

拆系数FFT

#pragma GCC optimize("Ofast")
#include <bits/stdc++.h>
using namespace std;

const int RLEN=1<<18|1;
inline char nc() {
    static char ibuf[RLEN],*ib,*ob;
    (ib==ob) && (ob=(ib=ibuf)+fread(ibuf,1,RLEN,stdin));
    return (ib==ob) ? -1 : *ib++;
}
inline int rd() {
    char ch=nc(); int i=0,f=1;
    while(!isdigit(ch)) {if(ch=='-')f=-1; ch=nc();}
    while(isdigit(ch)) {i=(i<<1)+(i<<3)+ch-'0'; ch=nc();}
    return i*f;
}
inline void W(int x) {
    static int buf[50];
    if(!x) {putchar('0'); return;}
    if(x<0) {putchar('-'); x=-x;}
    while(x) {buf[++buf[0]]=x%10; x/=10;}
    while(buf[0]) {putchar(buf[buf[0]--]+'0');}
}

typedef unsigned long long LL;
typedef __int128 IL; 
const int md=1e9+7; const LL LIM=(ULLONG_MAX/md-md)*md;
template <int mod> inline int add(int x,int y) {return (x+y>=mod) ? (x+y-mod) : (x+y);}
template <int mod> inline int dec(int x,int y) {return (x-y<0) ? (x-y+mod) : (x-y);}
template <int mod> inline int mul(int x,int y) {return (LL)x*y%mod;}
template <int mod> inline int power(int a,int b,int rs=1) {for(;b;b>>=1,a=mul<mod>(a,a)) if(b&1) rs=mul<mod>(rs,a); return rs;}
template <LL mod> inline LL quickmul (LL a,LL b) {return (IL)a*b%mod;}
const int N=55, K=8e4; int L=0;
int T,Q;
namespace FFT {
    const int mod1=998244353, mod2=1004535809, mod3=469762049, G=3;
    const int iM1=power <mod1> (mod2,mod1-2), iM2=power <mod2> (mod1,mod2-2);
    const LL M1=(LL)mod2*mod3, M2=(LL)mod1*mod3, M3=(LL)mod1*mod2;
    const int T3=power<mod3>(M3%mod3,mod3-2);
    int k,w[K],pos[K],tp[3][K],tp2[3][K];
    inline void init(int nn) {
        for(k=1;k<=nn;k<<=1);
        for(int i=1;i<k;i++) pos[i]=(i&1) ? ((pos[i>>1]>>1)^(k>>1)) : (pos[i>>1]>>1);
    }
    template <int mod> 
    inline void dft(int *a,int o=1) {
        for(int i=1;i<k;i++) if(pos[i]>i) swap(a[i],a[pos[i]]);
        for(int bl=1;bl<k;bl<<=1) {
            int tl=bl<<1, wn=power<mod>(G,(mod-1)/tl);
            w[0]=1; for(int i=1;i<bl;i++) w[i]=mul<mod>(w[i-1],wn);
            for(int bg=0;bg<k;bg+=tl)
                for(int j=0;j<bl;j++) {
                    int &t1=a[bg+j],&t2=a[bg+j+bl],t=mul<mod>(t2,w[j]);
                    t2=dec<mod>(t1,t); t1=add<mod>(t1,t);
                }
        }
        if(!~o) {
            const int inv=power <mod> (k,mod-2);
            reverse(a+1,a+k);
            for(int i=0;i<k;i++) a[i]=mul<mod>(a[i],inv);
        }
    } 
    template <int mod>
    inline void trans(int *a,int *b) {
        dft<mod>(a); dft<mod>(b);
        for(int i=0;i<k;i++) a[i]=mul<mod>(a[i],b[i]);
        dft<mod>(a,-1);
    }
    inline int merge(int a,int b,int c) {
        LL A=(quickmul<M3>((LL)iM1*mod2%M3,a)+quickmul<M3>((LL)iM2*mod1%M3,b))%M3; 
        LL ki=(c-A%mod3+mod3)%mod3; ki=T3%mod3*ki%mod3;
        A=(ki%md*(M3%md)%md+A%md)%md;
        return A;
    }
    inline void mul(int *a,int *b,int *c) {
        for(int i=0;i<=L;i++) 
            for(int j=0;j<3;j++) tp[j][i]=a[i], tp2[j][i]=b[i];
        for(int i=L+1;i<k;i++)
            for(int j=0;j<3;j++) tp[j][i]=tp2[j][i]=0;
        trans<mod1>(tp[0],tp2[0]);
        trans<mod2>(tp[1],tp2[1]);
        trans<mod3>(tp[2],tp2[2]);
        for(int i=0;i<=L;i++) c[i]=merge(tp[0][i],tp[1][i],tp[2][i]);
    }
}

int n,ans[K*3],G[N][N],F[N][K/7],Z[K];
int f[K],g[K],yc[K],fac[K],ifac[K]; 
struct Query {int l,r,k,id; Query(int l=0,int r=0,int k=0,int id=0):l(l),r(r),k(k),id(id){}} qry[K*3];
inline bool cmp(const Query &a,const Query &b) {return (a.l<b.l) || (a.l==b.l && a.r<b.r);}
struct Matrix {
    LL a[N][N];
    Matrix(int t=0) {
        memset(a,0,sizeof(a));
        for(int i=1;i<=n;i++) a[i][i]=t;
    }
    Matrix(int b[][N]) {
        for(int i=1;i<=n;i++)
            for(int j=1;j<=n;j++)
                a[i][j]=b[i][j];
    }
    friend inline Matrix operator *(const Matrix &a,const Matrix &b) {
        Matrix c;
        for(int i=1;i<=n;i++)
            for(int k=1;k<=n;k++) if(a.a[i][k])
                for(int j=1;j<=n;j++)
                    ((c.a[i][j]+=a.a[i][k]*b.a[k][j]) >= LIM) ? (c.a[i][j]%=md) : 0;
        for(int i=1;i<=n;i++)
            for(int j=1;j<=n;j++)
                if(c.a[i][j]>=md) c.a[i][j]%=md;
        return c;
    }
    inline int getval(int rs=0) {for(int i=1;i<=n;i++) rs=add<md>(rs,a[i][i]); return rs;}
} A[N];

inline int det(int a[][N],int v) {
    static int u[N][N];
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            u[i][j]=dec <md> (a[i][j],(i==j) ? v : 0);
    int rs=1, sgn=1;
    for(int i=1;i<=n;i++) {
        int l=i;
        for(;l<=n && !u[l][i];++l);
        if(l>n) return 0;
        if(l!=i) {
            sgn=md-sgn;
            for(int j=i;j<=n;j++) swap(u[l][j],u[i][j]);
        }
        const int inv=power<md>(u[i][i],md-2);
        for(int j=i+1;j<=n;j++) {
            int t=mul<md>(u[j][i],inv);
            for(int k=i;k<=n;k++)
                u[j][k]=dec<md>(u[j][k],mul<md>(u[i][k],t));
        }
    }
    for(int i=1;i<=n;i++) rs=mul<md>(rs,u[i][i]);
    return mul <md> (rs,sgn);
}

inline void getf() {
    memset(f,0,sizeof(f));
    for(int i=0;i<=n;i++) {
        static int tp[N]; 
        memset(tp,0,sizeof(tp)); tp[0]=1;
        for(int j=0;j<=n;j++) if(i!=j) {
            yc[i]=mul<md>(yc[i],power<md>(dec<md>(i,j),md-2));
            for(int k=n;k>=0;k--) {
                tp[k]=mul<md>(tp[k],md-j);
                if(k) tp[k]=add<md>(tp[k],tp[k-1]);
            }
        }
        for(int j=0;j<=n;j++) f[j]=add<md>(f[j],mul<md>(tp[j],yc[i]));
    }
}
inline void mul_g() {
    for(int i=n;i>=1;i--) g[i]=g[i-1]; g[0]=0;
    const int t=mul<md>(g[n],power<md>(f[n],md-2));
    for(int j=0;j<=n;j++) g[j]=dec<md>(g[j],mul<md>(t,f[j]));
}
struct Graph {
    int nn; int GG[N][N];
    inline void init() {
        nn=rd(); int m=rd();
        for(int i=1;i<=m;i++) {
            int x=rd(), y=rd();
            ++GG[x][y]; ++GG[y][x];
        }
    }
} gi[N];
inline void init(int id) {
    n=gi[id].nn; memset(G,0,sizeof(G));
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            G[i][j]=gi[id].GG[i][j];
    for(int i=1;i<=n;i++) ++G[i][i];
    A[0]=Matrix(1);
    for(int i=1;i<n;i++) A[i]=A[i-1]*Matrix(G);
    for(int i=0;i<n;i++) F[id][i]=A[i].getval();
    for(int i=0;i<=n;i++) yc[i]=det(G,i);
    getf(); 
    memset(g,0,sizeof(g)); g[0]=1;
    for(int i=1;i<n;i++) mul_g();
    for(int i=n;i<=L;i++) {
        mul_g(); int rs=0;
        for(int j=0;j<n;j++) rs=add <md> (rs,mul<md>(g[j],F[id][j]));
        F[id][i]=rs;
    }
}
inline void init(int l,int r) {
    static int gg[K/7],hh[K/7];
    for(int i=0;i<=L;i++) {
        gg[i]=1;
        for(int j=l;j<=r;j++) gg[i]=mul<md>(gg[i],F[j][i]);
        gg[i]=mul <md> (gg[i],ifac[i]);
    }
    hh[0]=1;
    for(int i=1;i<=L;i++) {
        hh[i]=ifac[i];
        if(i&1) hh[i]=md-hh[i];
    } 
    FFT::mul(gg,hh,gg);
    for(int i=0;i<=L;i++) Z[i]=mul <md> (gg[i],fac[i]);
    for(int i=2;i<=L;i++) Z[i]=add <md> (Z[i-1],Z[i]);
}
int main() {
    T=rd();
    for(int i=1;i<=T;i++) gi[i].init();
    Q=rd();
    for(int i=1;i<=Q;i++) {
        int l=rd(), r=rd(), k=rd(), id=i;
        qry[i]=Query(l,r,k,id); L=max(L,k);
    } 
    fac[0]=1;
    for(int i=1;i<=L;i++) fac[i]=mul <md> (fac[i-1],i);
    ifac[L]=power <md> (fac[L],md-2);
    for(int i=L-1;~i;i--) ifac[i]=mul <md> (ifac[i+1],i+1);
    for(int i=1;i<=T;i++) init(i);
    sort(qry+1,qry+Q+1,cmp);
    FFT::init(2*L);
    for(int i=1;i<=Q;i++) {
        if((qry[i].l!=qry[i-1].l) || (qry[i].r!=qry[i-1].r)) init(qry[i].l,qry[i].r);
        ans[qry[i].id]=Z[qry[i].k];
    }
    for(int i=1;i<=Q;i++) W(ans[i]), putchar('\n');
} 
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值