ZR 集训 Day11 – 矩阵、递推与动态 DP

ooliver 发布于 2026-07-27 391 次阅读 OI

AI 摘要

矩阵快速幂不止加速斐波那契——它还能碾过黑题CF1152F2，结合树剖实现动态DP，将O(n)优化至O(log²n)。本文带你一窥矩阵在递推与DP中的统治力。

ZR 集训 Day11 - 矩阵、递推与动态 DP

P1349 广义斐波那契数列

$$
\begin{bmatrix}
a_i & a_{i+1}
\end{bmatrix}

\times

\begin{bmatrix}
0 & q \\
1 & p
\end{bmatrix}

=

\begin{bmatrix}
a_{i+1} & q\times a_i+p\times a_{i+1}
\end{bmatrix}

=

\begin{bmatrix}
a_{i+1} & a_{i+2}
\end{bmatrix}
$$

使用矩阵快速幂加速即可。

代码：

C++
#include<bits/stdc++.h>
using namespace std;
#define int long long

int n,mod,p,q,a1,a2;

struct matrix{
    int c[3][3];
    matrix(){memset(c,0,sizeof c);}
};

matrix operator*(const matrix &x,const matrix &y){
    matrix tmp;
    for(int i=1;i<=2;i++)
        for(int j=1;j<=2;j++)
            for(int k=1;k<=2;k++)
                (tmp.c[i][j]+=x.c[i][k]*y.c[k][j]%mod)%=mod;
    return tmp;
}

matrix qpow(matrix A,int k){
    matrix res;
    for(int i=1;i<=2;i++) res.c[i][i]=1;
    while(k){
        if(k&1) res=res*A;
        A=A*A,k>>=1;
    }
    return res;
}

signed main(){
    matrix A,B;
    cin>>p>>q>>a1>>a2>>n>>mod;
    if(n==1) return cout<<a1%mod,0;
    if(n==2) return cout<<a2%mod,0;
    A.c[1][1]=0,A.c[1][2]=q,A.c[2][1]=1,A.c[2][2]=p;
    B.c[1][1]=a1,B.c[1][2]=a2;
    A=qpow(A,n-1);
    B=B*A;
    cout<<B.c[1][1];
    return 0;
}

CF1152F2 Neko Rules the Catniverse (Large Version)

也是做上黑题了。

这题跳来跳去是无序的，所以我们考虑枚举每一个星球，看是否能加入当前跳跃序列中。

我们用一个 $m$ 位二进制数 $S$ 维护一个滑动窗口，表示区间 $[i-m,i-1]$ 内选了哪些星球。

可以知道如果选择当前星球，可行的方案数应该是 $popcount(S)+1$，即可以加在整个跳跃序列的开头，或者每个星球的后面。

设 $dp_{i,j,S}$ 表示选到了第 $i$ 个星球，共 $j$ 次跳跃，滑动窗口状态为 $S$ 时的答案，转移方程分选 $i+1$ 和不选 $i+1$ 进行讨论，根据上面的方案书列出来即可，这里略。

会发现转移时和 $i$ 无关，并且转移只涉及乘法，而 $j$ 和 $S$ 的状态数又很小，所以考虑把 $j,S$ 压成一维，进行矩阵快速幂加速。

代码：

C++
#include<bits/stdc++.h>
using namespace std;
#define int long long

const int mod=1e9+7;
int n,k,m,sz;

struct matrix{
    int c[210][210];
    matrix(){memset(c,0,sizeof c);}
};

matrix operator*(const matrix &x,const matrix &y){
    matrix tmp;
    for(int i=1;i<=sz;i++)
        for(int k=1;k<=sz;k++)
            if(x.c[i][k])
                for(int j=1;j<=sz;j++)
                    (tmp.c[i][j]+=x.c[i][k]*y.c[k][j]%mod)%=mod;
    return tmp;
}

matrix qpow(matrix A,int k){
    matrix res;
    for(int i=1;i<=sz;i++) res.c[i][i]=1;
    while(k){
        if(k&1) res=res*A;
        A=A*A,k>>=1;
    }
    return res;
}

int get_id(int j,int S){
    return j*(1<<m)+S+1;
}

