2024.8.15

2024.8.15 【这雨生于天，死于地，中间的过程就是人生。】

Thursday 七月十二

---

动态DP （DDP）

我们选择从一道模板题讲起

[P4719 【模板】"动态 DP"&动态树分治](P4719 [模板]"动态 DP"&动态树分治 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn))

【模板】"动态 DP"&动态树分治

题目描述

给定一棵 \(n\) 个点的树，点带点权。

有 \(m\) 次操作，每次操作给定 \(x,y\)，表示修改点 \(x\) 的权值为 \(y\)。

你需要在每次操作之后求出这棵树的最大权独立集的权值大小。

输入格式

第一行有两个整数，分别表示结点个数 \(n\) 和操作个数 \(m\)。

第二行有 \(n\) 个整数，第 \(i\) 个整数表示节点 \(i\) 的权值 \(a_i\)。

接下来 \((n - 1)\) 行，每行两个整数 \(u, v\)，表示存在一条连接 \(u\) 与 \(v\) 的边。

接下来 \(m\) 行，每行两个整数 \(x,y\)，表示一次操作，修改点 \(x\) 的权值为 \(y\)。

输出格式

对于每次操作，输出一行一个整数表示答案。

样例 #1

样例输入 #1

10 10
-11 80 -99 -76 56 38 92 -51 -34 47 
2 1
3 1
4 3
5 2
6 2
7 1
8 2
9 4
10 7
9 -44
2 -17
2 98
7 -58
8 48
3 99
8 -61
9 76
9 14
10 93

样例输出 #1

186
186
190
145
189
288
244
320
258
304

提示

数据规模与约定

• 对于 \(30\%\) 的数据，保证 \(n,m\le 10\)。
• 对于 \(60\%\) 的数据，保证 \(n,m\le 10^3\)。
• 对于 \(100\%\) 的数据，保证 \(1\le n,m\le 10^5\)，\(1 \leq u, v , x \leq n\)，\(-10^2 \leq a_i, y \leq 10^2\)。

前置知识：

1. 树链剖分
2. 树上DP
3. 广义矩阵乘法
4. 最大权独立集

首先考虑当没有修改时怎么作，

我们发现基本就是一个树上dp的板子题

可以很快列出DP方程

设f[i][0]表示i不在独立集时以i为根节点的最大权值

f[i][1]则表示i在独立集中的最大权值

\[f[i][0] = \sum_{v\in i.son}max(f[v][0],f[v][1])\\ f[i][1] = \sum_{v\in i.son}f[v][0]+val(i) \]

但是现在他又加入了修改，就很难搞

我们考虑暴力修改

不难发现对于一个节点的修改，只有这个结点所在的链会受到影响

只对它的所有祖先修改一下就可以了

但是复杂度显然不是很好看

考虑使用树链剖分实现快速向上跳

itn fa[oo],dep[oo],siz[oo],son[oo];
__inline void dfs1(itn x,int f){
    fa[x] = f;siz[x] = 1;dep[x] = dep[f]+1;
    for (itn i=head[x];~i;i=st[i].nxt){
        itn v = st[i].v;
        if (v==f) continue;
        dfs1(v,x);
        if (siz[v]>siz[son[x]])son[x] = v;
    }
}
int tim,dfn[oo],nfd[oo],top[oo],end[oo];
__inline void dfs2(int x,itn t){
    top[x] = t;
    dfn[x] = ++tim;nfd[tim] = x;
    if (!son[x]){
        f[1][x] = val[x];
        g[x] = ident;
        end[x] = x;
        return void();
    }
    g[x](1,0) = val[x],g[x](1,1) = -inf;
    dfs2(son[x],t);
    end[x] = end[son[x]];
    for (itn i=head[x];~i;i=st[i].nxt){
        itn v = st[i].v;
        if (v==fa[x]||v==son[x])  continue;
        dfs2(v,v);
        g[x](0,0) = g[u](0,1)+=max(f[0][son[x]],f[1][son[x]]);
        g[u](1,0)+=f[0][v];
    }
    f[0][u]=g[u](0,0)+max(f[0][son[u]],f[1][son[u]]);
    f[1][u]=g[u](1,0)+f[0][son[u]];
}

