#loj3090 [BJOI2019] 勘破神机

Posted 2021-11-01 sweetphoenix

简单线性代数练习题
首先翻开具体数学生成函数一章,可以发现$F(n),G(n)$满足以下递推式

$$F(n)=F(n-1)+F(n-2),F(0)=1,F(1)=1$$
$$G(n)=4G(n-2)-G(n-4),G(2)=3,G(0)=1$$

我们发现$G$只有偶数项有值(证明的话可以把网格黑白染色什么的来证明)

那么我们可以设$T(n)=G(2n)$

那么$T$服从以下递推式

$$T(n)=4T(n-1)-T(n-2),T(1)=1,T(0)=1$$

那么我们发现$F$和$G$是两个二阶常系数线性递推数列

根据特征方程和生成函数那一套理论，我们可以知道$F,T$都有这样的通项公式

$$C×A^n+B×D^n$$

考虑题目让我们算什么

$$\sum_{i=1}^{n}{F(i) \choose k}$$

众所周知，组合数是个$O(k)$阶的多项式，那么我们现在的问题就变成了快速求

$$\sum_{i=1}^{n}F(i)^k$$

那么稍微画一下式子就是

$$\sum_{i=1}^{n}(CA^i+DB^i)^k$$

暴力二项式定理展开

$$\sum_{i=1}^{n}\sum_{j=0}^{k}C^{j}A^{ij}D^{k-j}B^{i(k-j)}$$

稍微交换一下求和号

$$\sum_{j=0}^{k}C^{j}D^{k-j}\sum_{i=1}^{n}(A^{j}B^{k-j})^i$$

后面显然是个等比数列求和，可以$O(logn)$的算出

那么这题就做完了，复杂度$O(Tk^2logn)$

注意当$m=3$的时候，你需要将输入的区间缩一下，因为询问的是$G$但是我们的做法只能处理$T$

另外你可能需要将输入的区间加加减减，因为我们的公式一般认为$F(0)=1$

下面是简单的线性代数内容，是关于如何解出$A,B,C,D$的

对于F显然是fib数列的通项公式，你应该熟练的背过它

$$A=\frac{\sqrt{5}+1}{2}$$

$$B=\frac{\sqrt{5}-1}{2}$$

$$C=\frac{\sqrt{5}}{5}$$

$$D=-\frac{\sqrt{5}}{5}$$

对于$T$稍微解一下特征方程可以得到

$$A=2+\sqrt{3}$$

$$B=2-\sqrt{3}$$

那么我们令$T(1)=1,T(2)=1$,就可以通过待定系数法解出$C,D$来

可能涉及到$\sqrt{3}$的方程组会有点复杂，推荐的解法是

令$x=C+D,y=\sqrt{3}(C-D)$

然后你会发现方程组干净了很多，从而我们可以解出x和y来，借助x和y就可以还原C,D了

解得

$$C=\frac{3+\sqrt{3}}{6},D=\frac{3-sqrt{3}}{6}$$

最后一个问题是5和3在膜998244353剩余系下统统没有二次剩余

把每个数字表示成$r+v\sqrt{k}$的形式即可实现模意义复数，然后就能计算了

等比数列求和涉及除法，使用这种复数实现除法的时候用分母有理化推个式子就行了

