好耶!来学新算法了(最近停课了就有时间学算法啦
因为 CSP-S 考了个什么 ddp 然后我不会
(CSP炸了)
学了一个晚上加一个上午才学会(我太菜了)
嗯嗯其实说起来是个很简单的东西.
前置知识:
树上dp,线段树,树链剖分,矩阵优化线性 dp
例题:
给你一棵 \(n\) 个节点的树,点有点权,有 \(m\) 次单点修改操作,每次求出一个独立点集(集合内任意两点没有边相连)的最大权值和.( \(1\le n,m \le 10^5\) )
\(O(nm)\) 暴力
首先 \(O(nm)\) 暴力很简单,直接对每次操作做一遍 dp 就好了.
记 \(f_{u,0/1}\) 表示 \(u\) 的子树内独立点集的最大权值和.
状态很好转移:
\[f_{u,0}=\sum_{v\in son_u} f_{v,1} \\ f_{u,1}=a_u+\sum_{v\in son_u} \max\{f_{v,0},f_{v,1}\} \]正解
那么我们来考虑怎么去优化这个东西吧!
我们看每次修改操作之后节点的 dp 值会怎么变化:
我们看式子就能发现,因为每个点的 dp 值是从他的子节点转移来的,所以当你调教了一个点,让它点权发生变化之后,只会对它到根节点路径上的这些点产生影响.
只有一条链上的 dp 值发生变化,但暴力算法每次把全部值都重新求了一遍,显然很暴力.(很优美bushi)
但是只改链上的 dp 值,这样做最坏情况也是 \(O(nm)\) 的.
诶那我们又注意到:有修改操作显然要用到数据结构呢(当然也不一定
想到之前各种树剖+线段树题(比如最经典的:重建)
我们想想怎么把这两个东西套起来(
最重要的一点就是树剖+线段树维护的信息必须要是满足幺半群性质的,也就是他要支持可合并性和结合律之类的.
而我们的 dp 式子是递推,不太好在线段树上维护.
那我们联想到:
形如 \(f_i=Af_{i-1}+Bf_{i-2}\) (可能有更多项)的线性递推式可以使用矩阵乘法优化,这样调教一下式子,就可以让每一次递推变成一次矩阵乘法,而矩阵乘法性质是很好的(没有时间出矩阵乘法的各种应用啦 摆烂!)
\[\left[\begin{array}{m} A & B\\ 1 & 0 \end{array}\right] \left[\begin{array}{m} f_{i-1}\\ f_{i-2} \end{array}\right] = \left[ \begin{array}{m} f_{i}\\ f_{i-1} \end{array} \right] \]再这样调教一下式子:
\[\left[\begin{array}{m} A & B\\ 1 & 0 \end{array}\right]^n \left[\begin{array}{m} f_{1}\\ f_{0} \end{array}\right] = \left[ \begin{array}{m} f_{n+1}\\ f_{n} \end{array} \right] \]诶然后矩阵快速幂不就完了
那我们可以轻松类比这题,让每一个递推变成一个矩阵!
(显然一点都不轻松不然这题就不会作为模板题形成一个单独的算法让我学了快一天 谁没见过还想得到这个啊啊啊.jpg)
那我们就通过查看题解的方式了解到:我们这个线性递推一个很重要一点就是每一次的 dp 式子是一样的,矩阵互不干扰,你如果修改了上一个节点的信息,下一个矩阵值也不一样了,显然是比较不好的一个东西
所以我们重新设计一下状态(设计状态是 dp 一个很重要的东西,可惜我不会)
我们令 \(g_{u,0/1}\) 表示 \(u\) 子树,除去了 \(hson[u]\) 的子树,这些点的独立点集的最大权值和,这样做呢,是不是更新的时候就可以让重链上节点的更新互不影响.
(\(hson[u]\) 是重儿子)
于是我们就有一个新的 dp 式子
\[f_{u,0}=\max\{f_{v,0},f_{v,1}\}+g_{u,0}\\ f_{u,1}=f_{v,0}+g_{u,1} \](其中 \(v=hson[u]\) )
呃呃是不是假的啊,这个 \(\max\) 根本去不掉,怎么写成矩阵乘法的形式啊?!不太可能吧ww
这又是一个套路(记下来!虽然快退役了但是感觉记下来总是好的)
类比加法与乘法,\(\max/\min\) 和加法也可以定义一种类似矩阵乘法的东西
怎么定义呢:
这是原版矩阵乘法:
\[C=AB\Rightarrow C_{i,j}=\sum_k A_{i,k}B_{k,j} \]那我们就把加法换成 \(\max\),把乘法换成加法
得到这个高仿 \(A\otimes B\) 运算( 电灯泡 / 雾 ):
\[C=A\otimes B\Rightarrow C_{i,j}=\max_k\{ A_{i,k}+B_{k,j}\} \]那他为什么满足矩阵乘法的各种性质呢(群论知识,我不会,长大后再学习!)
简单证明:
\[\begin{align*} & A \otimes (B \otimes C)=(A \otimes B) \otimes C\\ \Leftrightarrow\ & \forall i,j\ \max_k\{ A_{i,k}+\max_t\{ B_{k,t}+C_{t,j}\}\}= \max_k\{ \max_t\{ A_{i,t}+B_{t,k}\}+C_{k,j}\}\\ \Leftrightarrow\ & \forall i,j\ \max_k\{ \max_t\{A_{i,k}+ B_{k,t}+C_{t,j}\}\}= \max_k\{ \max_t\{ A_{i,t}+B_{t,k}+C_{k,j}\}\} \end{align*} \]交换 \(k,t\) 完显然成立
(注意到这个证明依赖于 \(\max\) 对加法的分配率,由此可以看出它也构成一个环)
嗯嗯然后呢我们把上面那个式子改成矩阵 \(\otimes\) 法的形式
上面的式子:
\[f_{u,0}=\max\{f_{v,0},f_{v,1}\}+g_{u,0}\\ f_{u,1}=f_{v,0}+g_{u,1} \]调教:
\[f_{u,0}=\max\{f_{v,0}+g_{u,0},f_{v,1}+g_{u,0}\}\\ f_{u,1}=\max\{f_{v,0}+g_{u,1},-\infin\} \]继续调教:
\[\left[\begin{array}{m} g_{u,0} & g_{u,0}\\ g_{u,1} & -\infin \end{array}\right] \otimes \left[\begin{array}{m} f_{v,0}\\ f_{v,1} \end{array}\right] = \left[ \begin{array}{m} f_{u,0}\\ f_{u,1} \end{array} \right] \]嗯嗯这就成一个可以在线段树上维护的东西了.
怎么实现呢?
首先先一次 dp 求出每个点的 \(f\) 和 \(g\),把 \(g\) 存到矩阵里
树剖,建线段树,注意线段树上是按 dfs 序存
这时候每次怎么求出答案呢
答案就是 \(\max\{f_{1,0},f_{1,1}\}\)
其实就是以 \(1\) 号节点为顶的这条重链上所有矩阵的 \(\otimes\) 积(
那我们每次修改操作,单点修改完 \(u\),哪些点的矩阵会发生变化?
显然只有每次的 \(fa[top[u]]\),因为根据定义, \(g\) 只会从轻儿子转移来.
因为从 \(u\) 到 \(1\) 最多经过 \(\log n\) 条重链,所以总体的时间复杂度是 \(O(m\log^2n)\)
实现的细节看代码:
代码
//https://www.luogu.com.cn/problem/P4719
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N = 1e5+5, inf = 1e9;
struct Mat //Matrix
{
int r,c;
int d[2][2];
void print() { for(int i=0; i<r; i++) { for(int j=0; j<c; j++) if(d[i][j]>-inf) printf("%4d ",d[i][j]); else printf("-inf "); puts(""); } }
void st(vector<vector<int> > v) //初始化矩阵
{
r=v.size(); c=v[0].size();
for(int i=0; i<r; i++) for(int j=0; j<c; j++) d[i][j]=v[i][j];
}
int * operator[] (int x) { return d[x]; } //小科技,可以把A.d[i][j]变成A[i][j]
friend Mat operator ^ (Mat a, Mat b) //重定义矩阵乘法
{
assert(a.c==b.r); Mat c; c.r=a.r; c.c=b.c;
for(int i=0; i<a.r; i++) for(int j=0; j<b.c; j++)
{
c[i][j]=-inf;
for(int k=0; k<a.c; k++) c[i][j]=max(c[i][j],a[i][k]+b[k][j]);
}
return c;
}
};
int n,m;
int a[N];
vector<int> g[N];
int fa[N],sz[N],d[N],dfn[N],rk[N],hson[N],top[N],btm[N],ct; //树剖
Mat F[N],G[N];
struct Sagiri //Segment tree
{
int l,r;
Mat dat;
} t[N<<2];
#define ls (p<<1)
#define rs (p<<1)|1
void build(int p, int l, int r)
{
t[p].l=l,t[p].r=r;
if(l==r) { t[p].dat=G[rk[l]]; return ; }
int mid=(l+r)>>1;
build(ls,l,mid),build(rs,mid+1,r);
t[p].dat=t[ls].dat^t[rs].dat;
}
void change(int p, int x, const Mat &v) //结构体记得传实参
{
if(t[p].l==t[p].r) { t[p].dat=v; return ; }
int mid=(t[p].l+t[p].r)>>1;
if(x<=mid) change(ls,x,v);
else change(rs,x,v);
t[p].dat=t[ls].dat^t[rs].dat;
}
Mat query(int p, int l, int r)
{
if(t[p].l==l and t[p].r==r) return t[p].dat;
int mid=(t[p].l+t[p].r)>>1;
if(r<=mid) return query(ls,l,r);
else if(l>mid) return query(rs,l,r);
else return query(ls,l,mid)^query(rs,mid+1,r); /**/
}
void slpf1(int u, int prv)
{
fa[u]=prv; sz[u]=1;
for(int v:g[u])
{
if(v==prv) continue;
d[v]=d[u]+1; slpf1(v,u); sz[u]+=sz[v];
if(sz[v]>sz[hson[u]]) hson[u]=v;
}
}
void slpf2(int u)
{
dfn[u]=++ct; rk[dfn[u]]=u; btm[u]=u;
F[u].st({{0},{a[u]}});
G[u].st({{0,0},{a[u],-inf}});
if(hson[u])
{
int v=hson[u]; top[v]=top[u]; slpf2(v); btm[u]=btm[v];
F[u][0][0]+=max(F[v][0][0],F[v][1][0]);
F[u][1][0]+=F[v][0][0];
}
for(int v:g[u])
{
if(v==fa[u] or v==hson[u]) continue;
top[v]=v; slpf2(v);
F[u][0][0]+=max(F[v][0][0],F[v][1][0]);
F[u][1][0]+=F[v][0][0];
G[u][0][0]+=max(F[v][0][0],F[v][1][0]);
G[u][1][0]+=F[v][0][0]; //G只由轻儿子转移
}
G[u][0][1]=G[u][0][0];
}
void init()
{
scanf("%d%d",&n,&m);
for(int i=1; i<=n; i++) scanf("%d",&a[i]);
for(int i=1,u,v; i<n; i++) scanf("%d%d",&u,&v),g[u].push_back(v),g[v].push_back(u);
slpf1(1,0); top[1]=1; slpf2(1);
build(1,1,n);
}
void update(int u, int w)
{
G[u][1][0]+=w-a[u]; a[u]=w; //先修改u的矩阵
while(u)
{
Mat st=query(1,dfn[top[u]],dfn[btm[u]]);
change(1,dfn[u],G[u]); //求出修改前和修改后轻儿子的贡献,就可以改变fa[top[u]]的矩阵
Mat ed=query(1,dfn[top[u]],dfn[btm[u]]);
u=fa[top[u]];
G[u][0][0]+=max(ed[0][0],ed[1][0])-max(st[0][0],st[1][0]);
G[u][0][1]=G[u][0][0];
G[u][1][0]+=ed[0][0]-st[0][0];
}
}
void answer_query()
{
for(int i=1,u,w; i<=m; i++)
{
scanf("%d%d",&u,&w);
update(u,w);
Mat ans=query(1,dfn[1],dfn[btm[1]]); //答案即这条重链的矩阵之积
printf("%d\n",max(ans[0][0],ans[1][0]));
}
}
int main()
{
init();
answer_query();
return 0;
}
小练习:
这题没差多少
只是稍微变化了一下 dp 式子
#include <bits/stdc++.h>
#define ll long long
//#define debug
using namespace std;
const int N = 1e5+5;
const ll inf = 1e17;
struct Mat
{
static const int B = 2;
ll d[2][2];
ll * operator[](int x) { return d[x]; }
void init(vector<vector<ll> > v) { for(int i=0; i<B; i++) for(int j=0; j<B; j++) d[i][j]=v[i][j]; }
friend Mat operator^ (Mat a, Mat b)
{
Mat c;
for(int i=0; i<B; i++)
for(int j=0; j<B; j++)
{
c[i][j]=inf;
for(int k=0; k<B; k++)
c[i][j]=min(c[i][j],a[i][k]+b[k][j]);
}
return c;
}
};
int n,m;
ll a[N];
vector<int> g[N];
int fa[N],d[N],dfn[N],rk[N],sz[N],hson[N],top[N],btm[N],ct;
ll F[N][2];
Mat G[N];
struct Sagiri
{
int l,r;
Mat dat;
} t[N<<2];
#define ls (p<<1)
#define rs (p<<1)|1
void build(int p, int l, int r)
{
t[p].l=l,t[p].r=r;
if(l==r) { t[p].dat=G[rk[l]]; return ; }
int mid=(l+r)>>1;
build(ls,l,mid),build(rs,mid+1,r);
t[p].dat=t[ls].dat^t[rs].dat;
}
void update(int p, int x, const Mat &v)
{
if(t[p].l==t[p].r) { t[p].dat=v; return ; }
int mid=(t[p].l+t[p].r)>>1;
if(x<=mid) update(ls,x,v);
else update(rs,x,v);
t[p].dat=t[ls].dat^t[rs].dat;
}
Mat query(int p, int l, int r)
{
if(t[p].l==l and t[p].r==r) return t[p].dat;
int mid=(t[p].l+t[p].r)>>1;
if(r<=mid) return query(ls,l,r);
else if(l>mid) return query(rs,l,r);
else return query(ls,l,mid)^query(rs,mid+1,r);
}
void slpf1(int u, int prv)
{
fa[u]=prv; sz[u]=1;
for(int v:g[u])
{
if(v==prv) continue;
d[v]=d[u]+1; slpf1(v,u); sz[u]+=sz[v];
if(sz[v]>sz[hson[u]]) hson[u]=v;
}
}
void slpf2(int u)
{
dfn[u]=++ct,rk[dfn[u]]=u; btm[u]=u;
F[u][0]=0,F[u][1]=a[u];
G[u].init({{inf,0},{a[u],a[u]}});
if(hson[u])
{
int v=hson[u];
top[v]=top[u]; slpf2(v); btm[u]=btm[v];
F[u][0]+=F[v][1]; F[u][1]+=min(F[v][0],F[v][1]);
}
for(int v:g[u])
{
if(v==fa[u] or v==hson[u]) continue;
top[v]=v; slpf2(v);
F[u][0]+=F[v][1]; F[u][1]+=min(F[v][0],F[v][1]);
G[u][0][1]+=F[v][1]; G[u][1][0]+=min(F[v][0],F[v][1]);
}
G[u][1][1]=G[u][1][0];
}
void change(int u, ll w)
{
G[u][1][1]+=w-a[u]; G[u][1][0]=G[u][1][1]; a[u]=w;
while(u!=0)
{
Mat st=query(1,dfn[top[u]],dfn[btm[u]]);
update(1,dfn[u],G[u]);
Mat ed=query(1,dfn[top[u]],dfn[btm[u]]);
u=fa[top[u]];
G[u][1][0]+=min(ed[0][1],ed[1][1])-min(st[0][1],st[1][1]);
G[u][0][1]+=ed[1][1]-st[1][1];
G[u][1][1]=G[u][1][0];
}
}
int main()
{
string s; scanf("%d%d",&n,&m); cin>>s;
for(int i=1; i<=n; i++) scanf("%lld",&a[i]);
for(int i=1,u,v; i<n; i++) scanf("%d%d",&u,&v),g[u].push_back(v),g[v].push_back(u);
slpf1(1,0),top[1]=1,slpf2(1); build(1,1,n);
for(int i=1,qa,x,qb,y; i<=m; i++)
{
scanf("%d%d%d%d",&qa,&x,&qb,&y);
ll ta=a[qa],tb=a[qb];
change(qa,x?0:inf); change(qb,y?0:inf);
Mat res=query(1,dfn[1],dfn[btm[1]]);
ll ans=min(res[0][1],res[1][1]);
if(x) ans+=ta; if(y) ans+=tb;
printf("%lld\n",ans>(inf/2)?-1:ans);
change(qa,ta); change(qb,tb);
}
return 0;
}
标签:hson,int,max,矩阵,笔记,array,动态,dp
From: https://www.cnblogs.com/copper-carbonate/p/16886598.html