但是我们很快会发现一点点小问题。。。

树剖的话。。。需要使用线段树维护，

那么线段树。。。。

但是能使用线段树维护的数据有要求必须满足可叠加性，或者说，结合律

这就很难搞了对吧

显然DP过程需要连续性

我们并不能跳过前面的DP值直接推出后面的值

真的不能吗？？？？

那么我们就该考虑一下，什么东西，能够让DP具有加法，乘法一样类似的性质呢？

我们不知不觉就一定会用聪明的小脑袋想到斐波那契数列数列的矩阵快速幂加速计算

那么矩阵乘法是否满足结合率呢？

\[\begin{pmatrix}\begin{vmatrix} a\ b \\ c \ d \end{vmatrix} \times \begin{vmatrix} e \ f \\ g \ h \end{vmatrix} \end{pmatrix} \times \begin{vmatrix} i\ j \\ k\ l \end{vmatrix} \\=\begin{vmatrix} (ae+bg)i+(af+bh)k \ \ (ae+bg)j+(af+bh)l \\ (ce+dg)i+(cf+dh)k \ \ (ce+dg)j+(cf+dh)l \end{vmatrix} \]

\[\begin{vmatrix} a\ b \\ c \ d \end{vmatrix} \times\begin{pmatrix} \begin{vmatrix} e \ f \\ g \ h \end{vmatrix} \times \begin{vmatrix} i\ j \\ k\ l \end{vmatrix}\end{pmatrix} \\=\begin{vmatrix} a(ei+fk)+b(gi+hk) \ \ a(ej+fl)+b(gj+hl) \\ c(ei+fk)+d(gi+hk) \ \ c(ei+fl)+d(gj+hl) \end{vmatrix} \]

显然，矩阵能满足结合律，对其涉及的两种操作：

加法和乘法有一定=要求

1. 加法之间满足交换律
2. 乘法之间满足交换律
3. 加法对乘法有分配律

现在我们考虑一下DP的转移过程长什么样子：

\[f[i][0] = \sum_{v\in i.son}max(f[v][0],f[v][1])\\ f[i][1] = \sum_{v\in i.son}f[v][0]+val(i) \]

显然这样转移，碰到菊花图复杂度会直接爆掉

不过既然我们已经使用了树链剖分，倒不如直接将儿子分开来算，

分成轻儿子和重儿子

我们设g[i][0]为当i为根的时候，其属于轻子树的\(\sum_{v\in i.son}max(f[v][0],f[v][1])\)

g[i][1]为i为根时，属于轻子树的\(\sum_{v\in i.son}f[v][0]\)

我们就可以将DP的转移改变一下

\[\begin{bmatrix} g[i][0]+max(f[i+1])\\(g[i][1]+val[i])+f[i+1][0] \end{bmatrix} = \begin{bmatrix} f[i][0]\\f[i][1] \end{bmatrix} \]

转移过程肉眼可见的鬼畜起来了

下面考虑构造一个矩阵使其满足DP的递推

在这个矩阵的转移中涉及到的两种操作，分别是求和以及取max，

我们发现这两种操作都是满足交换律的，并且加法对于取max运算是有分配律的

所以我们便可顺理成章的使用矩阵了

现在，我们的矩阵运算变成了这样：

\[\begin{bmatrix} a \ b \\ c \ d \end{bmatrix} \otimes \begin{bmatrix} e \ f\\ g \ h \end{bmatrix} = \begin{bmatrix} max(a+e,b+g) \ max(a+f,b+h)\\ max(c+e,d+g) \ max(c+f,d+h) \end{bmatrix} \]

所以，我们可以做出如下构建:

