Codechef CLOWAY Future of draughts
Link
可以发现不走其实就是走自环，所以我们给每张图的每个点添加一个自环。
对于第\(i\)张图，设其邻接矩阵为\(\mathbf A_i\)，那么在第\(i\)张图上走\(k\)步回到原点的方案数就是\(\operatorname{tr}(\mathbf A_i^k)\)。
设\(g_{l,r,k}=\prod\limits_{i=l}^r\operatorname{tr}(\mathbf A_i^k)\)，\(f_{l,r,k}\)表示在\(G_l,\cdots,G_r\)进行\(k\)轮移动的方案数。
注意到\(g\)可能存在一轮全都走自环的情况，所以我们枚举没有全都走自环的轮数得到\(g_{l,r,k}=\sum\limits_{i=0}^k{k\choose i}f_{l,r,i}\)。
二项式反演得到\(f_{l,r,k}=\sum\limits_{i=0}^k(-1)^{k-1}{k\choose i}g_{l,r,i}\)。
注意到这是个卷积的形式，我们可以拆组合数将其化为\(f_{l,r,k}=k!\sum\limits_{i+j=k}\frac{g_{l,r,i}}{i!}\frac{(-1)^j}{j!}\)。
因此如果我们能够求出所有的\(g\)，那么我们就能在\(O(T^2K\log K)\)的时间复杂度内求出所有\(f\)，进而\(O(1)\)地回答每个询问。
要求所有的\(g_{l,r,k}\)就是要求所有的\(g_{i,i,k}=\operatorname{tr}(\mathbf A_i^k)\)，直接做是\(O(TN^3K)\)的，考虑优化。
根据Cayley-Hamilton定理，对于\(\mathbf A_i\)，我们求出其特征多项式\(f_i(\lambda)=|\lambda \mathbf E-\mathbf A_i|\)，同时求出所有的\(h_{i,k}(\lambda)=\lambda^k\bmod f_i(\lambda)\)，那么\(\operatorname{tr}(\mathbf A_i^k)=\sum\limits_{j=0}^{n_i-1}[x^j]h_{i,k}(\lambda)\operatorname{tr}(\mathbf A_i^j)\)。
特征多项式可以\(O(N^4)\)插值也可以用对角化+上Hessenberg矩阵做到\(O(N^3)\)，这里用的是插值。
因为我们是要求出所有\(k\)的\(h_{i,k}\)，所以可以\(O(NK)\)递推。
然后求\(j\in[0,n_i-1)\)的\(\operatorname{tr}(\mathbf A_i^j)\)直接\(O(N^4)\)暴力就好了。
总时间复杂度为\(O(T(N^4+NK)+T^2K\log K+Q)\)。
因为模数是\(10^9+7\)所以要用MTT，码量增加的同时常数也大大增加了，需要加一些小小的常数优化。
#include<cmath>
#include<cstdio>
#include<cctype>
#include<cstring>
#include<algorithm>
namespace IO
{
char ibuf[(1<<21)+1],obuf[(1<<21)+1],st[15],*iS,*iT,*oS=obuf,*oT=obuf+(1<<21);
char Get(){return (iS==iT? (iT=(iS=ibuf)+fread(ibuf,1,(1<<21)+1,stdin),(iS==iT? EOF:*iS++)):*iS++);}
void Flush(){fwrite(obuf,1,oS-obuf,stdout),oS=obuf;}
void Put(char x){*oS++=x;if(oS==oT)Flush();}
int read(){int x=0,c=Get();while(!isdigit(c))c=Get();while(isdigit(c))x=x*10+c-48,c=Get();return x;}
}
using IO::read;
using ld=long double;
const int P=1000000007;
const ld pi=acos(-1);
int inc(int a,int b){return a+=b-P,a+(a>>31&P);}
int dec(int a,int b){return a-=b,a+(a>>31&P);}
int mul(int a,int b){return 1ll*a*b%P;}
int pow(int a,int k){int r=1;for(;k;k>>=1,a=mul(a,a))if(k&1)r=mul(a,r);return r;}
int inv(int a){return pow(a,P-2);}
struct matrix{int a[51][51],n;matrix(){memset(a,0,sizeof a);}int*operator[](int x){return a[x];}};
matrix operator+(matrix a,matrix b){for(int i=1;i<=a.n;++i)for(int j=1;j<=a.n;++j)a[i][j]=inc(a[i][j],b[i][j]);return a;}
matrix operator-(matrix a,matrix b){for(int i=1;i<=a.n;++i)for(int j=1;j<=a.n;++j)a[i][j]=dec(a[i][j],b[i][j]);return a;}
matrix operator*(matrix a,matrix b){matrix c;c.n=a.n;for(int i=1;i<=a.n;++i)for(int j=1;j<=a.n;++j)for(int k=1;k<=a.n;++k)c[i][j]=inc(c[i][j],mul(a[i][k],b[k][j]));return c;}
matrix I(int n){matrix a;a.n=n;for(int i=1;i<=n;++i)a[i][i]=1;return a;}
int tr(matrix a){int r=0;for(int i=1;i<=a.n;++i)r=inc(r,a[i][i]);return r;}
int det(matrix a)
{
int r=1;
for(int i=1;i<=a.n;++i)
{
int f=0;
for(int j=i;j<=a.n;++j)
if(a[j][i])
{
if(!f)
{
if(f=1,i^j) r=dec(0,r),std::swap(a.a[j],a.a[i]);
r=mul(r,a[i][i]);
for(int k=i,x=inv(a[i][i]);k<=a.n;++k) a[i][k]=mul(a[i][k],x);
}
else for(int k=a.n;k>=i;--k) a[j][k]=dec(a[j][k],mul(a[i][k],a[j][i]));
}
if(!f) return 0;
}
return r;
}
void work(matrix a,int*f)
{
static int r[51],t[51];int n=a.n;memset(f,0,(n+1)<<2);
for(int i=1;i<=n;++i) for(int j=1;j<=n;++j) a[i][j]=dec(0,a[i][j]);
for(int i=0;i<=n;++i) r[i]=det(a=a+I(n));
for(int i=0,s;i<=n;++i)
{
memset(t,0,(n+1)<<2),t[0]=1,s=r[i];
for(int j=0,k;j<=n;++j) if(i^j) for(s=mul(s,inv(dec(i,j))),k=n;~k;--k) t[k]=inc(mul(P-j-1,t[k]),k? t[k-1]:0);
for(int j=0;j<=n;++j) f[j]=inc(f[j],mul(s,t[j]));
}
}
struct complex{ld x,y;complex(ld a=0,ld b=0){x=a,y=b;}}w[32769];
complex operator+(complex a,complex b){return {a.x+b.x,a.y+b.y};}
complex operator-(complex a,complex b){return {a.x-b.x,a.y-b.y};}
complex operator*(complex a,ld x){return {a.x*x,a.y*x};}
complex operator/(complex a,ld x){return {a.x/x,a.y/x};}
complex operator*(complex a,complex b){return {a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
complex conj(complex a){return {a.x,-a.y};}
int rev[32769],lim;
void init(int n)
{
static ld x;
lim=1<<(32-__builtin_clz(n));
for(int i=1;i<lim;++i) rev[i]=(rev[i>>1]>>1)|(i&1? lim>>1:0);
for(int i=0;i<lim;++i) x=pi*i/lim,w[i]={cos(x),sin(x)};
}
void FFT(complex*a,int f)
{
static complex x;
if(!~f) std::reverse(a+1,a+lim);
for(int i=1;i<lim;++i) if(i<rev[i]) std::swap(a[i],a[rev[i]]);
for(int i=1;i<lim;i<<=1) for(int l=i<<1,j=0;j<lim;j+=l) for(int k=0;k<i;++k) x=a[i+j+k]*w[lim/i*k],a[i+j+k]=a[j+k]-x,a[j+k]=a[j+k]+x;
if(!~f) for(int i=0;i<lim;++i) a[i]=a[i]/lim;
}
void DFT(complex*a,complex*b)
{
static complex t[2][32769],x,y;
for(int i=0;i<lim;++i) t[0][i]=a[i]+b[i]*complex{0,1};
FFT(t[0],1),memcpy(t[1],t[0],lim<<5),std::reverse(t[1]+1,t[1]+lim);
for(int i=0;i<lim;++i) x=t[0][i],y=conj(t[1][i]),a[i]=(x+y)/2,b[i]=(x-y)*complex{0,-1}/2;
}
void IDFT(complex*a,complex*b)
{
for(int i=0;i<lim;++i) a[i]=a[i]+b[i]*complex{0,1};
FFT(a,-1);
for(int i=0;i<lim;++i) b[i]={a[i].y,0},a[i].y=0;
}
void MTT(int*a,int*b,int*c,int n,int m)
{
static complex t[4][32769],A,B,C,D;static int r[4]={1<<30,1<<15,1<<15,1};
init(n+m-1),memset(t,0,sizeof t);
for(int i=0;i<n;++i) t[0][i].x=a[i]>>15,t[1][i].x=a[i]&32767;
for(int i=0;i<m;++i) t[2][i].x=b[i]>>15,t[3][i].x=b[i]&32767;
DFT(t[0],t[1]),DFT(t[2],t[3]);
for(int i=0;i<lim;++i) A=t[0][i]*t[2][i],B=t[0][i]*t[3][i],C=t[1][i]*t[2][i],D=t[1][i]*t[3][i],t[0][i]=A,t[1][i]=B,t[2][i]=C,t[3][i]=D;
IDFT(t[0],t[1]),IDFT(t[2],t[3]),memset(c,0,(n+m)<<2);
for(int i=0;i<4;++i) for(int j=0;j<n+m;++j) c[j]=inc(c[j],mul(r[i],(long long)(t[i][j].x+0.5)%P));
}
int f[51][51][10001],g[51][51][10001],fac[10001],ifac[10001],mx[51][51],T,Q;
struct query{int l,r,k;}q[200001];
int main()
{
T=read(),fac[0]=1;
for(int i=1;i<=10000;++i) fac[i]=mul(fac[i-1],i);
ifac[10000]=inv(fac[10000]);
for(int i=10000;i;--i) ifac[i-1]=mul(ifac[i],i);
for(int t=1;t<=T;++t)
{
int n=read(),m=read();static matrix A,E;static int r[51],s[51],f[51];
A=I(n),E=I(n);
for(int i=1,u,v;i<=m;++i) u=read(),v=read(),A[u][v]=A[v][u]=1;
work(A,f),memset(s,0,n<<2),s[0]=1;
for(int i=0;i<n;++i) r[i]=tr(E),E=E*A;
for(int i=0,x;i<=1e4;++i)
{
for(int j=0;j<n;++j) g[t][t][i]=inc(g[t][t][i],mul(s[j],r[j]));
x=s[n-1],memcpy(s+1,s,(n-1)<<2),s[0]=0;
if(x) for(int j=0;j<n;++j) s[j]=dec(s[j],mul(x,f[j]));
}
}
Q=read();
for(int i=1;i<=Q;++i) q[i]={read(),read(),read()},mx[q[i].l][q[i].r]=std::max(mx[q[i].l][q[i].r],q[i].k);
for(int i=1;i<=T;++i) for(int j=i+1;j<=T;++j) for(int k=0;k<=10000;++k) g[i][j][k]=mul(g[i][j-1][k],g[j][j][k]);
for(int i=1;i<=T;++i)
for(int j=i;j<=T;++j)
{
for(int k=0;k<=mx[i][j];++k) g[i][j][k]=mul(g[i][j][k],ifac[k]),f[i][j][k]=k&1? dec(0,ifac[k]):ifac[k];
MTT(g[i][j],f[i][j],f[i][j],mx[i][j]+1,mx[i][j]+1),f[i][j][0]=0;
for(int k=1;k<=mx[i][j];++k) f[i][j][k]=inc(f[i][j][k-1],mul(f[i][j][k],fac[k]));
}
for(int i=1;i<=Q;++i) printf("%d\n",f[q[i].l][q[i].r][q[i].k]);
}
