(m=2)
此时答案为(frac{sumlimits_{i=l+1}^{r+1}{f_ichoose k}}{r-l+1}),其中(f_i)为第(i)个Fibonacci数。
也就是说我们现在要考虑如何求出(sumlimits_{i=l}^r{f_ichoose k})。
我们知道(f_n=frac1{sqrt5}(frac{1+sqrt 5}2)^n-frac1{sqrt5}(frac{1-sqrt5}2)^n),设(x=frac{1+sqrt5}2,y=frac{1-sqrt5}2,a=frac1{sqrt5},b=-frac1{sqrt5}),那么(f_n=ax^n+by^n)。
因此我们有
[egin{aligned}
&sumlimits_{n=l}^r{f_ichoose k}\
=&frac1{k!}sumlimits_{n=l}^rf_i^{underline k}\
=&frac1{k!}sumlimits_{n=l}^rsumlimits_{i=0}^k(-1)^{k-i}left[katop i
ight]f_n^i\
=&frac1{k!}sumlimits_{n=l}^rsumlimits_{i=0}^k(-1)^{k-i}left[katop i
ight]sumlimits_{j=0}^i{ichoose j}a^jb^{i-j}(x^jy^{i-j})^n\
&=frac1{k!}sumlimits_{i=0}^k(-1)^{k-i}left[katop i
ight]sumlimits_{j=0}^i{ichoose j}sumlimits_{n=l}^r(x^jy^{i-j})^n
end{aligned}
]
注意到(sqrt5mod998244353)并不存在,因此我们在(mathbb{F_p}(sqrt5))下计算即可。
(sumlimits_{n=l}^r(x^jy^{i-j})^n)部分可以直接运用等比数列求和公式求出,这样做的时间复杂度为(O(k^2log n))。
(m=3)
很显然(forall2
mid n,g_n=0),因此我们用(g_n)来表示(g_{2n})。
可以发现(g)的递推公式为(g_n=4g_{n-1}-g_{n-2}),解得(g_n=frac{3+sqrt3}6(2+sqrt3)^n+frac{3-sqrt3}6(2-sqrt3)^n)。
然后套用(m=2)时的做法即可。
利用Vieta定理+分治NTT可以做到(O(klog n+klog^2k))。
#include<cstdio>
#include<cstring>
#include<algorithm>
using i64=long long;
const int N=507,P=998244353;
int k;i64 C[N][N],S[N][N],m,l,r;
i64 add(i64 a,i64 b){return a+=b-P,a+(a>>63&P);}
i64 sub(i64 a,i64 b){return a-=b,a+(a>>63&P);}
i64 pow(i64 a,i64 b){i64 r=1;for(;b;b>>=1,a=a*a%P)if(b&1)r=a*r%P;return r;}
struct comp{i64 a,b;comp(i64 x=0,i64 y=0):a(x),b(y){}}pw1[N],pw2[N],pw3[N],pw4[N];
i64 norm(comp a){return sub(a.a*a.a%P,m*a.b*a.b%P);}
comp conj(comp a){return comp{a.a,sub(0,a.b)};}
int operator==(comp a,comp b){return a.a==b.a&&a.b==b.b;}
comp operator+(comp a,comp b){return comp{add(a.a,b.a),add(a.b,b.b)};}
comp operator-(comp a,comp b){return comp{sub(a.a,b.a),sub(a.b,b.b)};}
comp operator*(comp a,comp b){return comp{(a.a*b.a+m*a.b*b.b)%P,(a.a*b.b+a.b*b.a)%P};}
comp operator/(comp a,comp b){return a*conj(b)*pow(norm(b),P-2);}
comp operator^(comp a,i64 b){comp r(1);for(;b;b>>=1,a=a*a)if(b&1)r=a*r;return r;}
comp sum(comp z,i64 l,i64 r){return z==comp(1,0)? (r-l+1)%P:((z^(r+1))-(z^l))/(z-1);}
int main()
{
i64 fac=1,len;comp phi,x,ans;scanf("%*d%lld%lld%lld%d",&m,&l,&r,&k),len=(r-l+1)%P,pw1[0]=pw2[0]=pw3[0]=pw4[0]={1,0};
for(int i=2;i<=k;++i) fac=fac*i%P;
for(int i=0;i<=k;++i) for(int j=C[i][0]=1;j<=i;++j) C[i][j]=add(C[i-1][j-1],C[i-1][j]);
for(int i=S[0][0]=1;i<=k;++i) for(int j=1;j<=i;++j) S[i][j]=sub(S[i-1][j-1],S[i-1][j]*(i-1)%P);
if(m==2) ++l,++r,m=5,phi=comp(0,1)/5,x=comp(1,1)/2; else if(l=(l+1)/2,r/=2,l>r) return puts("0"),0; else phi=comp(3,1)/6,x=comp(2,1);
for(int i=1;i<=k;++i) pw1[i]=pw1[i-1]*phi,pw2[i]=pw2[i-1]*conj(phi),pw3[i]=pw3[i-1]*x,pw4[i]=pw4[i-1]*conj(x);
for(int i=0;i<=k;++i) for(int j=0;j<=i;++j) ans=ans+(S[k][i]*C[i][j])%P*pw1[j]*pw2[i-j]*sum(pw3[j]*pw4[i-j],l,r);
printf("%lld",ans.a*pow(len*fac%P,P-2)%P);
}