\[\begin{bmatrix} f[i][0]\\f[i][1] \end{bmatrix}\\ = \begin{bmatrix} g[i][0]+max(f[i+1][0],f[i+1][1])\\(g[i][1]+val[i])+f[i+1][0] \end{bmatrix} \\= \begin{bmatrix} max(g[i][0]+f[i+1][0],g[i][0]+f[i+1][1]) \\max(g[i][1]+val[i]+f[i+1][0],-\infty+f[i+1][1]) \end{bmatrix} \\= \begin{bmatrix} g[i][0] \ \ \ \ \ \ \ \ \ \ \ \ \ g[i][0]\\ g[i][1]+val[i] \ -\infty \end{bmatrix} \otimes \begin{bmatrix} f[i+1][0] \\ f[i+1][1] \end{bmatrix} \]

那么现在非常好，我们已经搞定一个类似矩阵快速幂的操作了

template <int row,int col>
struct mat{
    int r,c;
    int ele[row][col];
    mat():r(row),c(col) {}
    int& operator()(int a,int b) { return ele[a][b]; }
};
template <int m,int n,int p>
auto operator*(mat<m,n> m1,mat<n,p> m2){
    mat<m,p> ret;
    memset(ret.ele,0xcf,sizeof(ret.ele));
    for(int i=0;i<m;i++)        
        for(int k=0;k<n;k++)
            for(int j=0;j<p;j++)
                ret(i,j)=max(ret(i,j),m1(i,k)+m2(k,j));
    return ret;
}
mat<2,2> ident,g[oo];template <int row,int col>
struct mat{
    int r,c;
    int ele[row][col];
    mat():r(row),c(col) {}
    int& operator()(int a,int b) { return ele[a][b]; }
};
template <int m,int n,int p>
auto operator*(mat<m,n> m1,mat<n,p> m2){
    mat<m,p> ret;
    memset(ret.ele,0xcf,sizeof(ret.ele));
    for(int i=0;i<m;i++)        
        for(int k=0;k<n;k++)
            for(int j=0;j<p;j++)
                ret(i,j)=max(ret(i,j),m1(i,k)+m2(k,j));
    return ret;
}
mat<2,2> ident,g[oo];

下面考虑如何查询

我们记录每个节点的值显然不现实，所以查询的时候在上面现算就好

叶子节点的\(f\)值是已知的，\(f[n][0] = 0,f[n][1] = val[n]\),使用重链上维护即可

我们对每个叶子节点，构建关于他的初始矩阵，即：

\[\begin{bmatrix} 0\\val[end(x)] \end{bmatrix} \]

同时，方便起见，我们将线段树所有叶子节点放到单位矩阵中，单位矩阵长这个样子：

\[\begin{bmatrix} 0 \ \ -\infty\\ -\infty \ \ \ 0 \end{bmatrix} \]

所以，我们每次查询，就可以查询从当前节点到重链底端，用$query(dfn[x],dfn[end[x]])\times\begin{bmatrix}0\val[0]\end{bmatrix} $

求出\(\begin{bmatrix}f[x][0]\\f[x][1]\end{bmatrix}\)

那么修改操作呢？

只修改x结点不难，但是要考虑到x所在重链顶端不一定为树根，所以还要考虑祖先的变化

所以父亲所在节点的\(g\)数组会改变，需要同时修改节点的矩阵，这是一个类似递归的过程

\(g\)是由轻儿子的\(f\)求出的，所以要先查出来这个\(f\)

修改的总流程长这样：

1. 将矩阵$\begin{bmatrix}
   g[i][0] \ \ \ \ \ \ \ \ \ \ \ \ \ g[i][0]\
   g[i][1]+val[i] \ -\infty
   \end{bmatrix} $ 左下角减去老的\(val[x]\)加上新的\(val[x]\),算出一个新矩阵tmp
