【XSY3399】枸杞（特征多项式，二项式反演）

显然有一个暴力的想法：我们先求出 \(F_{i,j}\) 表示第 \(i\) 张图内恰好走 \(j\) 步走回原地的方案数（相邻时刻位置可以相同），那么再用容斥得到答案。

分为两个 Part：一是加速 \(F_{i,j}\) 的求解过程，二是加速容斥和答案的求解过程。

先看第二个 Part：假设我们已经求出了所有的 \(F_{i,j}\)。

设 \(A_{j}=\prod_{i=l}^r F_{i,j}\)，根据二项式反演可以得到答案为：

\[\begin{aligned} ans&=\sum_{i=1}^{k}\sum_{j=0}^i(-1)^{j}\binom{i}{j}A_{i-j}\\ &=\sum_{i=1}^k i!\sum_{j=0}\frac{(-1)^{j}}{j!}\times \frac{A_{i-j}}{(i-j)!} \end{aligned} \]

那么我们可以先枚举 \(l,r\)（此时 \(A\) 数组是确定的），然后通过卷积求出 \(l,r\) 时的 \(ans\) 数组，询问的时候直接输出即可。

这部分时间复杂度为 \(O(T^2k\log k)-O(q)\)。

再看怎么求 \(F_{i,j}\)：

显然每一张图都要分开求，假设当前要求的图为第 \(i\) 张图。暴力的想法是设 \(g_{j,u}\) 表示恰好走 \(j\) 步走到 \(u\) 的方案数，然后枚举起点 \(st\)，初始化 \(g_{0,st}=1\)，转移 \(g_{j,u}=\sum_{(u,v)}g_{j-1,v}\)，每次将 \(g_{j,st}\) 统计进 \(F_{i,j}\) 中。时间复杂度 \(O(Tn^3k)\)。

转化到矩阵上：设邻接矩阵为 \(G\)，设 \(S_{st}\) 为除第 \(st\) 位为 \(1\) 以外的其他位全为 \(0\) 的列向量，\(g_{j,st}\) 即为 \(S_{st}\times G_j\) 第一列第 \(st\) 位。

那么不同起点就可以一起搞：把 \(S_i\) 并起来记为 \(S\)（发现这是单位矩阵，可以略去），我们要求的就是对于每一个 \(j\) 的 \(G^j\) 的对角线和（称为矩阵的迹）。

我们 \(O(n^3)\)（因为瓶颈不在这里，事实上可以暴力 \(O(n^4)\)：\(O(n^4)\) 暴力代入 \(\lambda =1,\cdots,n+1\) 行列式求值，再用带变元的拉格朗日插值 \(O(n^2)\) 插回 \(f(\lambda)\)）求出 \(G\) 的特征多项式 \(f(\lambda)=|\lambda I-G|\)，有 \(f(G)=0\)，于是 \(G^j=G^j\bmod f(G)\)。

现在假设我们已经求出了 \(x^{j-1}\bmod f(x)\)，那么我们可以 \(O(n)\) 求出 \(x^j\bmod f(x)\)，于是可以设 \(G^j=G^j\bmod f(G)=\sum\limits_{i=0}^{n-1}c_iG^i\)，于是得到 \(tr[G^j]=\sum\limits_{i=0}^{n-1}c_i\cdot tr[G^i]\)。

那么我们可以先 \(O(n^4)\) 把 \(G^0,G^1,\cdots,G^{n-1}\) 及它们的迹预处理出来， \(O(n)\) 代入即可得到 \(G^j\) 的迹。

时间复杂度 \(O(Tn^4+Tkn)\)。

总时间复杂度 \(O(Tn^4+Tkn+T^2k\log k+q)\)。

#include<bits/stdc++.h>

#define N 31
#define K 10010

using namespace std;

namespace modular
{
	const int mod=998244353;
	inline int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
	inline int dec(int x,int y){return x-y<0?x-y+mod:x-y;}
	inline int mul(int x,int y){return 1ll*x*y%mod;}
	inline void Add(int &x,int y){x=x+y>=mod?x+y-mod:x+y;}
	inline void Dec(int &x,int y){x=x-y<0?x-y+mod:x-y;}
	inline void Mul(int &x,int y){x=1ll*x*y%mod;}
}using namespace modular;

inline int poww(int a,int b)
{
	int ans=1;
	while(b)
	{
		if(b&1) ans=mul(ans,a);
		a=mul(a,a);
		b>>=1;
	}
	return ans;
}