signed main(){
    cin>>n>>k>>m;
    sz=(k+1)*(1<<m);
    matrix A,B;
    int mask=(1<<m)-1;
    for(int j=0;j<=k;j++){
        for(int S=0;S<(1<<m);S++){
            int u=get_id(j,S);
            int S0=(S<<1)&mask;
            (A.c[u][get_id(j,S0)]+=1)%=mod;
            if(j<k){
                int S1=S0|1;
                int w=__builtin_popcount(S)+1;
                (A.c[u][get_id(j+1,S1)]+=w)%=mod;
            }
        }
    }
    B.c[1][get_id(0,0)]=1;
    A=qpow(A,n);
    B=B*A;
    int ans=0;
    for(int S=0;S<(1<<m);S++) (ans+=B.c[1][get_id(k,S)])%=mod;
    cout<<ans;
    return 0;
}

CF1117D Magic Gems

很容易推出转移方程：

$$
dp_i=dp_{i-1}+dp_{i-m}
$$

由于 $n$ 很大，考虑矩阵快速幂优化。

构造以下矩阵：

[dpidpi+1dpi+2…dpi+m−2dpi+m−1]×[000…01100…00010…00001…00…000…11]=[dpi+1dpi+2dpi+3…dpi+m−1dpi+m] \begin{bmatrix} dp_i & dp_{i+1} & dp_{i+2} & \dots & dp_{i+m-2} & dp_{i+m-1} \end{bmatrix} \times \begin{bmatrix} 0 & 0 & 0 &\dots & 0 & 1 \\ 1 & 0 & 0 &\dots & 0 & 0 \\ 0 & 1 & 0 &\dots & 0 & 0 \\ 0 & 0 & 1 &\dots & 0 & 0 \\ \dots \\ 0 & 0 & 0 &\dots & 1 & 1 \\ \end{bmatrix} = \begin{bmatrix} dp_{i+1} & dp_{i+2} & dp_{i+3} & \dots & dp_{i+m-1} & dp_{i+m} \end{bmatrix}

代码：

C++
#include<bits/stdc++.h>
using namespace std;
#define int long long

const int mod=1e9+7;
int n,k,m;

struct matrix{
    int c[110][110];
    matrix(){memset(c,0,sizeof c);}
};

matrix operator*(const matrix &x,const matrix &y){
    matrix tmp;
    for(int i=1;i<=m;i++)
        for(int k=1;k<=m;k++)
            if(x.c[i][k])
                for(int j=1;j<=m;j++)
                    (tmp.c[i][j]+=x.c[i][k]*y.c[k][j]%mod)%=mod;
    return tmp;
}

matrix qpow(matrix A,int k){
    matrix res;
    for(int i=1;i<=m;i++) res.c[i][i]=1;
    while(k){
        if(k&1) res=res*A;
        A=A*A,k>>=1;
    }
    return res;
}

signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    matrix A,B;
    cin>>n>>m;
    if(n-m+1<=0) return cout<<1,0;
    for(int i=1;i<=m;i++) B.c[1][i]=1,A.c[i+1][i]=1;
    A.c[1][m]=1,A.c[m][m]=1;
    A=qpow(A,n-m+1);
    B=B*A;
    cout<<B.c[1][m];
    return 0;
}

P1912 [NOI2009] 诗人小G

列出状态转移方程：

$$
\begin{aligned}
dp_i&=\min\{dp_j+w(j,i)\} \\
w(j,i)&=|sum_i-sum_j+i-j-1-L|^P
\end{aligned}
$$

不难发现 $w(i,j)$ 是满足四边形不等式的，所以这个 DP 满足决策单调性。

对于有决策单调性的 DP，我们通常使用单调队列+二分来进行优化。

我们用单调队列维护当前最优决策点，用一个数组 $a$ 表示第 $i$ 个点作为决策点时，将会在第 $a_{i-1}$ 个点时优于第 $i-1$ 个点。向队列添加新点时，我们将其与队尾进行比较，若新的 $a_{tail}$ 要小于 $a_{tail-1}$，就弹出，最后把当前决策点加入。

详细实现看代码：

C++
#include<bits/stdc++.h>
using namespace std;

#define int long long
typedef long double ld;
const int N=1e5+5;
int T,n,L,P;
int sum[N],opt[N],q[N],a[N],pm[N];
ld dp[N];
string str[N];

