题解 P5320

洛谷题面传送门

神仙题（为什么就没能自己想出来呢/zk/zk）

这是我 AC 的第 (2 imes 10^3) 道题哦

首先考虑 (m=2) 的情况，我们首先可以想到一个非常 trivial 的 DP：(dp_i) 表示填好前 (i) 列的方案数，那么第 (i) 列显然有横着放和竖着放两种可能，方案数分别是 (dp_{i-2}) 和 (dp_{i-1})，因此我们有 (dp_i=dp_{i-2}+dp_{i-1})，边界条件 (dp_0=1)，显然这个递推式与斐波那契数列有非常紧密的联系，具体来说 (dp_i=f_{i+1})。因此我们可以将 (L,R) 分别加一，那么答案可以表示成：

[ans=sumlimits_{i=L}^Rdbinom{f_i}{k} ]

一脸不好用正常方法维护的样子。

不过注意到对于斐波那契数列而言我们有一个通项公式 (f_i=dfrac{1}{sqrt{5}}((dfrac{1+sqrt{5}}{2})^n-(dfrac{1-sqrt{5}}{2})^n))，因此我们考虑在这个通项公式上做点手脚。方便起见下文中记 (A=dfrac{1}{sqrt{5}},B=-dfrac{1}{sqrt{5}},X=dfrac{1+sqrt{5}}{2},Y=dfrac{1-sqrt{5}}{2})，那么 (f_i=AX^i+BY^i)。有人会说，这 (sqrt{5}) 无理数都出来了，况且 (sqrt{5}) 还在组合数上指标上，那不是更不好维护了吗？别急，(sqrt{5}) 这样的无理数确实和阶乘幂不太搭，不过别忘了我们还有个阶乘幂转下降幂的公式：(x^{underline{k}}=sumlimits_{i=0}^kegin{bmatrix}k\iend{bmatrix}(-1)^{k-i}x^i)，因此 (dbinom{f_i}{k}=dbinom{1}{k!}f_i^{underline{k}}=dfrac{1}{k!}sumlimits_{j=0}^kegin{bmatrix}k\jend{bmatrix}(-1)^{k-j}f_i^j)，因此我们有：

[egin{aligned} ans&=sumlimits_{i=L}^Rdbinom{f_i}{k}\ &=dfrac{1}{k!}sumlimits_{i=L}^Rsumlimits_{j=0}^kegin{bmatrix}k\jend{bmatrix}(-1)^{k-j}f_i^j\ &=dfrac{1}{k!}sumlimits_{i=L}^Rsumlimits_{j=0}^kegin{bmatrix}k\jend{bmatrix}(-1)^{k-j}(AX^i+BY^i)^j\ &=dfrac{1}{k!}sumlimits_{i=L}^Rsumlimits_{j=0}^kegin{bmatrix}k\jend{bmatrix}(-1)^{k-j}sumlimits_{p=0}^jdbinom{j}{p}(AX^i)^p(BY^i)^{j-p}\ &=dfrac{1}{k!}sumlimits_{j=0}^kegin{bmatrix}k\jend{bmatrix}(-1)^{k-j}sumlimits_{p=0}^jdbinom{j}{p}A^pB^{j-p}sumlimits_{i=L}^R(X^pY^{j-p})^i end{aligned} ]

前面都可以在 (mathcal O(k^2)) 时间内枚举得出，后面可以一波等比数列求和公式带走，时间复杂度 (k^2log R)。

有几个细节需要注意：

• 注意特判 (X^pY^{j-p}=1) 的情况
• 由于 (5) 在 (mod 998244353) 意义下无二次剩余（凉 心 出 题 人），因此我们需要像二次剩余那样将数域由 (mathbb{R}) 扩展到 (mathbb{Z}_5={a+bsqrt{5}|a,binmathbb{R}})（瞎起名字 ing），加减乘除照样定义，这样就可以避免这个问题了。

这个故事告诉我们，有时候遇到斐波那契数列有关的题，即使模数不那么友好也能套通项公式/ts

接下来考虑 (m=3)，首先 (n) 是奇数的情况答案肯定是 (0)，我们 duck 不必管他们。对于 (n) 是偶数的情况咱们还是先从最 naive 的 DP 入手，那么最后一列可以三个横着的 (1 imes 2) 多米诺骨牌，方案数 (dp_{i-2})，也可以是第一行放一个横着的骨牌，下面两列各放一个竖着的骨牌，当然也可以 reverse 一下，横着的骨牌放在第三行，方案数 (2dp_{i-2})，还可以第一行来一个横着的骨牌，紧接着下面跟着一个竖着的骨牌，紧接着左边又是 (2k(kin[1,infty)capmathbb{Z})) 个横着的骨牌，最后又一个竖着的，那么此时应从 (dp_{i-2k-2}) 转移过来，方案数 (2dp_{i-2k-2})，因此我们有

[dp_i=3dp_{i-1}+2sumlimits_{j=2k,kge 2}dp_{i-j} ]

前缀和搞搞可以变成：

[dp_i=4dp_{i-2}-dp_{i-4} ]

如果我们设 (g_n=dp_{2n})，那么上式可写作 (g_i=4g_{i-1}-g_{i-2})。

喜闻乐见的 (a_i=Aa_{i-1}+Ba_{i-2}) 的形式，列出特征根方程 (1=4x-x^2)，解出它的两个根为 (x_1=2+sqrt{3},x_2=2-sqrt{3})，然后套个待定系数法，设 (g_i=Ax_1^i+Bx_2^i)，那么：

