「BJOI 2019」勘破神机

传送门

---

problem

经过了一个月的艰苦尝试，你的研究团队破译了 $“2”$ 型奥术宝石和 $“3”$ 型奥术宝石的内部能量结构。

这两类结构有着一定的相似性，它们的内部具有 $k$ 个反应核心，$“2”$ 型奥术宝石的每个核心都可以看成是一个 $2\times n$ 的网格，而 $“3”$ 型奥术宝石的每个核心都可以看成是一个 $3\times n$ 的网格。（注意奥术宝石的 $k$ 和 $n$ 可能不同）当神力反应进行时，每个核心自动填充满神力颗粒。

形式化地描述，每个神力颗粒可以看成是一个 $1\times 2$ 横置或竖置的方格，核心填满的定义为每个网格都恰好被一某个方格覆盖。若在两种填满反应核心的方案中存在一个方格放置的位置或方式不同，就认为方案不同。

如填满 $2\times 4$ 的网格有 $5$ 种不同的方案，填满 $3\times 2$ 的网格有 $3$ 种不同的方案。

如果奥术宝石的 $k$ 个核心的填充方式互不相同，它们就会组合出强大的咒术。黑袍想知道对于某个宝石一共有多少种不同的咒术（对于两种咒术组合，如果第一种咒术中每个核心 $a$ 的填充方式都可以找到第二种咒术的某个核心 $b$，使得 $a$ 和 $b$ 的填充方式完全相同，则认为这两种咒术组合相同）。

对于宽度为 $n$ 、反应核心个数为 $k$ 的 $“2”$ 型奥术宝石，设不同的咒术为 $F(n,k)$ ；对于宽度为 $n$ 、反应核心个数为 $k$ 的 $“3”$ 型奥术宝石，设不同的咒术为 $G(n,k)$ 。

地灾军团的科技水平还不能精准测量反应核心的长度 $n$，只能确定出核心长度的大致范围 $[l,r]$。你需要计算出反应核心长度在此区间内的平均咒术数，即

$$an{s}_2=\frac{1}{r-l+1}\sum _{n=l}^{r}F(n,k)$$

$$an{s}_3=\frac{1}{r-l+1}\sum _{n=l}^{r}G(n,k)$$

数据范围：$1\le l\le r\le 10^{18}$，$1\le k\le 501$。

---

solution

$m=2$ 的情况

设 $f_i$ 表示斐波那契数列（$f_1=1,f_2=1$），那么用 $1\times 2$ 的砖铺满 $2\times n$ 的格子的方案数是 $f_{n+1}$（为了方便，我们一开始就令 $l,r$ 加 $1$，这样答案就是 $f_n$），所以我们要求的实际上就是 ${\sum }_{n=l}^r\left(\genfrac{}{}{0ex}{}{f_n}{k}\right)$。

直接带有斐波那契的话好像不好做，我们考虑斐波那契数列的通项公式为：

$$f_n=\frac{1}{\sqrt{5}}{\left(\frac{1+\sqrt{5}}{2}\right)}^n-\frac{1}{\sqrt{5}}{\left(\frac{1-\sqrt{5}}{2}\right)}^n$$

那不妨设 $A=\frac{1}{\sqrt{5}},B=-\frac{1}{\sqrt{5}},x=\frac{1+\sqrt{5}}{2},y=\frac{1-\sqrt{5}}{2}$，那么 $f_n=A{x}^n+B{y}^n$。

这个先留着，待会用于把 $f_n$ 拆开。

在学习第一类斯特林数时，我们貌似学习过这个公式：

$$x^{\underline{n}}=\sum _{i=0}^n(-1{)}^{n-i}\left[\begin{array}{c}n\\ i\end{array}\right]x^i$$

其中 $\left[\begin{array}{c}n\\ i\end{array}\right]$ 为第一类斯特林数。证明就在这里看啦。

又根据 $\left(\genfrac{}{}{0ex}{}{n}{m}\right)=\frac{n^{\underline{m}}}{m!}$，我们继续化简：

$$\begin{array}{rl}\sum _{n=l}^r\left(\genfrac{}{}{0ex}{}{f_n}{k}\right)& =\frac{1}{k!}\sum _{n=l}^r\sum _{i=0}^k(-1{)}^{k-i}\left[\begin{array}{c}k\\ i\end{array}\right]f_n^i\\ & =\frac{1}{k!}\sum _{n=l}^r\sum _{i=0}^k(-1{)}^{k-i}\left[\begin{array}{c}k\\ i\end{array}\right](A{x}^n+B{y}^n{)}^i\\ & =\frac{1}{k!}\sum _{n=l}^r\sum _{i=0}^k(-1{)}^{k-i}\left[\begin{array}{c}k\\ i\end{array}\right]\sum _{j=0}^i(\genfrac{}{}{0ex}{}{i}{j})A^j{B}^{i-j}(x^j{y}^{i-j}{)}^n\\ & =\frac{1}{k!}\sum _{i=0}^k(-1{)}^{k-i}\left[\begin{array}{c}k\\ i\end{array}\right]\sum _{j=0}^i(\genfrac{}{}{0ex}{}{i}{j})A^j{B}^{i-j}\sum _{n=l}^r(x^j{y}^{i-j}{)}^n\end{array}$$

最后面的那个 ${\sum }_{n=l}^r(x^j{y}^{i-j}{)}^n$ 是一个等比数列求和，要注意有公比 $=1$ 的情况。

