[BJOI2019]勘破神機(斯特林數,數論)

[BJOI2019]勘破神機(斯特林數,數論)

題面

洛谷ios

題解

先考慮\(m=2\)的狀況。
顯然方案數就是\(f_i=f_{i-1}+f_{i-2}\),即斐波那契數,雖然這裏求出來是斐波那契的第\(n+1\)項,可是本質上沒什麼區別,就默認是斐波那契數列了。
斐波那契數列的特徵根是\(\alpha=\frac{1+\sqrt 5}{2},\beta=\frac{1-\sqrt 5}{2}\),而後大力設一下通項是\(f_n=A\alpha^n+B\beta^n\),能夠解出\(f_n=\frac{1}{\sqrt 5}\alpha^n-\frac{1}{\sqrt 5}\beta^n\)
而要求的東西就是\(\displaystyle \sum_{i=l}^r {f_i\choose k}\),本質上就是要求\(\frac{1}{K!}f_i^{\underline k}\)
這個東西直接拿第一類斯特林數展開\(\displaystyle \sum_j\sum_{i=0}^k \begin{bmatrix}k\\ i\end{bmatrix}f_j^k(-1)^{k-i}\)
把斐波那契的通項直接帶進去,獲得:\(\displaystyle \sum_{i=0}^k\begin{bmatrix}k\\ i\end{bmatrix}(-1)^{k-i}\sum_{j}(A\alpha^j+B\alpha^j)^i\)
後面那個東西再拆一步就是\(\displaystyle \sum_{j=0}^i{i\choose j}(A\alpha^n)^j(B\beta^n)^{i-j}\)
帶回去交換求和順序以後獲得:\(\displaystyle \sum_{k=0}^K(-1)^{K-k}\begin{bmatrix}K\\k\end{bmatrix}\sum_{i=0}^k{k\choose i}A^iB^{k-i}\sum_{j=L}^R \alpha^{ji}\beta^{j(k-i)}\)
後面一個部分顯然是等比數列能夠直接計算,這樣子咱們就把複雜度作到了\(O(k^2)\)spa

如今考慮\(m=3\)的狀況,如今要作的顯然就是把通項給求出來。
首先\(n\)爲奇數的時候顯然無解。
考慮爲偶數的時候,令\(g_i\)表示第\(i\)個偶數的時候的答案,那麼爲了防止重複計算,咱們須要找到一個放置方法知足在偶數位置必定不會剛好放滿,發現這樣子必定是在第一行豎着放一個橫着放一個,而後後面每次都強制橫着放跨過偶數位置,這裏的貢獻是\(\displaystyle \sum_{j=0}^{i-2}2g_j\),而後\(g_{i-1}\)的貢獻比較特殊是\(3\),由於這個位置能夠不用強制跨過奇數位置,因此能夠任意的放滿\(3*2\)的格子,所以轉移能夠寫成:\(\displaystyle g_i=g_{i-1}+2\sum_{j=0}^{i-1}g_j\)
而後寫兩個出來差分:
\[\begin{cases}g_i=g_{i-1}+2\sum_{j=0}^{i-1}g_j\\ g_{i+1}=g_i+2\sum_{j=0}^i g_j\end{cases}\]
推出\(g_{i+1}-g_i=3g_i-g_{i-1}\),推出遞推式\(g_{i+1}=4g_i-g_{i-1}\)
而後特徵方程:\(x^2-4x+1=0\),兩個解:\(2\pm \sqrt 3\)
而後待定係數解出來通項是:\(\displaystyle g_n=\frac{3+\sqrt 3}{6}(2+\sqrt 3)^n+\frac{3-\sqrt 3}{6}(2-\sqrt 3)^n\)
而後就和前面的東西相似了。code

#include<iostream>
#include<cstdio>
using namespace std;
#define ll long long
const int MOD=998244353;
#define inv2 499122177
#define inv5 598946612
#define inv6 166374059
inline ll read()
{
    ll x=0;bool t=false;char ch=getchar();
    while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
    if(ch=='-')t=true,ch=getchar();
    while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
    return t?-x:x;
}
int fpow(int a,ll b){int s=1;while(b){if(b&1)s=1ll*s*a%MOD;a=1ll*a*a%MOD;b>>=1;}return s;}
int S[505][505],C[505][505];ll V;
struct Num{int a,b;}z;
Num operator+(Num a,Num b){return (Num){(a.a+b.a)%MOD,(a.b+b.b)%MOD};}
Num operator-(Num a,Num b){return (Num){(a.a+MOD-b.a)%MOD,(a.b+MOD-b.b)%MOD};}
Num operator*(Num a,Num b){return (Num){(1ll*a.a*b.a+V*a.b*b.b)%MOD,(1ll*a.b*b.a+1ll*a.a*b.b)%MOD};}
Num operator*(Num a,int b){return (Num){1ll*a.a*b%MOD,1ll*a.b*b%MOD};}
Num fpow(Num a,ll b){Num s=(Num){1,0};while(b){if(b&1)s=s*a;a=a*a;b>>=1;}return s;}
int m,K;ll L,R;
Num Inv(Num a){return (Num){a.a,(MOD-a.b)%MOD}*fpow((1ll*a.a*a.a+MOD-V*a.b*a.b%MOD)%MOD,MOD-2);}
int main()
{
    int T=read();m=read();if(m==2)V=5;else V=3;
    while(T--)
    {
        L=read();R=read();K=read();
        S[0][0]=1;
        for(int i=1;i<=K;++i)
            for(int j=1;j<=i;++j)
                S[i][j]=(S[i-1][j-1]-1ll*S[i-1][j]*(i-1)%MOD+MOD)%MOD;
        for(int i=0;i<=K;++i)C[i][0]=1;
        for(int i=1;i<=K;++i)
            for(int j=1;j<=i;++j)C[i][j]=(C[i-1][j]+C[i-1][j-1])%MOD;
        Num alpha,beta,A,B;ll QwQ=R-L+1;
        if(m==2)
        {
            ++L;++R;
            alpha=(Num){inv2,inv2},beta=(Num){inv2,MOD-inv2};
            A=(Num){0,inv5},B=(Num){0,MOD-inv5};
        }
        else
        {
            R=R>>1;L=(L+1)>>1;
            alpha=(Num){2,1},beta=(Num){2,MOD-1};
            A=(Num){inv2,inv6},B=(Num){inv2,MOD-inv6};
        }
        int ans=0;ll t=R-L;
        for(int i=0;i<=K;++i)
        {
            Num ret=z;
            for(int j=0;j<=i;++j)
            {
                Num s=fpow(A,j)*fpow(B,i-j)*C[i][j];
                Num p=fpow(fpow(alpha,L),j)*fpow(fpow(beta,L),(i-j));
                Num q=fpow(alpha,j)*fpow(beta,i-j);
                Num w=p*((fpow(q,t+1)-(Num){1,0})*Inv(q-(Num){1,0}));
                if(q.a==1&&q.b==0)w=p*((R-L+1)%MOD);
                ret=ret+s*w;
            }
            ans=(ans+1ll*ret.a*S[K][i])%MOD;
        }
        int Ans=1ll*ans*fpow(QwQ%MOD,MOD-2)%MOD;
        for(int i=1;i<=K;++i)Ans=1ll*Ans*fpow(i,MOD-2)%MOD;
        printf("%d\n",Ans);
    }
    return 0;
}
相關文章
相關標籤/搜索