DP 套 DP

引入

本文将介绍 DP 套 DP 的思想,并通过两道例题展示它如何应用到具体问题中.

思想

所谓「DP 套 DP」,实际上指的是在动态规划的过程中,把一个子问题的求解过程(通常是一个 DP)抽象成一个自动机(DFA),然后在这个自动机的基础上再设计一层新的 DP 的方法.

这种技巧主要应用于一类 序列计数、概率 或 期望 问题.一个典型的问题具有如下结构:

  • 给定字符集 Σ\Sigma 和它上面一个长度为 nn 的「合法序列」的集合 A⊆ΣnA\subseteq\Sigma^n.根据字符集不同,序列可以是二进制串、数字串、状态序列等.
  • 对于每一个具体的序列 s∈Σns\in\Sigma^n,可以通过动态规划来判断它是否合法(即 s∈As\in A),计算它的权值或求出相关值.
  • 最终,我们希望统计集合 AA 中所有序列的数目、总权值、期望值等.

这时,枚举所有序列是不可行的.于是我们考虑将「判断一个序列是否合法」的过程(即内层 DP)抽象为一个 确定有限状态自动机(DFA).一般地,对于固定的序列 s∈Σns\in\Sigma^n,内层 DP 的状态函数可以表示为 g(i,x;s)g(i,x;s),即已经处理完序列 ss 的长度为 ii 的前缀,且其他状态分量为 xx 时某个量的取值.相应地,内层 DP 的状态转移方程为

g(i,⋅;s)=G(g(i−1,⋅;s),si).g(i,\cdot;s) = G(g(i-1,\cdot;s),s_i).

也就是说,函数 g(i,⋅;s)g(i,\cdot;s) 由之前的函数 g(i−1,⋅;s)g(i-1,\cdot;s) 和当前的字符 sis_i 唯一确定.如果将函数 g(i,⋅;s)g(i,\cdot;s) 看作是自动机的一个状态,那么,内层 DP 的状态转移方程就给出了自动机的一个转移.因此,内层 DP 对应的自动机 (Q,Σ,δ,q0,F)(Q,\Sigma,\delta,q_0,F) 结构如下:

  • 状态集合 QQ 就是所有可能的 s∈Σns\in\Sigma^n 和 i=0,1,⋯ ,ni=0,1,\cdots,n 对应的函数 g(i,⋅;s)g(i,\cdot;s) 的集合;
  • 转移函数 δ:Q×Σ→Q\delta:Q\times\Sigma\to Q 就是内层 DP 的状态转移方程中的 GG;
  • 起始状态 q0q_0 通常是显然的,即内层 DP 的初始状态;
  • 接受状态集合 FF 就对应着所有的合法序列 s∈As\in A.

函数 g(i,⋅;s)g(i,\cdot;s) 本身可能相当复杂,因此,在处理具体问题时,通常需要进行 状态压缩 或结合 DFA 最小化的技巧来压缩状态空间.这也是 DP 套 DP 相较于暴力 DP 能够显著降低时空复杂度的主要原因.

将内层 DP 抽象为 DFA 后,就可以在这个 DFA 上设计一个新的 DP 用于求解原问题,即外层 DP.为方便表述,以简单的计数问题为例.外层 DP 的状态函数定义为 f(i,q)f(i,q),即处理到长度为 ii 的前缀时,到达 DFA 中状态 q∈Qq\in Q 的前缀的数目.它的状态转移方程为

