luogu P4593 [TJOI2018]教科书般的亵渎

Posted 2020-11-25 smyjr

luogu

这题怎么没人用矩乘啊

首先可以发现,这题的 (k) ,也就是亵渎使用次数为 (m+1) ，然后给出的 (a_i) 又会把 ([1,n]) 划分成至多 (m+1) 个连续段。所以对每次亵渎，一个段 ([l,r]) 给答案加上 (sum_{i=l}^r i^k) ，并且所有连续段整体 (-1) ；此外如果某次亵渎后左端点最小的段的左端点为 (0) ，记它的右端点为 (r_1) ，那么会删掉这个段，并且其他段整体 (-r_1)。

那么我们模拟上述过程，然后每次对段操作前统计 (sum_{i=l}^r i^k) 即可。考虑矩乘求这个东西，先把他拆成 (sum_{i=1}^r i^k-sum_{i=1}^{l-1} i^k) ，然后我们构造行向量 (vec{a}={x^0,x^1,x^2...x^{m+1},s}) ，即维护当前 (i) 的 (0) 到 (m+1) 次幂和答案，构造转移矩阵 (M) 就考虑 ((i+1)^k=sum_{j=0}^{k}inom{k}{j}i^j) ，所以可以得到

• (forall 0le i,jle m+1,M_{i,j}=inom{j}{i})
• (M_{m+2,m+2}=1)
• (forall 0le ile m+1,M_{i,m+2}=inom{m+1}{i})

(vec{a}M^{n}) 的第 (m+2) 项就是 (sum i^k) 。不过注意到要做 (m^2) 次这样的求值，所以复杂度为 (O(Tm^5logn))。

现在考虑倒着模拟亵渎的过程，可以发现相当于是要维护一个行向量集合，每次操作先给所有行向量乘上 (M) ，然后加入一个初始行向量 ({1,0,0...}) ，再给所有行向量乘 (M^b) (其中 (b) 这次亵渎删掉的段的右端点值)接着加入一个行向量 ({-1,0,0...})，最后把所有行向量的第 (m+2) 项加入答案。由于矩阵乘法具有分配律，所以可以直接维护所有行向量的和，然后进行一些矩阵加法和乘法即可做到 (O(Tm^4logn)) 。注意到我们矩乘做的是行向量乘矩阵，所以单次矩乘可以 (O(m^2)) ，然后我们再把转移矩阵的 (2^k) 次幂预处理出来，就可以做到 (O(Tm^3logn))。

大 力 出 奇 迹

#include<bits/stdc++.h>
#define LL long long

using namespace std;
const int N=53,mod=1e9+7;
const LL inf=8e18;
LL rd()
{
    LL x=0,w=1;char ch=0;
    while(ch<'0'||ch>'9'){if(ch=='-') w=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=x*10+(ch^48);ch=getchar();}
    return x*w; 
}
LL n,b[N],hp[N][2];
int m=51,t,p;
struct matrix
{
    int a[N][N];
    matrix(){memset(a,0,sizeof(a));}
    void clr(){memset(a,0,sizeof(a));}
    matrix operator + (const matrix &bb) const
    {
        matrix an;
        for(int i=0;i<=m+1;++i)
            for(int j=0;j<=m+1;++j)
            an.a[i][j]=(a[i][j]+bb.a[i][j])%mod;
        return an;
    }
    matrix operator * (const matrix &bb) const
    {
        matrix an;
        for(int i=0;i<=m+1;++i)
            for(int j=0;j<=m+1;++j)
            {
                LL nw=0;
                for(int k=0;k<=m+1;++k)
                {
                    nw+=1ll*a[i][k]*bb.a[k][j];
                    if(nw>inf) nw%=mod;
                }
                an.a[i][j]=nw%mod;
            }
        return an;
    }
    matrix operator & (const matrix &bb) const
    {
        matrix an;
        for(int j=0;j<=m+1;++j)
        {
            LL nw=0;
            for(int k=0;k<=m+1;++k)
            {
                nw+=1ll*a[0][k]*bb.a[k][j];
                if(nw>inf) nw%=mod;
            }
            an.a[0][j]=nw%mod;
        }
        return an;
    }
}aa,m1,bb[N][N];
int c[N][N];
void inii(int h)
{
    for(int i=0;i<=h;++i)
    for(int j=0;j<=h;++j)
        bb[h][0].a[j][i]=c[i][j];
    bb[h][0].a[h+1][h+1]=1;
    for(int j=0;j<=h;++j)
    bb[h][0].a[j][h+1]=c[h][j];
    for(int i=1;i<N;++i) bb[h][i]=bb[h][i-1]*bb[h][i-1];
}

int main()
{
    for(int i=0;i<N;++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;
    }
    aa.a[0][0]=1;
    int T=rd();
    while(T--)
    {
        n=rd(),m=rd();
        b[0]=0;
        for(int i=1;i<=m;++i) b[i]=rd();
        sort(b+1,b+m+1),m=unique(b+1,b+m+1)-b;
        while(m>1&&b[m-1]==n) --m,--n;
        if(!bb[m][0].a[0][0]) inii(m);
        b[m]=n+1,t=0;
        for(int i=1;i<=m;++i)
        {
            if(b[i]==b[i-1]+1) continue;
            hp[++t][0]=b[i-1]+1,hp[t][1]=b[i]-b[i-1]-1;
        }
        p=0;
        while(t)
        {
            b[++p]=0;
            for(int i=1;i<=t;++i) --hp[i][0];
            if(!hp[1][0])
            {
                b[p]=hp[1][1];
                for(int i=2;i<=t;++i) hp[i][0]-=hp[1][1];
                for(int i=2;i<=t;++i) hp[i-1][0]=hp[i][0],hp[i-1][1]=hp[i][1];
                --t;
            }
        }
        m1.clr();
        int ans=0;
        while(p)
        {
            m1=m1*bb[m][0];
            if(b[p])
            {
                m1.a[0][0]=(m1.a[0][0]+1)%mod;
                for(int i=0;b[p];++i,b[p]>>=1)
                    if(b[p]&1) m1=m1&bb[m][i];
                m1.a[0][0]=(m1.a[0][0]-1+mod)%mod;
            }
            ans=(ans+m1.a[0][m+1])%mod;
            --p;
        }
        printf("%d
",ans);
    }
    return 0;
}