题目描述
定义\(f_n\)为用\(1\times 2\)骨牌填满\(2\times n\)网格的方案数，\(g_n\)为填满\(3\times n\)网格的方案数。
求：
\[\frac{1}{r-l+1}\sum_{i=l}^rC_{f_i}^k/\frac{1}{r-l+1}\sum_{i=l}^rC_{g_i}^k\texttt{mod 998244353} \]\[l\leq r\leq 10^{18},k\leq 500 \]解题思路
\[C_{f_i}^k=\frac{f_i^{\underline k}}{k!} \]\[f_n^{\underline k}=\sum_{i=0}^k(-1)^{k-i}S_k^if_n^i \]其中 \(S_n^k\)代表第一类斯特林数。
\[\begin{aligned} \sum_{i=0}^nC_{f_i}^k&=\frac{1}{k!}\sum_{i=0}^n\sum_{j=0}^k(-1)^{k-j}S_k^jf_i^j\\ &=\frac{1}{k!}\sum_{j=0}^k(-1)^{k-j}S_k^j\sum_{i=0}^nf_i^j\\ \end{aligned} \]问题变为求解 \(\sum_{i=0}^nf_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{aligned} \sum_{i=0}^n(A\alpha^{i+1}+B\beta^{i+1})^k&=\sum_{i=0}^n\sum_{j=0}^kC_k^j(A\alpha^{i+1})^j(B\beta^{i+1})^{k-j}\\ &=\sum_{i=0}^n\sum_{j=0}^kC_k^jA^j\alpha^{(i+1)j}B^{k-j}\beta^{(i+1)(k-j)}\\ &=\sum_{i=0}^n\sum_{j=0}^kC_k^jA^jB^{k-j}(\alpha^j\beta^{k-j})^{i+1}\\ &=\sum_{j=0}^kC_k^jA^jB^{k-j}\sum_{i=0}^n(\alpha^j\beta^{k-j})^{i+1} \end{aligned} \]前面枚举后直接算，后边是等比数列求和可以\(O(1)\)算，需要特判公比为\(1\)的情况。
因为\(\sqrt 5\)不存在\(\texttt{mod 998244353}\)下的二次剩余，因此需要扩域，即把数字表示成\(a+b\sqrt 5\)的形式，类似复数的运算。
\(g_n\)也可以找出通项公式然后表示成和\(f_n\)相同的形式，可以类似算。
复杂度\(O(Tk^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;
}