f(i,q)=∑c∈Σ∑q′∈Q:δ(q′,c)=qf(i−1,q′).f(i,q) = \sum_{c\in\Sigma}\sum_{q'\in Q:\delta(q',c)=q} f(i-1,q').

起始状态当然是 f(0,q0)f(0,q_0),而最终要求的答案通常可以根据 {f(n,q):q∈F}\{f(n,q):q\in F\} 简单计算得到.外层 DP 实际上是 DAG 上 DP 的特殊情形.

例题

接下来的两个例题会详细说明 DP 套 DP 的一般做法.

例一

Hero meet devil

给定一个字符集为 ACGT 的字符串 SS,且 ∣S∣≤15|S|\le 15.对于每个 0≤i≤∣S∣0\leq i \leq |S|,求有多少个长度为 mm,字符集 ACGT 的字符串 TT,满足它与 SS 的最长公共子序列长度为 ii.

题解

我们首先会想到一个 DP:设 fi,jf_{i,j} 表示长度为 ii 的 TT 中,和 SS 的最长公共子序列的长度为 jj 的方案数.但是这样无法转移,我们发现主要的问题是,我们不知道这个最长公共子序列对应的是哪些字符.

考虑朴素求最长公共子序列的过程.设 gi,jg_{i,j} 表示 TT 的前 ii 位和 SS 的前 jj 位,它们的最长公共子序列的长度,就有

gi,j=max⁡{gi−1,j,gi,j−1,gi−1,j−1+[Ti=Sj]}.g_{i,j} = \max\{g_{i-1,j},g_{i,j-1},g_{i-1,j-1}+[T_i=S_j]\}.

我们发现,对于一个 ii,只需要记录 gig_i 这个一维数组每一位的值,就能准确维护当前 SS 与 TT 前 ii 位最长公共子序列的状态.因为 SS 长度只有 1515,我们发现这个思想是可行的.

于是重新设状态 fi,xf_{i,x} 表示对于长度为 ii 的 TT,与 SS 的 DP 数组(就是 gig_i)状态为 xx 的方案数.这个 DP 看起来状态数很多,然而我们发现 gi,j−gi,j−1∈{0,1}g_{i,j}-g_{i,j-1}\in\{0,1\},就可以维护 gig_i 的差分数组,状态数是 2∣S∣2^{|S|} 的.

现在思考怎么转移.容易发现,如果我们知道了 gig_i 这个数组,也知道了 Ti+1T_{i+1},就能通过朴素 LCS 转移(即前文的 DP 方程)求出 gi+1g_{i+1}.于是朴素的 LCS 就成为了帮助 ff 转移的内层 DP.

因此,我们枚举 Ti+1T_{i+1},计算出 xx 转移后的状态 x′x',再将 fi+1,x′f_{i+1,x'} 加上 fi,xf_{i,x} 就可以完成外层 DP 的状态转移.最后,我们记录 ansi\textit{ans}_i 为 LCS 长度为 ii 的答案,枚举每个状态 SS,anspopcount⁡(S)\textit{ans}_{\operatorname{popcount}(S)} 加上 fm,Sf_{m,S} 即可.

参考代码
#include <cstring>
#include <iostream>
#include <string>
using namespace std;

const int MOD = 1e9 + 7;
const int MAXN = 100;

int T, a[MAXN + 1], h[MAXN + 1];
int g[MAXN + 1], m, n;
int nxt[1 << 16]
       [4];  // 状态转移表:nx[mask][ch] 表示 mask 状态下追加 ch 后的新状态
int f[1001]
     [1 << 16];  // DP 数组:f[i][mask] 表示长度为 i,处于状态 mask 的方案数
int ans[MAXN +
        1];  // 最终答案:ans[i] 表示长度为 m 时,包含 i 个匹配位置的方案数

void add(int &x, int y) { x = (x + y) % MOD; }

int calc(int mask, int x) {  // 给定当前 mask 状态和当前字符
                             // x(0~3),返回新状态
  int res = 0;
  for (int i = 0; i < n; ++i) g[i + 1] = g[i] + ((mask >> i) & 1);

  for (int i = 1; i <= n; ++i) {
    if (a[i] == x)
      h[i] = g[i - 1] + 1;
    else
      h[i] = max(h[i - 1], g[i]);
  }

  for (int i = 1; i <= n; ++i)
    if (h[i] > h[i - 1]) res |= (1 << (i - 1));

  return res;
}

int popcount(int x) {  // 计算 mask 中有多少个 1(即匹配了多少位置)
  int res = 0;
  while (x) {
    res += x & 1;
    x >>= 1;
  }
  return res;
}

