[BJOI2019]勘破神机（斯特林数，数论）

[BJOI2019]勘破神机（斯特林数，数论）

题面

洛谷

题解

先考虑(m=2)的情况。
显然方案数就是(f_i=f_{i-1}+f_{i-2})，即斐波那契数，虽然这里求出来是斐波那契的第(n+1)项，但是本质上没什么区别，就默认是斐波那契数列了。
斐波那契数列的特征根是(alpha=frac{1+sqrt 5}{2},eta=frac{1-sqrt 5}{2})，然后大力设一下通项是(f_n=Aalpha^n+Beta^n)，可以解出(f_n=frac{1}{sqrt 5}alpha^n-frac{1}{sqrt 5}eta^n)。
而要求的东西就是(displaystyle sum_{i=l}^r {f_ichoose k})，本质上就是要求(frac{1}{K!}f_i^{underline k})。
这个东西直接拿第一类斯特林数展开(displaystyle sum_jsum_{i=0}^k egin{bmatrix}k\ iend{bmatrix}f_j^k(-1)^{k-i})。
把斐波那契的通项直接带进去，得到：(displaystyle sum_{i=0}^kegin{bmatrix}k\ iend{bmatrix}(-1)^{k-i}sum_{j}(Aalpha^j+Balpha^j)^i)。
后面那个东西再拆一步就是(displaystyle sum_{j=0}^i{ichoose j}(Aalpha^n)^j(Beta^n)^{i-j})，
带回去交换求和顺序之后得到：(displaystyle sum_{k=0}^K(-1)^{K-k}egin{bmatrix}K\kend{bmatrix}sum_{i=0}^k{kchoose i}A^iB^{k-i}sum_{j=L}^R alpha^{ji}eta^{j(k-i)})。
后面一个部分显然是等比数列可以直接计算，这样子我们就把复杂度做到了(O(k^2))。

现在考虑(m=3)的情况，现在要做的显然就是把通项给求出来。
首先(n)为奇数的时候显然无解。
考虑为偶数的时候，令(g_i)表示第(i)个偶数的时候的答案，那么为了防止重复计算，我们需要找到一个放置方法满足在偶数位置一定不会恰好放满，发现这样子一定是在第一行竖着放一个横着放一个，然后后面每次都强制横着放跨过偶数位置，这里的贡献是(displaystyle sum_{j=0}^{i-2}2g_j)，然后(g_{i-1})的贡献比较特殊是(3)，因为这个位置可以不用强制跨过奇数位置，所以可以任意的放满(3*2)的格子，因此转移可以写成：(displaystyle g_i=g_{i-1}+2sum_{j=0}^{i-1}g_j)。
然后写两个出来差分：

[egin{cases}g_i=g_{i-1}+2sum_{j=0}^{i-1}g_j\ g_{i+1}=g_i+2sum_{j=0}^i g_jend{cases}]

推出(g_{i+1}-g_i=3g_i-g_{i-1})，推出递推式(g_{i+1}=4g_i-g_{i-1})。
然后特征方程：(x^2-4x+1=0)，两个解：(2pm sqrt 3)。
然后待定系数解出来通项是：(displaystyle g_n=frac{3+sqrt 3}{6}(2+sqrt 3)^n+frac{3-sqrt 3}{6}(2-sqrt 3)^n)。
然后就和前面的东西类似了。

#include<iostream>
#include<cstdio>
using namespace std;
#define ll long long
const int MOD=998244353;
#define inv2 499122177
#define inv5 598946612
#define inv6 166374059
inline ll read()
{
	ll x=0;bool t=false;char ch=getchar();
	while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
	if(ch=='-')t=true,ch=getchar();
	while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
	return t?-x:x;
}
int fpow(int a,ll b){int s=1;while(b){if(b&1)s=1ll*s*a%MOD;a=1ll*a*a%MOD;b>>=1;}return s;}
int S[505][505],C[505][505];ll V;
struct Num{int a,b;}z;
Num operator+(Num a,Num b){return (Num){(a.a+b.a)%MOD,(a.b+b.b)%MOD};}
Num operator-(Num a,Num b){return (Num){(a.a+MOD-b.a)%MOD,(a.b+MOD-b.b)%MOD};}
Num operator*(Num a,Num b){return (Num){(1ll*a.a*b.a+V*a.b*b.b)%MOD,(1ll*a.b*b.a+1ll*a.a*b.b)%MOD};}
Num operator*(Num a,int b){return (Num){1ll*a.a*b%MOD,1ll*a.b*b%MOD};}
Num fpow(Num a,ll b){Num s=(Num){1,0};while(b){if(b&1)s=s*a;a=a*a;b>>=1;}return s;}
int m,K;ll L,R;
Num Inv(Num a){return (Num){a.a,(MOD-a.b)%MOD}*fpow((1ll*a.a*a.a+MOD-V*a.b*a.b%MOD)%MOD,MOD-2);}
int main()
{
	int T=read();m=read();if(m==2)V=5;else V=3;
	while(T--)
	{
		L=read();R=read();K=read();
		S[0][0]=1;
		for(int i=1;i<=K;++i)
			for(int j=1;j<=i;++j)
				S[i][j]=(S[i-1][j-1]-1ll*S[i-1][j]*(i-1)%MOD+MOD)%MOD;
		for(int i=0;i<=K;++i)C[i][0]=1;
		for(int i=1;i<=K;++i)
			for(int j=1;j<=i;++j)C[i][j]=(C[i-1][j]+C[i-1][j-1])%MOD;
		Num alpha,beta,A,B;ll QwQ=R-L+1;
		if(m==2)
		{
			++L;++R;
			alpha=(Num){inv2,inv2},beta=(Num){inv2,MOD-inv2};
			A=(Num){0,inv5},B=(Num){0,MOD-inv5};
		}
		else
		{
			R=R>>1;L=(L+1)>>1;
			alpha=(Num){2,1},beta=(Num){2,MOD-1};
			A=(Num){inv2,inv6},B=(Num){inv2,MOD-inv6};
		}
		int ans=0;ll t=R-L;
		for(int i=0;i<=K;++i)
		{
			Num ret=z;
			for(int j=0;j<=i;++j)
			{
				Num s=fpow(A,j)*fpow(B,i-j)*C[i][j];
				Num p=fpow(fpow(alpha,L),j)*fpow(fpow(beta,L),(i-j));
				Num q=fpow(alpha,j)*fpow(beta,i-j);
				Num w=p*((fpow(q,t+1)-(Num){1,0})*Inv(q-(Num){1,0}));
				if(q.a==1&&q.b==0)w=p*((R-L+1)%MOD);
				ret=ret+s*w;
			}
			ans=(ans+1ll*ret.a*S[K][i])%MOD;
		}
		int Ans=1ll*ans*fpow(QwQ%MOD,MOD-2)%MOD;
		for(int i=1;i<=K;++i)Ans=1ll*Ans*fpow(i,MOD-2)%MOD;
		printf("%d
",Ans);
	}
	return 0;
}