ld qpow(ld b,int p){
    ld a=1;
    while(p){
        if(p&1) a*=b;
        b*=b;
        p>>=1;
    }
    return a;
}

ld calc(int i,int j){
    return dp[j]+qpow(fabsl((ld)(sum[i]-sum[j]+i-j-1-L)),P);
}

int bound(int x,int y){
    int l=y,r=n+1;
    while(l<r){
        int mid=(l+r)>>1;
        if(calc(mid,x)>=calc(mid,y)) r=mid;
        else l=mid+1;
    }
    return l;
}

void solve(){
    cin>>n>>L>>P;
    for(int i=1;i<=n;i++){
        cin>>str[i];
        sum[i]=sum[i-1]+str[i].size();
    }
    int l=1,r=1;
    for(int i=1;i<=n;i++){
        while(l<r&&a[l]<=i) l++;
        dp[i]=calc(i,q[l]),opt[i]=q[l];
        while(l<r&&a[r-1]>=bound(q[r],i)) r--;
        a[r]=bound(q[r],i),q[++r]=i;
    }
    if(dp[n]>1e18) cout<<"Too hard to arrange\n";
    else{
        cout<<fixed<<setprecision(0)<<dp[n]<<"\n";
        int idx=0;
        for(int i=n;i;i=opt[i]) pm[++idx]=i;
        pm[++idx]=0;
        for(int i=idx;i>1;i--){
            int st=pm[i]+1,ed=pm[i-1];
            for(int j=st;j<ed;j++) cout<<str[j]<<" ";
            cout<<str[ed]<<"\n";
        }
    }
    cout<<"--------------------\n";
}

signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    cin>>T;
    while(T--) solve();
    return 0;
}

P3502 [POI 2010] CHO-Hamsters

用 $dis_{i,j}$ 表示第 $j$ 个字符串要接到第 $i$ 个字符串后面还需要几个字符，这里可以用字符串哈希先预处理一下。

设 $dp_{i,j}$ 表示当前串内共 $i$ 个名字，最后一个名字为 $j$，就得到了转移方程：

$$
dp_{i,j}=\min\{dp_{i-1,k}+dis_{k,j}\}(k\in [1,n])
$$

发现如果把 $\min$ 换成求和，把加法换成乘法，这就变成了矩阵乘法。

所以我们使用广义矩阵乘法进行矩阵快速幂优化，即用取最小值代替加法、用加法代替乘法。

代码：

C++
#include<bits/stdc++.h>
using namespace std;

#define int long long
typedef unsigned long long ull;
const int N=505,M=10005,base=101,inf=1e18;
int n,m,len[N];
string s[N];
ull mi[M],ha[M],hb[M];

struct mat{
    int a[N][N];
    mat(){for(int i=0;i<=n;i++) for(int j=0;j<=n;j++) a[i][j]=inf;}
};

mat operator*(mat x,mat y){
    mat z;
    for(int i=0;i<=n;i++) for(int k=0;k<=n;k++) if(x.a[i][k]!=inf)
        for(int j=0;j<=n;j++) if(x.a[i][k]+y.a[k][j]<z.a[i][j]) z.a[i][j]=x.a[i][k]+y.a[k][j];
    return z;
}

mat qpow(mat x,int k){
    mat res;
    for(int i=0;i<=n;i++) res.a[i][i]=0;
    while(k){
        if(k&1) res=res*x;
        x=x*x,k>>=1;
    }
    return res;
}

int calc(int x,int y,int same){
    for(int i=1;i<=len[x];i++) ha[i]=ha[i-1]*base+s[x][i-1];
    for(int i=1;i<=len[y];i++) hb[i]=hb[i-1]*base+s[y][i-1];
    for(int i=len[y]-same;i>=1;i--){
        if(len[x]<i) continue;
        ull h1=hb[i],h2=ha[len[x]]-ha[len[x]-i]*mi[i];
        if(h1==h2) return len[y]-i;
    }
    return len[y];
}

signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    cin>>n>>m;
    for(int i=1;i<=n;i++) cin>>s[i],len[i]=s[i].size();
    mi[0]=1;
    for(int i=1;i<M;i++) mi[i]=mi[i-1]*base;
    mat A;
    for(int i=0;i<=n;i++) for(int j=0;j<=n;j++){
        if(j==0) A.a[i][j]=inf;
        else if(i==0) A.a[i][j]=len[j];
        else A.a[i][j]=calc(i,j,i==j);
    }
    A=qpow(A,m);
    int ans=inf;
    for(int i=1;i<=n;i++) if(A.a[0][i]<ans) ans=A.a[0][i];
    cout<<ans<<"\n";
    return 0;
}

P4719 【模板】动态 DP

对于动态树上最大权独立集问题，我们首先考虑静态树 DP 的方程。设 $f_{u,0}$ 表示不选 $u$ 时 $u$ 子树的最大权独立集，$f_{u,1}$ 表示选择 $u$ 时 $u$ 子树的最大权独立集，易得转移方程为 $f_{u,0} = \sum_{v \in son(u)} \max(f_{v,0}, f_{v,1})$ 以及 $f_{u,1} = a_u + \sum_{v \in son(u)} f_{v,0}$。

为了在动态修改点权时高效维护 DP 值，我们需要借助树链剖分将子树贡献拆分为重儿子 $hs(u)$ 和轻儿子集合 $Light(u)$。我们定义辅助数组 $g_{u,0}$ 和 $g_{u,1}$ 分别存储轻儿子对 $u$ 的贡献：$g_{u,0} = \sum_{v \in Light(u)} \max(f_{v,0}, f_{v,1})$，$g_{u,1} = a_u + \sum_{v \in Light(u)} f_{v,0}$。此时 DP 转移方程可以重写为 $f_{u,0} = g_{u,0} + \max(f_{hs(u),0}, f_{hs(u),1})$ 和 $f_{u,1} = g_{u,1} + f_{hs(u),0}$。

观察上述方程，若我们重新定义广义矩阵乘法 $(A \otimes B)_{i,j} = \max_k (A_{i,k} + B_{k,j})$（即加法映射为 $\max$，乘法映射为 $+$），则转移方程可以完美写成矩阵形式：


$$
\begin{bmatrix} f_{u,0} \\ f_{u,1} \end{bmatrix} = \begin{bmatrix} g_{u,0} & g_{u,0} \\ g_{u,1} & -\infty \end{bmatrix} \otimes \begin{bmatrix} f_{hs(u),0} \\ f_{hs(u),1} \end{bmatrix}
$$


令 $M_u = \begin{bmatrix} g_{u,0} & g_{u,0} \ g_{u,1} & -\infty \end{bmatrix}$，由于树链剖分保证同一条重链上的节点在 DFS 序中连续，因此一条重链上的 DP 值对应着链上各节点转移矩阵区间乘积的连续复合。我们使用线段树维护重链上节点的矩阵乘积，线段树每个节点维护区间矩阵乘积 $tr[p] = tr[lc] \otimes tr[rc]$。

当单点 $x$ 的权值发生修改时，我们先更新 $g_{x,1}$ 和 $x$ 对应的矩阵 $M_x$，并在线段树中更新点 $x$。随后我们沿着重链向上跳（$x \leftarrow fa[top[x]]$），每次跳之前查询修改前该重链链顶的 DP 值，修改线段树后再查询新的链顶 DP 值，用两者的差值增量更新父节点 $fa[top[x]]$ 的 $g$ 数组，如此循环往复直至根节点。由于点到根节点最多经过 $O(\log n)$ 条重链，每次更新重链矩阵需 $O(\log n)$ 的线段树查询与修改，因此单次修改的时间复杂度为 $O(\log^2 n)$，整体时间复杂度为 $O((n+m) \log^2 n)$，空间复杂度为 $O(n)$。

C++
#include<bits/stdc++.h>
using namespace std;
#define lc p<<1
#define rc p<<1|1

const int N=100005,inf=1e9;

int n,m;
int a[N];
int dp[N],fa[N],sz[N],sn[N];
int id[N],wt[N],tp[N],ed[N],idx;
int f[N][2],g[N][2];
vector<int> t[N];

struct mat{
    int c[2][2];
    mat(){
        c[0][0]=c[0][1]=c[1][0]=c[1][1]=-inf;
    }
}vl[N],tr[N<<2];

