BJOI2019勘破神机（斯特林数+二项式定理+数学）

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=600,mod=998244353;
int c[N][N],s[N][N];
int qpow(int a,ll b)
{
    int ret=1;
    while(b)
    {
        if(b&1)ret=1ll*ret*a%mod;
        a=1ll*a*a%mod,b>>=1;
    }
    return ret;
}
struct num{
    int a,b,c;
    num operator+(num x){return (num){(a+x.a)%mod,(b+x.b)%mod,c};}
    num operator-(num x){return (num){(a-x.a+mod)%mod,(b-x.b+mod)%mod,c};}
    num operator*(num x)
    {return (num){(1ll*a*x.a+1ll*c*b%mod*x.b)%mod,(1ll*a*x.b+1ll*b*x.a)%mod,c};}
    num inv()
    {
        int f=qpow((1ll*a*a-1ll*b*b%mod*c%mod+mod)%mod,mod-2);
        return (num){1ll*a*f%mod,1ll*(mod-b)*f%mod,c};
    }
    bool operator==(num x){return a==x.a&&b==x.b&&c==x.c;}
    num operator/(num x){return (*this)*x.inv();}
};
num qpow(num a,ll b)
{
    num ret=(num){1,0,a.c};
    while(b)
    {
        if(b&1)ret=ret*a;
        a=a*a,b>>=1;
    }
    return ret;
}
int F(ll n,int k)
{
    if(!k)return n%mod;
    n++;
    num ans=(num){0,0,5},A=(num){1,1,5}/(num){2,0,5},B=(num){1,mod-1,5}/(num){2,0,5};
    num x=(num){1,0,5},y=(num){1,0,5};
    for(int j=1;j<=k;j++)y=y*B;
    for(int j=0;j<=k;j++)
    {
        int ret=c[k][j];
        if(k-j&1)ret=(mod-ret)%mod;
        num t=x*y;
        if(t==(num){1,0,5})ans=ans+(num){n%mod*ret%mod,0,5};
        else ans=ans+(num){ret,0,5}*(qpow(t,n+1)-t)/(t-(num){1,0,5});
        x=x*A,y=y/B;
    }
    for(int j=1;j<=k;j++)ans=ans/(num){0,1,5};
    return (ans.a+mod-1)%mod;
}
int calf(ll n,int k)
{
    if(!n)return 0;
    int ret=0,sum;
    for(int j=0;j<=k;j++)
    sum=1ll*F(n,j)*s[k][j]%mod,ret=(k-j&1)?(ret-sum+mod)%mod:(ret+sum)%mod;
    for(int i=1;i<=k;i++)ret=1ll*ret*qpow(i,mod-2)%mod;
    return ret;
}
int G(ll n,int k)
{
    if(!k)return n%mod;
    num ans=(num){0,0,3},x1=(num){2,1,3},x2=(num){2,mod-1,3};
    num A=(num){3,1,3},B=(num){3,mod-1,3},x=(num){1,0,3},y=x,a=x,b=x;
    for(int j=1;j<=k;j++)y=y*x2,b=b*B;
    for(int j=0;j<=k;j++)
    {
        int ret=c[k][j];
        num t=x*y;
        if(t==(num){1,0,3})ans=ans+a*b*(num){n%mod*ret%mod,0,3};
        else ans=ans+(num){ret,0,3}*a*b*(qpow(t,n+1)-t)/(t-(num){1,0,3});
        x=x*x1,y=y/x2,a=a*A;b=b/B;
    }
    for(int j=1;j<=k;j++)ans=ans/(num){6,0,3};
    return ans.a;
}
int calg(ll n,int k)
{
    if(!n)return 0;
    n/=2;
    int ret=0,sum;
    for(int j=0;j<=k;j++)
    sum=1ll*G(n,j)*s[k][j]%mod,ret=(k-j&1)?(ret-sum+mod)%mod:(ret+sum)%mod;
    for(int i=1;i<=k;i++)ret=1ll*ret*qpow(i,mod-2)%mod;
    return ret;
}
int main()
{
    for(int i=0;i<=590;i++)
    {
        c[i][0]=1;
        for(int j=1;j<=i;j++)c[i][j]=(c[i-1][j]+c[i-1][j-1])%mod; 
    }
    s[0][0]=1;
    for(int i=1;i<=590;i++)
    for(int j=1;j<=i;j++)
    s[i][j]=(s[i-1][j-1]+1ll*(i-1)*s[i-1][j])%mod;
    int T,m;scanf("%d%d",&T,&m);
    while(T--)
    {
        ll l,r;int k;cin>>l>>r>>k;
        int ans=m==2?(calf(r,k)-calf(l-1,k)+mod)%mod:(calg(r,k)-calg(l-1,k)+mod)%mod;
        ans=1ll*ans*qpow((r-l+1)%mod,mod-2)%mod;
        printf("%d\n",ans);
    }
}