int main() {
  ios::sync_with_stdio(false);
  cin.tie(nullptr);

  cin >> T;
  while (T--) {
    string S;
    cin >> S >> m;
    n = S.length();

    memset(f, 0, sizeof(f));
    memset(ans, 0, sizeof(ans));
    memset(a, 0, sizeof(a));

    for (int i = 0; i < n; ++i) {
      if (S[i] == 'A') a[i + 1] = 0;
      if (S[i] == 'C') a[i + 1] = 1;
      if (S[i] == 'G') a[i + 1] = 2;
      if (S[i] == 'T') a[i + 1] = 3;
    }

    int lim = (1 << n);
    for (int i = 0; i < lim; ++i)
      for (int j = 0; j < 4; ++j) nxt[i][j] = calc(i, j);

    f[0][0] = 1;
    for (int i = 0; i < m; ++i)
      for (int j = 0; j < lim; ++j)
        if (f[i][j])
          for (int k = 0; k < 4; ++k) add(f[i + 1][nxt[j][k]], f[i][j]);

    for (int i = 0; i < lim; ++i) add(ans[popcount(i)], f[m][i]);

    for (int i = 0; i <= n; ++i) cout << ans[i] << '\n';
  }

  return 0;
}

例二

[ZJOI2019] 麻将

假设麻将牌有 nn 种大小的牌,每种大小有 44 张牌.定义面子为三张相邻大小的麻将牌 i,i+1,i+2i,i+1,i+2(顺子)或三种相同大小的麻将牌 i,i,ii,i,i(刻子),对子为两张相同大小的麻将牌 i,ii,i.定义一个麻将牌的序列是胡的,当且仅当它(看作多重集合)可以拆成四个面子和一个对子,或者七个不同的对子.给定 1313 张麻将牌,问期望再摸多少张牌可以满足存在一个胡牌的子序列的条件.

题解

首先,对于一副牌,我们只需考虑每种牌的数量,而不必关心它们的顺序.因此,对于任何一副牌的任意前缀,我们都可以将它转化为一个长度为 nn,每个位置上为 0∼40\sim 4 的序列.初始时,第 ii 张牌的数量为 aia_i,就相当于限制了序列中第 ii 个数字 xix_i 的取值为 [ai,4][a_i,4] 内的整数.注意,转化后的序列没有考虑摸牌的顺序,但是题目中的麻将牌的序列是考虑顺序的.

设 XX 表示可以胡牌的最小摸牌次数.直接计算期望 E[X]\mathbf E[X] 较为困难,可以考虑进行如下转化.设 hih_i 表示摸了 ii 张牌后 没有胡 的序列数.因为它对应的麻将牌序列中,这 ii 张牌必然排在剩下的 (4n−13−i)(4n-13-i) 张牌前方,但是这 ii 张牌和剩下的 (4n−13−i)(4n-13-i) 张牌的顺序是任意的,所以,只摸前 ii 张牌无法胡牌的麻将牌序列的数目就是

hi⋅i!(4n−13−i)!.h_i\cdot i!(4n-13-i)!.

因为麻将牌序列的总数为 (4n−13)!(4n-13)!,所以,只摸前 ii 张牌无法胡牌的概率为

P[X>i]=hi⋅i!(4n−13−i)!(4n−13)!.\mathbf P[X>i] = \dfrac{h_i\cdot i!(4n-13-i)!}{(4n-13)!}.

利用尾求和公式,就可以得到所求期望

E[X]=∑i=0∞P[X>i]=1+∑i=14n−13hi⋅i!(4n−13−i)!(4n−13)!.\mathbf E[X] = \sum_{i=0}^\infty\mathbf P[X>i] = 1 + \sum_{i=1}^{4n-13}\dfrac{h_i\cdot i!(4n-13-i)!}{(4n-13)!}.

至此,问题转化为如何计算 hih_i.我们采用 DP 套 DP 的方法来解决这一问题.

首先,考虑内层 DP,也就是要用动态规划判断一个(转化后的)序列是否对应一副能胡的牌.七对子的情形较为容易,重点讨论第一种胡牌的形式.为此,设 g0/1,i,j,kg_{0/1,i,j,k} 表示处理完前 ii 种牌,还剩 jj 组 (i−1,i)(i−1,i) 以及 kk 张 ii,且存在/不存在对子(即 0/10/1)时最多的面子数.如果对于一个序列进行 DP,最后得到的 g1,ng_{1,n} 中包括一个大于等于 44 的数字,那么这个序列就是能胡的.