mat operator*(const mat &x,const mat &y){
    mat tmp;
    for(int i=0;i<2;i++){
        for(int k=0;k<2;k++){
            for(int j=0;j<2;j++){
                tmp.c[i][j]=max(tmp.c[i][j],x.c[i][k]+y.c[k][j]);
            }
        }
    }
    return tmp;
}

void df1(int u,int p,int d){
    dp[u]=d,fa[u]=p,sz[u]=1;
    int mx=-1;
    for(int i=0;i<t[u].size();i++){
        int v=t[u][i];
        if(v==p) continue;
        df1(v,u,d+1);
        sz[u]+=sz[v];
        if(sz[v]>mx) mx=sz[v],sn[u]=v;
    }
}

void df2(int u,int tf){
    id[u]=++idx,wt[idx]=u,tp[u]=tf,ed[tf]=u;
    f[u][0]=0,f[u][1]=a[u];
    g[u][0]=0,g[u][1]=a[u];
    if(sn[u]){
        df2(sn[u],tf);
        f[u][0]+=max(f[sn[u]][0],f[sn[u]][1]);
        f[u][1]+=f[sn[u]][0];
    }
    for(int i=0;i<t[u].size();i++){
        int v=t[u][i];
        if(v==fa[u]||v==sn[u]) continue;
        df2(v,v);
        f[u][0]+=max(f[v][0],f[v][1]);
        f[u][1]+=f[v][0];
        g[u][0]+=max(f[v][0],f[v][1]);
        g[u][1]+=f[v][0];
    }
    vl[id[u]].c[0][0]=g[u][0],vl[id[u]].c[0][1]=g[u][0];
    vl[id[u]].c[1][0]=g[u][1],vl[id[u]].c[1][1]=-inf;
}

void pup(int p){
    tr[p]=tr[lc]*tr[rc];
}

void bld(int p,int l,int r){
    if(l==r){
        tr[p]=vl[l];
        return;
    }
    int md=l+r>>1;
    bld(lc,l,md);
    bld(rc,md+1,r);
    pup(p);
}

void chg(int p,int l,int r,int ps,const mat &k){
    if(l==r){
        tr[p]=k;
        return;
    }
    int md=l+r>>1;
    if(ps<=md) chg(lc,l,md,ps,k);
    else chg(rc,md+1,r,ps,k);
    pup(p);
}

mat que(int p,int l,int r,int ql,int qr){
    if(ql<=l&&r<=qr) return tr[p];
    int md=l+r>>1;
    if(qr<=md) return que(lc,l,md,ql,qr);
    if(ql>md) return que(rc,md+1,r,ql,qr);
    return que(lc,l,md,ql,qr)*que(rc,md+1,r,ql,qr);
}

void upd(int x,int y){
    g[x][1]+=y-a[x];
    a[x]=y;
    while(x){
        mat om=que(1,1,n,id[tp[x]],id[ed[tp[x]]]);
        int o0=max(om.c[0][0],om.c[0][1]);
        int o1=max(om.c[1][0],om.c[1][1]);

        vl[id[x]].c[0][0]=g[x][0],vl[id[x]].c[0][1]=g[x][0];
        vl[id[x]].c[1][0]=g[x][1],vl[id[x]].c[1][1]=-inf;
        chg(1,1,n,id[x],vl[id[x]]);

        mat nm=que(1,1,n,id[tp[x]],id[ed[tp[x]]]);
        int n0=max(nm.c[0][0],nm.c[0][1]);
        int n1=max(nm.c[1][0],nm.c[1][1]);

        x=fa[tp[x]];
        if(!x) break;

        g[x][0]+=max(n0,n1)-max(o0,o1);
        g[x][1]+=n0-o0;
    }
}

int main(){
    ios::sync_with_stdio(0);
    cin.tie(0);
    cin>>n>>m;
    for(int i=1;i<=n;i++) cin>>a[i];
    for(int i=1;i<n;i++){
        int u,v;
        cin>>u>>v;
        t[u].push_back(v);
        t[v].push_back(u);
    }
    df1(1,0,1);
    df2(1,1);
    bld(1,1,n);
    while(m--){
        int x,y;
        cin>>x>>y;
        upd(x,y);
        mat res=que(1,1,n,id[1],id[ed[1]]);
        cout<<max(max(res.c[0][0],res.c[0][1]),max(res.c[1][0],res.c[1][1]))<<"\n";
    }
    return 0;
}