2. 修改前查询出一个$old = \begin{bmatrix}f[top[x]][0]\f[top[x]][1]\end{bmatrix} $
3. 修改节点为tmp，在查询出一个$new =\begin{bmatrix}f[x][0]\f[x][1]\end{bmatrix} $
4. 将x置为\(fa[top[x]]\),如果x为根节点的父亲（0），结束
5. x的矩阵$\begin{bmatrix}
   g[x][0] \ \ \ \ \ \ \ \ \ \ \ \ \ g[x][0]\
   g[x][1]+val[x] \ -\infty
   \end{bmatrix} \(,第一行减去max(old),加max(new),左下减去\)old.f[top[x]][0]\(加\)new.f[top[x]][0]$,记录到tmp
6. 返回步骤2

代码长这样：

struct Node{
    int l,r;
    mat<2,2> val;
}sgt[oo<<2];
#define ls(x) (x<<1)
#define rs(x) (x<<1|1)
#define pushup(x) sgt[x].val=sgt[ls(x)].val*sgt[rs(x)].val
void build(int l,int r,int k=1){
    sgt[k].l=l,sgt[k].r=r;
    if(l==r){
        sgt[k].val=g[nfd[l]];
        return;
    }
    int m = (l+r)>>1;
    build(l,m,ls(k));
    build(m+1,r,rs(k));
    pushup(k);
}
auto query(int x,int y,int k=1){
    int l=sgt[k].l,r=sgt[k].r;
    if(x<=l&&y>=r) return sgt[k].val;
    int m = (l+r)>>1;
    mat<2,2> ret = ident;
    if(x<=m) ret = query(x,y,ls(k));
    if(y>m) ret = ret * query(x,y,rs(k));
    return ret;
}
void modify(int x,int p,int k=1){
    int l=sgt[k].l,r=sgt[k].r;
    if(l==r){
        sgt[k].val=g[p];
        return;
    }
    int m = (l+r)>>1;
    if(x<=m) modify(x,p,ls(k));
    else modify(x,p,rs(k));
    pushup(k);
}
auto queryt(int x){
    mat<2,1> tmp;
    tmp(0,0)=0,tmp(1,0)=val[edf[x]];
    return query(dfn[x],dfn[edf[x]]) * tmp;
}
void modifyt(int x,int z){
    mat<2,1> od,nw;
    g[x](1,0)+=z-val[x];
    od = queryt(top[x]);
    val[x]=z;
    while(x){
        modify(dfn[x],x);
        nw = queryt(top[x]);
        x = fa[top[x]];
        g[x](0,0)=g[x](0,1)+=max(nw(0,0),nw(1,0))-max(od(0,0),od(1,0));
        g[x](1,0)+=nw(0,0)-od(0,0);
        od = queryt(top[x]);
    }
}

下面直接上所有代码咯）

//2024.8.15
//by white_ice
//P4719 【模板】"动态 DP"&动态树分治
#include<bits/stdc++.h>
using namespace std;
constexpr int oo = 1e5+5;
constexpr int inf = 0x3f3f3f3f;
struct nod {int to,next;}st[oo<<1];
int tot,head[oo];
inline void add(int u,int v){
    st[tot]=(nod){v, head[u]};
    head[u]=tot++;
}

template <int row,int col>
struct mat{
    int r,c;
    int ele[row][col];
    mat():r(row),c(col) {}
    int& operator()(int a,int b) { return ele[a][b]; }
};
template <int m,int n,int p>
auto operator*(mat<m,n> m1,mat<n,p> m2){
    mat<m,p> ret;
    memset(ret.ele,0xcf,sizeof(ret.ele));
    for(int i=0;i<m;i++)        
        for(int k=0;k<n;k++)
            for(int j=0;j<p;j++)
                ret(i,j)=max(ret(i,j),m1(i,k)+m2(k,j));
    return ret;
}
mat<2,2> ident,g[oo];
int val[oo],f[2][oo];