这个 DP 的状态转移较为复杂.我们分两步讨论.第一步,考虑 g0/1,ig_{0/1,i} 的转移.这相当于说,如果要向现在的牌型中,添加 xix_i 张大小为 ii 的牌,但是不组成新的对子时,面子数如何转移.显然,如果希望在添加 xix_i 张大小为 ii 的牌后,要得到 ℓ\ell 个顺子,jj 个 (i−1,i)(i-1,i) 和 kk 个单独的 ii,那么,就应该从 (g0/1,i−1)ℓ,j(g_{0/1,i-1})_{\ell,j} 中转移过来(这里的选择最大程度避免了浪费),并将剩余的牌 (xi−ℓ−j−k)(x_i-\ell-j-k) 用于组成尽可能多的刻子.穷举所有的可能性,就得到如下转移方程:

G~(g0/1,i−1,xi)j,k=max⁡{(g0/1,i−1)ℓ,j+ℓ+⌊xi−ℓ−j−k3⌋:ℓ+j+k≤xi}.\tilde G(g_{0/1,i-1}, x_i)_{j,k} = \max\left\{(g_{0/1,i-1})_{\ell,j} + \ell + \left\lfloor\dfrac{x_i-\ell-j-k}{3}\right\rfloor:\ell+j+k\le x_i\right\}.

第二步,再考虑需要组成对子的情形.加入 xix_i 张大小为 ii 的牌时,有如下三种转移:

  • 将 g0,i−1g_{0,i-1} 加 xix_i 张牌转移到 g0,ig_{0,i};
  • 将 g1,i−1g_{1,i-1} 加 xix_i 张牌转移到 g1,ig_{1,i};
  • 若 xi≥2x_i\ge 2,将 g0,i−1g_{0,i-1} 加 xi−2x_i-2 张牌转移到 g1,ig_{1,i}.

由此,就得到了从 gi−1g_{i-1} 添加 xix_i 张牌得到 gig_i 的全部转移.

解决了内层 DP 的状态转移,就可以构建 胡牌自动机.自动机的转移就是上述内层 DP 的转移,还需要考虑如何表示自动机的每个状态.每个状态都对应 gig_i 的一种可能的取值.它有三个维度 (0/1,j,k)(0/1,j,k).因为 jj 和 kk 对应的维度中保留的 (i−1,i)(i-1,i) 和 ii 都是用于将来组成顺子的,而三个相同顺子总是可以重组为三个刻子,所以,只需要考虑组成不超过 22 个相同顺子的需求就行,每种牌型也就只要保留不超过 22 个,即 j,k∈{0,1,2}j,k\in\{0,1,2\}.因此,gig_i 可以表示为一个 2×3×32\times 3\times 3 的数组.另外,为了维护七对子的胡牌牌型,还需要为每个状态添加一个计数器,用于表示当前最多可组成的对子数目.

数组 gig_i 中每个元素的取值范围可能是 {−∞}∪N\{-\infty\}\cup\mathbf N,但是,因为面子数目大于等于 44 都是胡牌,所以,可以限制每个元素取值不超过 44.由于胡牌序列再添加任何牌都是胡牌序列,所以,可以利用 DFA 最小化的思想,将全体胡牌状态压缩为一个状态.因此,对于非胡牌状态,实际上每个位置的取值只要考虑 {−∞}∪{0,1,2,3}\{-\infty\}\cup\{0,1,2,3\} 就可以了.实现时,−∞-\infty 用 −1-1 表示.

尽管如此,所有可能的状态依然相当地多,共有 1+7×5181+7\times 5^{18} 种.穷举它们并不现实.实际上,绝大多数这些可能性都不会真的出现在一个胡牌自动机中.为了避免考虑实际不存在的状态,可以利用 BFS 的思想,从初始状态开始,一步一步扩展状态,直到胡牌状态处停止.这样得到的自动机中有 N=2092N = 2092 个状态.