// luogu-judger-enable-o2
#include<cstdio>
#include<algorithm>
using namespace std;const int N=1e5+10;
typedef unsigned long long ll;const ll mod=998244353;
inline ll po(ll a,ll p){ll r=1;for(;p;p>>=1,a=a*a%mod)if(p&1)r=r*a%mod;return r;}
# define md(x) (x=(x>=mod)?x-mod:x)
ll sqr;
struct cmp
{
    ll r;ll v;
    cmp(ll R=0,ll V=0){r=R;v=V;}
    friend cmp operator +(cmp a,cmp b)
    {
        cmp c=cmp(a.r+b.r,a.v+b.v);
        md(c.r);md(c.v);return c;
    }
    friend cmp operator -(cmp a,cmp b)
    {
        cmp c=cmp(a.r+mod-b.r,a.v+mod-b.v);
        md(c.r);md(c.v);return c;
    }
    friend cmp operator *(cmp a,cmp b)
    {
        return cmp((a.r*b.r+sqr*a.v%mod*b.v)%mod,(a.r*b.v+b.r*a.v)%mod);
    }
    friend cmp operator /(cmp a,cmp b)
    {
        cmp c=a*cmp(b.r,mod-b.v);
        ll iv=((b.r*b.r+(mod-sqr)*b.v%mod*b.v))%mod;
        iv=po(iv,mod-2);
        (c.r*=iv)%=mod;(c.v*=iv)%=mod;
        return c;
    }
}A,B,C,D;int T;int m;
// fib(n) =(A^n-B^n)*C
inline cmp cpo(cmp a,ll p)
{
    cmp r=cmp(1,0);
    for(;p;p>>=1,a=a*a)if(p&1)r=r*a;return r;
}
ll ifac[N];ll c[1010][1010];ll sr1[1010][1010];
inline void pre()
{   
    ifac[0]=1;ifac[1]=1;
    for(int i=2;i<N;i++)
        ifac[i]=(mod-mod/i)*ifac[mod%i]%mod;
    for(int i=1;i<N;i++)
        (ifac[i]*=ifac[i-1])%=mod;
    for(int i=0;i<1006;i++)
    {
        c[i][0]=c[i][i]=1;
        for(int j=1;j<i;j++)
            c[i][j]=(c[i-1][j-1]+c[i-1][j])%mod;
    }
    sr1[0][0]=1;
    for(int i=0;i<1006;i++)
        for(int j=1;j<=i;j++)
            sr1[i][j]=(sr1[i-1][j-1]+(mod-i+1)*sr1[i-1][j])%mod;    
    if(m==2)
    {   
        ll iv2=(mod+1)/2;sqr=5;
        A=cmp(iv2,iv2);B=cmp(iv2,mod-iv2);
        C=cmp(0,po(5,mod-2));
        D=cmp(C.r,mod-C.v);
    }else 
    {
        ll iv6=po(6,mod-2);sqr=3;
        A=cmp(2,1);B=cmp(2,mod-1);
        C=cmp(3*iv6%mod,(mod-1)*iv6%mod);
        D=cmp(C.r,mod-C.v); 
    } 
}
//mar <cmp> st;mar <cmp> trs;
inline cmp csum(cmp q,ll n)
{
    if(q.r==1&&q.v==0)
        return (cmp){(n+1)%mod,0};
    cmp res=cpo(q,n+1)-cmp(1,0);
    res=res/(q-cmp(1,0));
//  printf("q=%lld,%lld,res=%lld,%lld\n",q.r,q.v,res.r,res.v);
    return res;
}
inline ll cfibk(ll n,int k)
{
    cmp res=cmp(0,0);
    for(int j=0;j<=k;j++)
    {
        cmp q=cpo(A,k-j)*cpo(B,j);
        cmp xs=(cmp){c[k][j],0}*cpo(C,k-j)*cpo(D,j);
        res=res+xs*csum(q,n);
    }
    return res.r;
}
namespace solver2
{
    inline void solve()
    {
        ll l;ll r;int k;
        scanf("%lld%lld%d",&l,&r,&k);
        l++;r++;ll res=0;
        for(int i=1;i<=k;i++)
            (res+=(cfibk(r,i)+mod-cfibk(l-1,i))*sr1[k][i])%=mod;
        (res*=ifac[k])%=mod;
        (res*=po((r-l+1)%mod,mod-2))%=mod;
        printf("%lld\n",res);
    }
}
namespace solver3
{
    inline void solve()
    {
        ll l;ll r;int k;
        scanf("%lld%lld%d",&l,&r,&k);
        ll len=r-l+1;
        if(l&1)l++;l/=2;
        if(r&1)r--;r/=2;
        l++;r++;
        ll res=0;
        for(int i=1;i<=k;i++)
            (res+=(cfibk(r,i)+mod-cfibk(l-1,i))*sr1[k][i])%=mod;
        (res*=ifac[k])%=mod;
        (res*=po(len%mod,mod-2))%=mod;
        printf("%lld\n",res);
    }
}
int main()
{
    scanf("%d%d",&T,&m);
    pre();
    if(m==2)
        for(int z=1;z<=T;z++)
            solver2::solve();
    else 
        for(int z=1;z<=T;z++)
            solver3::solve();
    return 0;
}