int fa[oo],dep[oo],siz[oo],son[oo];
void dfs1(int u,int f){
    fa[u]=f,siz[u]=1,dep[u]=dep[f]+1;
    for(int i=head[u];~i;i=st[i].next){
        int v = st[i].to;
        if(v==f) continue;
        dfs1(v,u);
        siz[u]+=siz[v];
        if(siz[v]>siz[son[u]])
            son[u]=v;
    }
}

int tim,dfn[oo],nfd[oo],top[oo],edf[oo];
void dfs2(int u,int t){
    top[u]=t,dfn[u]=++tim,nfd[tim]=u;
    if(!son[u]){
        f[1][u]=val[u];
        g[u]=ident;
        edf[u]=u;
        return;
    }
    g[u](1,0)=val[u],g[u](1,1)=-inf;
    dfs2(son[u],t);
    edf[u]=edf[son[u]];
    for(int i=head[u];~i;i=st[i].next){
        int v = st[i].to;
        if(v==fa[u]||v==son[u]) continue;
        dfs2(v,v);
        g[u](0,0)=g[u](0,1)+=max(f[0][v],f[1][v]);
        g[u](1,0)+=f[0][v];
    }
    f[0][u]=g[u](0,0)+max(f[0][son[u]],f[1][son[u]]);
    f[1][u]=g[u](1,0)+f[0][son[u]];
}
struct Node{
    int l,r;
    mat<2,2> val;
}sgt[oo<<2];
#define ls(x) (x<<1)
#define rs(x) (x<<1|1)
#define pushup(x) sgt[x].val=sgt[ls(x)].val*sgt[rs(x)].val
void build(int l,int r,int k=1){
    sgt[k].l=l,sgt[k].r=r;
    if(l==r){
        sgt[k].val=g[nfd[l]];
        return;
    }
    int m = (l+r)>>1;
    build(l,m,ls(k));
    build(m+1,r,rs(k));
    pushup(k);
}
auto query(int x,int y,int k=1){
    int l=sgt[k].l,r=sgt[k].r;
    if(x<=l&&y>=r) return sgt[k].val;
    int m = (l+r)>>1;
    mat<2,2> ret = ident;
    if(x<=m) ret = query(x,y,ls(k));
    if(y>m) ret = ret * query(x,y,rs(k));
    return ret;
}
void modify(int x,int p,int k=1){
    int l=sgt[k].l,r=sgt[k].r;
    if(l==r){
        sgt[k].val=g[p];
        return;
    }
    int m = (l+r)>>1;
    if(x<=m) modify(x,p,ls(k));
    else modify(x,p,rs(k));
    pushup(k);
}
auto queryt(int x){
    mat<2,1> tmp;
    tmp(0,0)=0,tmp(1,0)=val[edf[x]];
    return query(dfn[x],dfn[edf[x]]) * tmp;
}
void modifyt(int x,int z){
    mat<2,1> od,nw;
    g[x](1,0)+=z-val[x];
    od = queryt(top[x]);
    val[x]=z;
    while(x){
        modify(dfn[x],x);
        nw = queryt(top[x]);
        x = fa[top[x]];
        g[x](0,0)=g[x](0,1)+=max(nw(0,0),nw(1,0))-max(od(0,0),od(1,0));
        g[x](1,0)+=nw(0,0)-od(0,0);
        od = queryt(top[x]);
    }
}
int main(int argc, char const *argv[]){
    clock_t c1 = clock();
    memset(head,-1,sizeof(head));
    ident(1,0)=ident(0,1)=-inf;
    int n,m;
    cin>>n>>m;
    for(int i=1;i<=n;i++) cin>>val[i];
    for(int i=1;i<n;i++){
        int u,v;
        cin>>u>>v;
        add(u,v),add(v,u);
    }
    dfs1(1,0);
    dfs2(1,1);
    build(1,n);
    while(m--){
        int x,z;
        cin>>x>>z;
        modifyt(x,z);
        auto ans = queryt(1);
        cout<< max(ans(0,0),ans(1,0)) << endl;
    }
    cerr << "Time:" << clock() - c1 << "ms" << endl;
    return 0;
}