最后,考虑如何在胡牌自动机上 DP(即外层 DP).设 fi,j,kf_{i,j,k} 表示处理到第 ii 张牌,共摸了 jj 张牌,走到了胡牌自动机上的 kk 号状态的序列数.转移时,枚举摸牌数 0≤t≤4−ai0\leq t\leq 4-a_i,其中 aia_i 为初始 1313 张牌中用掉的 ii 的张数,将之前的序列数乘以 4−ai4−a_i 张牌中选 tt 张牌的方案数 (4−ait)\dbinom{4-a_i}{t},再累加到一起.形式化地,有:

fi+1,j+t,k′=∑t=04−ai(4−ait)fi,j,k.f_{i+1,j+t,k'} = \sum_{t=0}^{4-a_i}\dbinom{4-a_i}{t}f_{i,j,k}.

其中,k′=δ(k,ai+t)k'=\delta(k,a_i+t),表示向自动机的状态 kk 中加入 ai+ta_i+t 张牌后的状态.外层 DP 结束后,就可以计算出摸了 ii 张牌仍然没有胡牌的序列数目,即

hi=∑j=1Nfn,i,j.h_i=\sum_{j=1}^{N} f_{n,i,j}.

代入前文所述表达式,就可以得到所要求的期望.

参考代码
#include <algorithm>
#include <cstdio>
#include <iostream>
#include <map>
#include <unordered_map>
using namespace std;

const int N = 100, M = 400, mod = 998244353;

int n, m, a[N + 5], fact[M + 5], inv[M + 5];

inline int C(int x, int y) {
  return 1LL * fact[x] * inv[y] % mod * inv[x - y] % mod;
}

class HuAutomation {  // 胡牌自动机
 private:
  class Mat {
   private:
    int f[3][3];

   public:
    Mat() {
      for (int i = 0; i <= 2; i++)
        for (int j = 0; j <= 2; j++) f[i][j] = -1;
    }

    int* operator[](const int& x) { return f[x]; }

    inline bool operator==(Mat x) const {
      for (int i = 0; i <= 2; i++)
        for (int j = 0; j <= 2; j++)
          if (f[i][j] != x[i][j]) return 0;
      return 1;
    }

    inline bool operator<(Mat x) const {
      for (int i = 0; i <= 2; i++)
        for (int j = 0; j <= 2; j++)
          if (f[i][j] != x[i][j]) return f[i][j] < x[i][j];
      return 0;
    }

    inline bool Check() {
      for (int i = 0; i <= 2; i++)
        for (int j = 0; j <= 2; j++)
          if (f[i][j] > 3) return 1;
      return 0;
    }

    inline void Upd(Mat x, int t) {
      // 将已有矩阵 x 的值更新到当前矩阵中,模拟加入 t 张当前牌后的状态
      for (int i = 0; i <= 2; i++)
        for (int j = 0; j <= 2; j++)
          if (x[i][j] != -1)
            for (int k = 0; k < 3 && i + j + k <= t; k++)
              f[j][k] = max(f[j][k], min(i + x[i][j] + (t - i - j - k) / 3, 4));
      // i,j,k,(t-i-j-k)
      // 表示分别枚举用于拼面子、用于保留(i-1,i)、用于保留i和直接拼面子的牌数
      //  更新时限制最多 4 个对子
    }
  };

  struct node {
    int t, state[5];
    Mat F[2];

    node() {
      for (int i = 0; i <= 4; i++) state[i] = 0;
      t = 0;
      for (int i = 0; i <= 1; i++) F[i] = Mat();
    }

    inline bool Check() { return t == -1 || t >= 7 || F[1].Check(); }

    node Hu() {
      node x;
      x.t = -1;
      return x;
    }

    bool operator<(const node& x) const {
      if (t == x.t) {
        if (F[0] == x.F[0]) return F[1] < x.F[1];
        return F[0] < x.F[0];
      }
      return t < x.t;
    }

    node operator+(int x) {      // 加上 x 张新牌
      if (Check()) return Hu();  // 胡了直接返回
      node res;
      res.F[0].Upd(F[0], x), res.F[1].Upd(F[1], x);  // 更新
      res.t = t;
      if (x > 1) res.F[1].Upd(F[0], x - 2), ++res.t;
      if (res.Check()) res = Hu();  // 判断是否胡了
      return res;
    }
  } A[2100];

  map<node, int> mp;

  inline int Get(node x) {
    if (mp.count(x)) return mp[x];
    mp[x] = ++tot;
    A[tot] = x;
    return tot;
  }

  inline void Expand(int x) {
    // 构建状态转移图:
    // 对于当前状态 A[x],尝试加入 0~4 张同种牌,得到 5 个新状态
    for (int i = 0; i <= 4; i++) A[x].state[i] = Get(A[x] + i);
  }

  inline node Hu() {
    node x;
    x.t = -1;
    return x;
  }

  inline node Emp() {
    node x;
    x.F[0][0][0] = 0;
    return x;
  }

 public:
  int tot, f[N + 5][M + 5][2100];

  inline void Build() {
    A[1] = Emp(), A[2] = Hu();  // 状态1是初始合法状态,2 是胡牌状态
    mp[A[1]] = 1, mp[A[2]] = 2;
    tot = 2;
    Expand(1);
    for (int i = 3; i <= tot; i++)
      Expand(i);  // 枚举所有可达状态,构建状态图(因为 2
                  // 是胡牌状态所以不需要拓展)
  }

  void DP() {
    f[0][0][1] = 1;
    for (int i = 1; i <= n; i++)
      for (int j = m; j >= 0; j--)
        for (int k = 1; k <= tot; k++)
          if (f[i - 1][j][k])
            for (int t = 0; t <= 4 - a[i]; t++)
              // 枚举这张牌能加多少个(不能超过4,总共已有a[i]个),做组合转移
              (f[i][j + t][A[k].state[a[i] + t]] +=
               1LL * f[i - 1][j][k] * C(4 - a[i], t) % mod) %= mod;
  }
} Hfu;

inline long long qpow(long long x, int y) {
  long long res = 1;
  while (y) {
    if (y & 1) res = res * x % mod;
    x = x * x % mod;
    y >>= 1;
  }
  return res;
}

void init(int n) {
  fact[0] = 1;
  for (int i = 1; i <= n; i++) fact[i] = 1LL * fact[i - 1] * i % mod;
  inv[n] = qpow(fact[n], mod - 2);
  for (int i = n - 1; i >= 0; i--) inv[i] = 1LL * inv[i + 1] * (i + 1) % mod;
}

int main() {
  ios::sync_with_stdio(false);
  cin.tie(nullptr);

  int ans = 0;
  Hfu.Build();  // 状态图构建
  cin >> n;
  for (int i = 1; i <= 13; i++) {
    int x, y;
    cin >> x;
    cin >> y;
    ++a[x];
  }
  m = (n << 2) - 13;
  init(m);
  Hfu.DP();

  for (int i = 1; i <= m; i++) {
    (ans += 1LL * Hfu.f[n][i][1] * fact[i] % mod * fact[m - i] % mod) %= mod;
    for (int j = 3; j <= Hfu.tot; j++)  // 因为 2 节点是胡牌节点所以不统计
      (ans += 1LL * Hfu.f[n][i][j] * fact[i] % mod * fact[m - i] % mod) %= mod;
    // f[i][j][k]:表示处理到第 i 张牌,共摸了 j 张牌,走到了胡牌自动机上的 k
    // 号节点的方案数 每个方案乘以 i!(用于排列这 i 张) × (m - i)!
    // (剩下的牌排列),求和后取模
  }

  cout << (1LL * ans * inv[m] + 1) % mod;
  // 最终答案乘上 inv[m] 是除以 m!,即去掉总排列数,最后 +1 是 dp 的定义是摸完 i
  // 张还没有胡
  return 0;
}

习题

本文内容来自 OI Wiki《DP 套 DP》;作者与修订历史见原文。AIOJ 对排版、站内链接和资源路径做了格式适配。

内容依 CC BY-SA 4.0 与 SATA 许可共享;附加条款以原文声明为准。感谢 OI Wiki 作者与贡献者;请按 SATA 约定前往 OI Wiki 仓库 Star 支持。