那我们就可以枚举 $i,j$ 然后 $O(k^2\log n)$ 算答案了。

还有一个细节是，就是 $\sqrt{5}$ 在模 $998244353$ 意义下是不存在的，那么我们类似于表示虚数的方法，把每个数表示成 $a+b\sqrt{5}$ 就可以了。

$m=3$ 的情况

以下内容部分来自 ViXBob 的博客，实在不想写写不动了。

先考虑求出递推式，比较显然的是当 $n$ 是奇数时答案是 $0$，所以设 $g_n$ 表示 $3\times 2n$ 的答案。

如下图：

枚举一下第二种情况跨越红线的数量就有递推式：

$$\begin{array}{rl}{g}_n& =3\times g_{n-1}+2\times \sum _{i=1}^{n-1}{g}_{n-1-i}\\ & =3\times g_{n-1}+2\times \sum _{i=0}^{n-2}{g}_i\end{array}$$

则有：

$$\{\begin{array}{l}{g}_n=3\times g_{n-1}+2\times {\sum }_{i=0}^{n-2}{g}_i\\ g_{n+1}=3\times g_n+2\times g_{n-1}+2\times {\sum }_{i=0}^{n-2}{g}_i\end{array}$$

差分一下：

$$g_{n+1}=4\times g_n-g_{n-1}$$

也就是说 $g_n$ 的递推式为：

$$\{\begin{array}{ll}{g}_n=1,& n=0\\ g_n=3,& n=1\\ g_n=4\times g_{n-1}-g_{n-2},& n>1\end{array}$$

用特征方程方程解得：

$$g_n=\frac{3+\sqrt{3}}{6}(2+\sqrt{3}{)}^n+\frac{3-\sqrt{3}}{6}(2-\sqrt{3}{)}^n$$

然后就是一模一样的求解方法了。

时间复杂度 $O(k^2\log n)$。

---

code

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=505,P=998244353,inv2=499122177,inv6=166374059;
int T,m,k,V,s[N][N],fac[N],ifac[N];
int add(int x,int y)  {return x+y>=P?x+y-P:x+y;}
int dec(int x,int y)  {return x-y< 0?x-y+P:x-y;}
int mul(int x,int y)  {return 1ll*x*y%P;}
int power(int a,int b){
	int ans=1;
	for(;b;b>>=1,a=mul(a,a))  if(b&1)  ans=mul(ans,a);
	return ans;
}
struct num{
	int x,y;
	num(int _x=0,int _y=0):x(_x),y(_y){}
	friend num operator+(const num &a,const num &b)  {return num(add(a.x,b.x),add(a.y,b.y));}
	friend num operator-(const num &a,const num &b)  {return num(dec(a.x,b.x),dec(a.y,b.y));}
	friend num operator-(const num &a,int b)  {return num(dec(a.x,b),a.y);}
	friend num operator*(const num &a,int b)  {return num(mul(a.x,b),mul(a.y,b));}
	friend num operator*(const num &a,const num &b){
		return num(add(mul(a.x,b.x),1ll*V*mul(a.y,b.y)%P),add(mul(a.x,b.y),mul(a.y,b.x)));
	}
	friend num Inv(const num &a){
		return num(a.x,P-a.y)*power(dec(mul(a.x,a.x),1ll*V*mul(a.y,a.y)%P),P-2);
	}
	friend num power(num a,ll b){
		num ans=num(1,0);
		for(;b;b>>=1,a=a*a)  if(b&1)  ans=ans*a;
		return ans;
	}
}A,B,X,Y;
void prework(){
	fac[0]=fac[1]=1;
	for(int i=2;i<N;++i)  fac[i]=mul(fac[i-1],i);
	ifac[N-1]=power(fac[N-1],P-2);
	for(int i=N-2;~i;--i)  ifac[i]=mul(ifac[i+1],i+1);
	s[0][0]=1;
	for(int i=1;i<N;++i)
		for(int j=1;j<N;++j)
			s[i][j]=add(s[i-1][j-1],mul(i-1,s[i-1][j]));
}
int C(int n,int m)  {return n<m?0:mul(fac[n],mul(ifac[m],ifac[n-m]));}
void init(){
	if(m==2){
		V=5;
		X=num(inv2,inv2),Y=num(inv2,P-inv2),A=Inv(num(0,1)),B=Inv(num(0,P-1));
	}
	else{
		V=3;
		X=num(2,1),Y=num(2,P-1),A=num(inv2,inv6),B=num(inv2,P-inv6);
	}
}
ll l,r,L,R;
int solve(){
	int ans=0;
	for(int i=0;i<=k;++i){
		int res=0,tmp=((k-i)&1)?(P-s[k][i]):s[k][i];
		for(int j=0;j<=i;++j){
			num val1=power(A,j)*power(B,i-j)*C(i,j),val2=power(X,j)*power(Y,i-j);
			val2=(val2.x==1&&val2.y==0)?num((r-l+1)%P,0):(power(val2,r+1)-power(val2,l))*Inv(val2-1);
			res=add(res,(val1*val2).x);
		}
		ans=add(ans,mul(res,tmp));
	}
	return mul(mul(ans,ifac[k]),power((R-L+1)%P,P-2));
}
int main(){
	prework();
	scanf("%d%d",&T,&m),init();
	while(T--){
		scanf("%lld%lld%d",&L,&R,&k);
		(m==2)?(l=L+1,r=R+1):(l=(L+1)/2,r=R/2);
		printf("%d\n",solve());
	}
	return 0;
}