[BJOI2019]勘破神机

题目描述

定义$f_n$为用$1\times 2$骨牌填满$2\times n$网格的方案数，$g_n$为填满$3\times n$网格的方案数。

求：

$$\frac{1}{r-l+1}\sum _{i=l}^r{C}_{f_i}^k/\frac{1}{r-l+1}\sum _{i=l}^r{C}_{g_i}^k\text{mod 998244353}$$

$$l\le r\le 10^{18},k\le 500$$

解题思路

$$C_{f_i}^k=\frac{f_i^{\underline{k}}}{k!}$$

$$f_n^{\underline{k}}=\sum _{i=0}^k(-1{)}^{k-i}{S}_k^i{f}_n^i$$

其中 $S_n^k$代表第一类斯特林数。
$$\begin{array}{rl}\sum _{i=0}^n{C}_{f_i}^k& =\frac{1}{k!}\sum _{i=0}^n\sum _{j=0}^k(-1{)}^{k-j}{S}_k^j{f}_i^j\\ & =\frac{1}{k!}\sum _{j=0}^k(-1{)}^{k-j}{S}_k^j\sum _{i=0}^n{f}_i^j\end{array}$$

问题变为求解 ${\sum }_{i=0}^n{f}_i^j$

众所周知 $f_n$ 是斐波那契第$n+1$项。

斐波那契通项公式：
$$f_n=\frac{1}{\sqrt{5}}((\frac{1+\sqrt{5}}{2}{)}^{n+1}-(\frac{1-\sqrt{5}}{2}{)}^{n+1})$$

可以化为$f_n=A{\alpha }^{n+1}+B{\beta }^{n+1}$

带入原式并展开：
$$\begin{array}{rl}\sum _{i=0}^n(A{\alpha }^{i+1}+B{\beta }^{i+1}{)}^k& =\sum _{i=0}^n\sum _{j=0}^k{C}_k^j(A{\alpha }^{i+1}{)}^j(B{\beta }^{i+1}{)}^{k-j}\\ & =\sum _{i=0}^n\sum _{j=0}^k{C}_k^j{A}^j{\alpha }^{(i+1)j}{B}^{k-j}{\beta }^{(i+1)(k-j)}\\ & =\sum _{i=0}^n\sum _{j=0}^k{C}_k^j{A}^j{B}^{k-j}({\alpha }^j{\beta }^{k-j}{)}^{i+1}\\ & =\sum _{j=0}^k{C}_k^j{A}^j{B}^{k-j}\sum _{i=0}^n({\alpha }^j{\beta }^{k-j}{)}^{i+1}\end{array}$$

前面枚举后直接算，后边是等比数列求和可以$O(1)$算，需要特判公比为$1$的情况。

因为$\sqrt{5}$不存在$\text{mod 998244353}$下的二次剩余，因此需要扩域，即把数字表示成$a+b\sqrt{5}$的形式，类似复数的运算。

$g_n$也可以找出通项公式然后表示成和$f_n$相同的形式，可以类似算。

复杂度$O(T{k}^2\log n)$

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int mod = 998244353;
inline LL Pow(LL a,LL b)
{
	a%=mod;
	LL res=1;
	while(b)
	{
		if(b&1)res=res*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	return res;
}
LL D;
struct Complex
{
	LL a,b;
	Complex(){a=b=0;}
	Complex(LL x,LL y){a=(x%mod+mod)%mod;b=(y%mod+mod)%mod;}
	Complex operator +(Complex x){return Complex(a+x.a,b+x.b);}
	Complex operator -(Complex x){return Complex(a-x.a,b-x.b);}
	Complex operator *(Complex x){return Complex(a*x.a+b*x.b%mod*D,a*x.b+b*x.a);}
};
Complex operator *(LL k,Complex x){return Complex(k*x.a,k*x.b);}
Complex inv(Complex x)
{
	LL P=(x.a*x.a%mod-x.b*x.b%mod*D%mod+mod)%mod;
	P=Pow(P,mod-2);
	return Complex(x.a*P%mod,mod-x.b*P%mod);
}
Complex Pow(Complex a,LL b)
{
	Complex res=Complex(1,0);
	a=Complex(a.a,a.b);
	while(b)
	{
		if(b&1)res=res*a;
		a=a*a;
		b>>=1;
	}
	return res;
}
const int N = 520;
LL S[N][N],C[N][N]; 
LL I2=Pow(2,mod-2),I6=Pow(6,mod-2);
int m;
Complex alpha,beta,A,B;
LL L,R,k;
LL pw(LL x){return (x&1)?mod-1:1;}
LL calc(LL n,LL k)
{
	LL coef=1;
	for(int i=1;i<=k;i++)coef=1ll*coef*i%mod;
	coef=Pow(coef,mod-2);
	Complex One=Complex(1,0),Zero=Complex(0,0);
	LL ans=0;
	for(int i=0;i<=k;i++)
	{
		LL P=1ll*S[k][i]*pw(k-i)%mod;
		Complex res=Zero;
		for(int j=0;j<=i;j++)
		{
			Complex u=Pow(alpha,i-j)*Pow(beta,j);
			Complex v=Pow(u,n+1);
			if(u.a==1&&u.b==0) v.a=(n+1)%mod,v.b=0;
			else u=inv(One-u),v=(One-v)*u;
			v=v*Pow(A,i-j)*Pow(B,j);
			v=C[i][j]*v;
			res=res+v;
		}
		ans=(ans+1ll*P*res.a%mod)%mod;
	}
	
	return 1ll*ans*coef%mod;
}
void solve()
{
	scanf("%lld %lld %lld",&L,&R,&k);
	LL coef=Pow(R-L+1,mod-2);
	if(m==2)R++;
	else L=(L-1)/2,R/=2;
	LL res=(calc(R,k)-calc(L,k)+mod)%mod*coef%mod;
	printf("%lld\n",res);
}
void put(Complex x)
{
	cout<<x.a<<' '<<x.b<<endl;
}
int main()
{
	C[0][0]=1;S[0][0]=1;
	for(int i=1;i<N;i++)
	{
		C[i][0]=1;
		for(int j=1;j<=i;j++)
		{
			C[i][j]=(C[i-1][j-1]+C[i-1][j])%mod;
			S[i][j]=(S[i-1][j-1]+1ll*(i-1)*S[i-1][j]%mod)%mod;
		}
	}
	int T;
	cin>>T>>m;
	if(m==2)
	{
		D=5;
		alpha=Complex(I2,I2);
		beta=Complex(I2,-I2);
		A=Complex(1,0)*inv(Complex(0,1));
		B=Complex(0,0)-A;
	}
	if(m==3)
	{
		D=3;
		alpha=Complex(2,-1);
		beta=Complex(2,1);
		A=Complex(I2,-I6);
		B=Complex(I2,I6);
	}
	while(T--)
	{
		solve();
	}
	return 0;
}