inline int read()
{
	int x=0,f=1;
	char ch=getchar();
	while(ch<'0'||ch>'9')
	{
		if(ch=='-') f=-1;
		ch=getchar();
	}
	while(ch>='0'&&ch<='9')
	{
		x=(x<<1)+(x<<3)+(ch^'0');
		ch=getchar();
	}
	return x*f;
}

const int maxk=10000;

int T;
int fac[K],ifac[K];
int A[K<<2],B[K<<2];

namespace lpoly
{
	typedef vector<int> poly;
	poly pmul(poly a,int b)
	{
		for(int &x:a) x=mul(x,b);
		return a;
	}
	poly pmul(const poly &a,const poly &b)
	{
		const int sa=a.size(),sb=b.size();
		poly res(sa+sb-1);
		for(int i=0;i<sa;i++)
			for(int j=0;j<sb;j++)
				Add(res[i+j],mul(a[i],b[j]));
		return res;
	}
	poly padd(const poly &a,const poly &b)
	{
		const int sa=a.size(),sb=b.size();
		poly res(max(sa,sb));
		for(int i=0;i<min(sa,sb);i++)
			res[i]=add(a[i],b[i]);
		if(sa<sb) for(int i=sa;i<sb;i++) res[i]=b[i];
		else for(int i=sb;i<sa;i++) res[i]=a[i];
		return res;
	}
	poly pdiv(poly a,const poly &b)
	{
		const int sa=a.size(),sb=b.size();
		poly res(sa+1-sb);
		int tmp=poww(b[0],mod-2);
		for(int i=0;i<sa+1-sb;i++)
		{
			res[i]=mul(a[i],tmp);
			for(int j=0;j<sb;j++)
				Dec(a[i+j],mul(b[j],res[i]));
		}
		return res;
	}
}

namespace poly
{
	int rev[K<<2],w[16][K<<2][2];
	void init(int limit)
	{
		for(int bit=0,mid=1;mid<limit;bit++,mid<<=1)
		{
			int len=mid<<1;
			int gn=poww(3,(mod-1)/len);
			int ign=poww(gn,mod-2);
			int g=1,ig=1;
			for(int j=0;j<mid;j++,g=mul(g,gn),ig=mul(ig,ign))
				w[bit][j][0]=g,w[bit][j][1]=ig;
		}
	}
	void NTT(int *a,int limit,int opt)
	{
		opt=(opt<0);
		for(int i=0;i<limit;i++)
			rev[i]=(rev[i>>1]>>1)|((i&1)*(limit>>1));
		for(int i=0;i<limit;i++)
			if(i<rev[i]) swap(a[i],a[rev[i]]);
		for(int bit=0,mid=1;mid<limit;bit++,mid<<=1)
		{
			for(int i=0,len=mid<<1;i<limit;i+=len)
			{
				for(int j=0;j<mid;j++)
				{
					int x=a[i+j],y=mul(w[bit][j][opt],a[i+mid+j]);
					a[i+j]=add(x,y),a[i+mid+j]=dec(x,y);
				}
			}
		}
		if(opt)
		{
			int tmp=poww(limit,mod-2);
			for(int i=0;i<limit;i++)
				a[i]=mul(a[i],tmp);
		}
	}
}

struct Matrix
{
	int a[N][N];
	Matrix(){memset(a,0,sizeof(a));}
	void init(int n){for(int i=1;i<=n;i++)a[i][i]=1;}
};

Matrix Mmul(Matrix a,Matrix b,int n)
{
	Matrix c;
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
			for(int k=1;k<=n;k++)
				Add(c.a[i][j],mul(a.a[i][k],b.a[k][j]));
	return c;
}

int Det(Matrix M,int n)
{
	int ans=1;
	for(int i=1;i<=n;i++)
	{
		int p=i;
		for(int j=i+1;j<=n;j++)
			if(M.a[j][i]) p=j;
		if(p!=i)
		{
			swap(M.a[i],M.a[p]);
			ans=mod-ans;
		}
		int inv=poww(M.a[i][i],mod-2);
		for(int j=i+1;j<=n;j++)
		{
			int tmp=mul(M.a[j][i],inv);
			for(int k=i;k<=n;k++)
				Dec(M.a[j][k],mul(M.a[i][k],tmp));
		}
	}
	for(int i=1;i<=n;i++)
		ans=mul(ans,M.a[i][i]);
	return ans;
}

