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

ooliver 发布于 4 小时前 48 次阅读 OI

AI 摘要

从广义斐波那契到动态 DP，矩阵快速幂能将递推加速至 O(log n)。本文将揭秘如何用状态压缩与矩阵转移攻克黑题 CF1152F2，甚至用广义矩阵乘法解决 Hamsters 问题。掌握这些矩阵技巧，让你的算法瞬间起飞！

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;
}