• 「BJOI2019」勘破神机


    Description

    \(F_n\) 表示用 \(1\times 2\) 的骨牌填满 \(2\times n\) 的矩阵的方案数,\(G_n\) 表示用 \(1\times 2\) 的骨牌填满 \(3\times n\) 的矩阵的方案数。给出 \(l,r,k\),分别求出:

    \[\frac{1}{r-l+1}\sum_{n=l}^r\binom{F_n}{k} \]
    \[\frac{1}{r-l+1}\sum_{n=l}^r\binom{G_n}{k} \]

    \(998244353\) 取模。\(1\leq T\leq 5\)\(1\leq l\leq r\leq 10^{18}\)\(k\leq 501\)

    Solution

    m=2

    考虑第 \(i\) 列是竖着放还是横着放,\(F_i=F_{i-1}+F_{i-2}\)\(F_0=1\)。发现 \(F_i=f_{i+1}\)(其中 \(f_i\) 为斐波那契数列,\(f_1=f_2=1\))。为了方便,一开始就令 \(l,r\)\(1\),我们实际上要求 \(\sum_{n=l}^r\binom{f_n}{k}\)

    直接带 \(f\) 不太好做,考虑通项公式:

    \[f_n=\frac{1}{\sqrt 5}(\frac{1+\sqrt 5}{2})^n-\frac{1}{\sqrt 5}(\frac{1-\sqrt 5}{2})^n \]

    不妨设 \(A=\frac{1}{\sqrt 5}\)\(B=-\frac{1}{\sqrt 5}\)\(x=\frac{1+\sqrt 5}{2}\)\(y=\frac{1-\sqrt 5}{2}\),则 \(f_n=Ax^n+By^n\)。注意到 \(\sqrt 5\) 在组合数上指标上不太好维护。根据 \(\binom n m=\frac{n^{\underline m}}{m!}\),考虑下降幂转普通幂的公式 \(x^{\underline k}=\sum_{i=0}^k\begin{bmatrix}k\\i\end{bmatrix}(-1)^{k-i}x^i\)

    \[\begin{aligned} \sum_{n=l}^r\binom{f_n}{k} &=\frac{1}{k!}\sum_{n=l}^r f_n^{\underline k} =\frac{1}{k!}\sum_{n=l}^r\sum_{i=0}^k(-1)^{k-i}\begin{bmatrix}k\\i\end{bmatrix}f_n^i\\ &=\frac{1}{k!}\sum_{n=l}^r\sum_{i=0}^k(-1)^{k-i}\begin{bmatrix}k\\i\end{bmatrix}(Ax^n+By^n)^i\\ &=\frac{1}{k!}\sum_{n=l}^r\sum_{i=0}^k(-1)^{k-i}\begin{bmatrix}k\\i\end{bmatrix}\sum_{j=0}^i\binom i j A^jB^{i-j}(x^jy^{i-j})^n\\ &=\frac{1}{k!}\sum_{i=0}^k(-1)^{k-i}\begin{bmatrix}k\\i\end{bmatrix}\sum_{j=0}^i\binom i j A^jB^{i-j}\sum_{n=l}^r (x^jy^{i-j})^n \end{aligned} \]

    最后面的 \(\sum_{n=l}^r(x^jy^{i-j})^n\) 是个等比数列求和。这样就能 \(\mathcal O(k^2\log n)\) 计算了。注意:

    • 特判公比为 \(1\) 的情况。

    • \(\sqrt 5\)\(\bmod 998244353\) 意义下不存在,类似表示虚数的方法,将每个数表示成 \(A+B\sqrt 5\),加减乘除照样定义。最后算出来的答案一定满足 \(B=0\),不必担心。

      除法 \(\frac{a+b\sqrt 5}{c+d\sqrt 5}=\frac{(a+b\sqrt 5)(c-d\sqrt 5)}{(c+d\sqrt 5)(c-d\sqrt 5)}=\frac{ac-5bd}{c^2-5d^2}+\frac{bc-ad}{c^2-5d^2}\sqrt 5\)

    m=3

    先考虑求出递推式。\(G_i=3G_{i-2}+2\sum_{k\geq 1}G_{i-2k-2}\)

    image

    \(n\) 为奇数时答案为 \(0\),设 \(g_i=G_{2i}\),那么 \(g_i=3g_{i-1}+2\sum_{k=1}^{i-1}g_{i-k-1}=3g_{i-1}+2\sum_{k=0}^{i-2}g_k\)

    \[\begin{cases} g_i=3g_{i-1}+2\sum_{k=0}^{i-2}g_k\\ g_{i+1}=3g_i+2g_{i-1}+2\sum_{k=0}^{i-2}g_k \end{cases} \Rightarrow g_{i+1}=3g_i+(g_i-g_{i-1}) \Rightarrow g_{i+1}=4g_i-g_{i-1} \]

    特征方程为 \(x^2=4x-1\),解得 \(x_1=2+\sqrt 3,x_2=2-\sqrt 3\)。则 \(g_i=Ax_1^n+Bx_2^n\)。再根据 \(g_0=1,g_1=3\)

    \[\begin{cases} A+B=1\\ (2+\sqrt 3)A+(2-\sqrt 3)B=3 \end{cases} \Rightarrow \begin{cases} A=\frac{3+\sqrt 3}{6}\\ B=\frac{3-\sqrt 3}{6} \end{cases} \]

    然后就和 \(m=2\) 一样了。只是这里 \(l,r\) 没有 \(+1\) 而变成了 \(/2\)

    #include<bits/stdc++.h>
    #define ll long long
    using namespace std;
    const int N=510,mod=998244353;
    int t,m,w,c[N][N],s[N][N],inv[N],k,iv2=(mod+1)/2,tmp;
    ll l,r;
    int qpow(int x,int n){
        int ans=1;
        for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
        return ans;
    }
    struct num{
        int x,y;
        num operator+(num a){return {(x+a.x)%mod,(y+a.y)%mod};}
        num operator-(num a){return {(x-a.x+mod)%mod,(y-a.y+mod)%mod};}
        num operator*(int a){return {1ll*x*a%mod,1ll*y*a%mod};}
        num operator*(num a){return {(1ll*x*a.x%mod+1ll*y*a.y%mod*w%mod)%mod,(1ll*x*a.y%mod+1ll*y*a.x%mod)%mod};}
        num operator/(num a){
            int iv=qpow((1ll*a.x*a.x%mod-1ll*w*a.y%mod*a.y%mod+mod)%mod,mod-2);
            return {1ll*(1ll*x*a.x%mod-1ll*w*y%mod*a.y%mod+mod)%mod*iv%mod,1ll*(1ll*y*a.x%mod-1ll*x*a.y%mod+mod)%mod*iv%mod};
        }
    }a,b,x,y;
    num qpow(num x,ll n){
        num ans={1,0};
        for(;n;n>>=1,x=x*x) if(n&1) ans=ans*x;
        return ans;
    }
    int calc(){
        int ans=0;
        for(int i=0;i<=k;i++){
            int cnt=0;
            for(int j=0;j<=i;j++){
                num p=qpow(a,j)*qpow(b,i-j)*c[i][j],q=qpow(x,j)*qpow(y,i-j);
                q=q.x==1&&q.y==0?(num){(r-l+1)%mod,0}:(qpow(q,r+1)-qpow(q,l))/(q-(num){1,0});
                cnt=(cnt+(p*q).x)%mod;
            }
            ans=(ans+1ll*((k-i)&1?mod-1:1)*s[k][i]%mod*cnt%mod)%mod;
        }
        return 1ll*inv[k]%mod*ans%mod;
    }
    signed main(){
        inv[0]=inv[1]=1;
        for(int i=2;i<N;i++) inv[i]=1ll*inv[mod%i]*(mod-mod/i)%mod;
        c[0][0]=s[0][0]=1;
        for(int i=1;i<N;i++){
            c[i][0]=1,inv[i]=1ll*inv[i-1]*inv[i]%mod;
            for(int j=1;j<=i;j++){
                c[i][j]=(c[i-1][j]+c[i-1][j-1])%mod; 
                s[i][j]=(1ll*s[i-1][j]*(i-1)%mod+s[i-1][j-1])%mod;
            }
        }
        scanf("%d%d",&t,&m);
        if(m==2) w=5,a=(num){0,1}/(num){5,0},b=(num){0,-1}/(num){5,0},x={iv2,iv2},y={iv2,mod-iv2};
        else w=3,a=(num){3,1}/(num){6,0},b=(num){3,mod-1}/(num){6,0},x={2,1},y={2,mod-1};
        while(t--){
            scanf("%lld%lld%d",&l,&r,&k),tmp=qpow((r-l+1)%mod,mod-2);    //注意后面 r-l+1 改变了
            m==2?(l++,r++):(l=(l+1)/2,r/=2),printf("%lld\n",1ll*tmp*calc()%mod);
        }
        return 0;
    }
  • 相关阅读:
    如何保持页脚始终在页面底部
    CSS自适应宽度圆角按钮
    ACM1004
    java输出格式
    北大ACM1001题Exponentiation(求高精度幂)
    深入理解sizeof
    java之类BigDecimal
    ACM1003
    ACM1005
    C的输出格式printf
  • 原文地址:https://www.cnblogs.com/maoyiting/p/16415793.html
Copyright © 2020-2023  润新知