void Lagrange(int *y,lpoly::poly &res,int n)
{
	using lpoly::poly;
	using namespace lpoly;
	poly prod{1};
	for(int i=1;i<=n+1;i++)
		prod=pmul(prod,poly{mod-i,1});
	res.resize(0);
	for(int i=1;i<=n+1;i++)
	{
		int tmp=mul(y[i],mul(ifac[i-1],mul(((n+1-i)&1)?mod-1:1,ifac[n+1-i])));
		poly now=pdiv(prod,poly{mod-i,1});
		now=pmul(now,tmp);
		res=padd(res,now);
	}
}

struct Graph
{
	int n,m;
	int F[K];
	int tr[N];
	Matrix G[N];
	lpoly::poly f;
	void initG()
	{
		G[0].init(n);
		for(int i=1;i<=m;i++)
		{
			int u=read(),v=read();
			G[1].a[u][v]=G[1].a[v][u]=1;
		}
		for(int i=1;i<=n;i++) G[1].a[i][i]=1;
		for(int i=2;i<n;i++)
			G[i]=Mmul(G[i-1],G[1],n);
		for(int i=0;i<n;i++)
		{
			for(int j=1;j<=n;j++)
				Add(tr[i],G[i].a[j][j]);
			F[i]=tr[i];
		}
	}
	void getf()
	{
		int val[N];
		for(int i=1;i<=n+1;i++)
		{
			Matrix now;
			for(int j=1;j<=n;j++)
				for(int k=1;k<=n;k++)
					now.a[j][k]=dec(0,G[1].a[j][k]);
			for(int j=1;j<=n;j++)
				Add(now.a[j][j],i);
			val[i]=Det(now,n);
		}
		Lagrange(val,f,n);
		int inv=poww(f[n],mod-2);
		for(int i=0;i<n;i++) f[i]=dec(0,mul(f[i],inv));
		f.resize(n);
	}
	void work()
	{
		n=read(),m=read();
		initG();
		getf();
		using lpoly::poly;
		using namespace lpoly;
		poly now(n);
		now[n-1]=1;
		for(int i=n;i<=maxk;i++)
		{
			now=pmul(now,poly{0,1});
			now=padd(now,pmul(f,now[n]));
			now.resize(n);
			for(int j=0;j<n;j++) Add(F[i],mul(now[j],tr[j]));
		}
	}
}g[21];

int pf[K][21],ipf[K][21];
int ans[21][21][K];

int main()
{
//	freopen("ex1.in","r",stdin);
//	freopen("ex1.out","w",stdout);
	fac[0]=1;
	for(int i=1;i<=maxk;i++) fac[i]=mul(fac[i-1],i);
	ifac[maxk]=poww(fac[maxk],mod-2);
	for(int i=maxk;i>=1;i--) ifac[i-1]=mul(ifac[i],i);
	T=read();
	for(int i=1;i<=T;i++)
	{
		g[i].work();
	}
	for(int k=0;k<=maxk;k++)
	{
		pf[k][0]=ipf[k][0]=1;
		for(int i=1;i<=T;i++)
		{
			pf[k][i]=mul(pf[k][i-1],g[i].F[k]);
			ipf[k][i]=poww(pf[k][i],mod-2);
		}
	}
	int limit=1;
	while(limit<=(maxk<<1)) limit<<=1;
	poly::init(limit);
	for(int l=1;l<=T;l++)
	{
		for(int r=l;r<=T;r++)
		{
			for(int i=0;i<=maxk;i++) A[i]=mul(ifac[i],mul(pf[i][r],ipf[i][l-1]));
			for(int i=0;i<=maxk;i++) B[i]=mul((i&1)?mod-1:1,ifac[i]);
			poly::NTT(A,limit,1),poly::NTT(B,limit,1);
			for(int i=0;i<limit;i++) A[i]=mul(A[i],B[i]);
			poly::NTT(A,limit,-1);
			for(int i=1;i<=maxk;i++) ans[l][r][i]=add(ans[l][r][i-1],mul(fac[i],A[i]));
			for(int i=0;i<limit;i++) A[i]=B[i]=0;
		}
	}
	int q=read();
	while(q--)
	{
		int l=read(),r=read(),k=read();
		printf("%d\n",ans[l][r][k]);
	}
	return 0;
}
/*
2
3 3
1 2
2 3
1 3
3 2
1 2
2 3
3
1 1 4
2 2 4
1 2 2
*/