[egin{cases} A+B=1\ (2+sqrt{3})A+(2-sqrt{3})B=3 end{cases} ]

解得

[egin{cases} A=dfrac{3+sqrt{3}}{6}\ B=dfrac{3-sqrt{3}}{6} end{cases} ]

然后就和 (m=2) 的情况一样了。

然鹅由于太久没有接触普通幂与阶乘幂之间的转换，一开始甚至把第一类斯特林数敲成了第二类斯特林数……这里稍微澄清一下：

• 用阶乘幂表示普通幂用第二类斯特林数，以等式 (x^k=sumlimits_{i=0}^kegin{Bmatrix}k\iend{Bmatrix}x^{underline{i}}) 为例，如果我们把等号左右两边叫做一区和二区，那么下划线在哪个区就用第几类斯特林数。
• 用大幂表示小幂需要加 ((-1)^{k-i})

const int MAXK=501;
const int MOD=998244353;
int qpow(int x,int e){
	int ret=1;
	for(;e;e>>=1,x=1ll*x*x%MOD) if(e&1) ret=1ll*ret*x%MOD;
	return ret;
}
int fac[MAXK+5],ifac[MAXK+5],s[MAXK+5][MAXK+5];
void init_fac(int n){
	for(int i=(fac[0]=ifac[0]=ifac[1]=1)+1;i<=n;i++) ifac[i]=1ll*ifac[MOD%i]*(MOD-MOD/i)%MOD;
	for(int i=1;i<=n;i++) fac[i]=1ll*fac[i-1]*i%MOD,ifac[i]=1ll*ifac[i-1]*ifac[i]%MOD;
	s[0][0]=1;for(int i=1;i<=n;i++) for(int j=1;j<=i;j++) s[i][j]=(s[i-1][j-1]+1ll*(i-1)*s[i-1][j])%MOD;
}
int sq;
struct comp{
	int x,y;
	comp(int _x=0,int _y=0):x(_x),y(_y){}
	comp operator +(const comp &rhs){return comp((x+rhs.x)%MOD,(y+rhs.y)%MOD);}
	comp operator -(const comp &rhs){return comp((x-rhs.x+MOD)%MOD,(y-rhs.y+MOD)%MOD);}
	comp operator *(const comp &rhs){return comp((1ll*x*rhs.x+1ll*y*rhs.y%MOD*sq)%MOD,(1ll*x*rhs.y+1ll*y*rhs.x)%MOD);}
	comp operator /(const comp &rhs){
		ll inv=qpow((1ll*rhs.x*rhs.x%MOD-1ll*sq*rhs.y%MOD*rhs.y%MOD+MOD)%MOD,MOD-2);
		return comp(1ll*(1ll*x*rhs.x%MOD-1ll*y*rhs.y%MOD*sq%MOD+MOD)*inv%MOD,1ll*(1ll*y*rhs.x%MOD-1ll*x*rhs.y%MOD+MOD)*inv%MOD);
	}
	bool operator ==(const comp &rhs){return (!(x^rhs.x)&&!(y^rhs.y));}
} A,B,X,Y;
int binom(int x,int y){return 1ll*fac[x]*ifac[y]%MOD*ifac[x-y]%MOD;}
comp qpow_c(comp x,ll e){
	comp res(1);
	for(;e;e>>=1,x=x*x) if(e&1) res=res*x;
	return res;
}
comp getsum(comp bs,ll t){
	if(t<0) return comp();
	if(bs==comp(1)) return comp((t+1)%MOD,0);
	return (qpow_c(bs,t+1)-comp(1))/(bs-comp(1));
}
comp getsum(comp bs,ll l,ll r){return getsum(bs,r)-getsum(bs,l-1);}
int solve(ll l,ll r,int k){
	if(l>r) return 0;comp res;
	for(int j=0;j<=k;j++){
		comp sum=comp();
		for(int p=0;p<=j;p++){
			comp mul(1);
			mul=mul*comp(binom(j,p))*qpow_c(A,p)*qpow_c(B,j-p);
			mul=mul*getsum(qpow_c(X,p)*qpow_c(Y,j-p),l,r);
			sum=sum+mul;
		} sum=sum*comp(s[k][j]);
//		printf("%d %d
",sum.x,sum.y);
		if((k-j)&1) res=res-sum;
		else res=res+sum;
	} return 1ll*res.x*ifac[k]%MOD;
}
int main(){
	int qu,opt;scanf("%d%d",&qu,&opt);init_fac(MAXK);
	if(opt==2){
		sq=5;A=comp(0,1)/comp(5);B=comp(0,MOD-1)/5;
		X=comp(1,1)/comp(2);Y=comp(1,MOD-1)/comp(2);
	} else {
		sq=3;A=comp(3,1)/comp(6);B=comp(3,MOD-1)/comp(6);
		X=comp(2,1);Y=comp(2,MOD-1);
	} while(qu--){
		ll l,r;int k;scanf("%lld%lld%d",&l,&r,&k);
		if(opt==2) printf("%d
",1ll*qpow((r-l+1)%MOD,MOD-2)*solve(l+1,r+1,k)%MOD);
		else printf("%d
",1ll*qpow((r-l+1)%MOD,MOD-2)*solve(l+1>>1,r>>1,k)%MOD);
	}
	return 0;
}