Skip to content

DP of DP

Introduction

This article introduces the concept of DP over DP and demonstrates its application through two example problems.

Concept

The so-called "DP over DP" actually refers to a technique where, during the dynamic programming process, the solution process of a subproblem (usually a DP) is abstracted into an automaton (DFA), and then a new layer of DP is designed based on this automaton.

This technique is mainly applied to a class of sequence counting, probability, or expectation problems. A typical problem has the following structure:

  • Given a character set \(\Sigma\) and a set \(A\subseteq\Sigma^n\) of "valid sequences" of length \(n\) over it. Depending on the character set, sequences can be binary strings, digit strings, state sequences, etc.
  • For each specific sequence \(s\in\Sigma^n\), a dynamic programming algorithm can be used to determine whether it is valid (i.e., \(s\in A\)), compute its weight, or find related values.
  • Ultimately, we want to count the number of all sequences in set \(A\), their total weight, expected value, etc.

In such cases, enumerating all sequences is infeasible. Therefore, we consider abstracting the process of "determining whether a sequence is valid" (i.e., the inner DP) into a Deterministic Finite Automaton (DFA). Generally, for a fixed sequence \(s\in\Sigma^n\), the state function of the inner DP can be expressed as \(g(i,x;s)\), which represents the value of some quantity after processing the prefix of length \(i\) of sequence \(s\), with other state components being \(x\). Correspondingly, the state transition equation of the inner DP is

\[ g(i,\cdot;s) = G(g(i-1,\cdot;s),s_i). \]

That is, the function \(g(i,\cdot;s)\) is uniquely determined by the previous function \(g(i-1,\cdot;s)\) and the current character \(s_i\). If we view the function \(g(i,\cdot;s)\) as a state of the automaton, then the state transition equation of the inner DP gives a transition of the automaton. Therefore, the structure of the automaton \((Q,\Sigma,\delta,q_0,F)\) corresponding to the inner DP is as follows:

  • The state set \(Q\) is the set of all possible functions \(g(i,\cdot;s)\) for \(s\in\Sigma^n\) and \(i=0,1,\cdots,n\);
  • The transition function \(\delta:Q\times\Sigma\to Q\) is the \(G\) in the inner DP's state transition equation;
  • The start state \(q_0\) is usually obvious, i.e., the initial state of the inner DP;
  • The set of accepting states \(F\) corresponds to all valid sequences \(s\in A\).

The function \(g(i,\cdot;s)\) itself can be quite complex. Therefore, when dealing with specific problems, it is often necessary to perform state compression or combine DFA minimization techniques to compress the state space. This is the main reason why DP over DP can significantly reduce time and space complexity compared to brute-force DP.

After abstracting the inner DP into a DFA, a new DP can be designed on this DFA to solve the original problem, i.e., the outer DP. For ease of explanation, take a simple counting problem as an example. The state function of the outer DP is defined as \(f(i,q)\), which represents the number of prefixes of length \(i\) that reach state \(q\in Q\) in the DFA. Its state transition equation is

\[ f(i,q) = \sum_{c\in\Sigma}\sum_{q'\in Q:\delta(q',c)=q} f(i-1,q'). \]

The start state is naturally \(f(0,q_0)\), and the final answer is usually computed simply from \(\{f(n,q):q\in F\}\). The outer DP is actually a special case of DP on DAGs.

Examples

The next two examples will detail the general approach of DP over DP.

Example 1

Hero meet devil

Given a string \(S\) over the character set ACGT, with \(|S|\le 15\). For each \(0\leq i \leq |S|\), find the number of strings \(T\) of length \(m\) over the character set ACGT such that the length of their longest common subsequence with \(S\) is \(i\).

Solution

We first think of a DP: let \(f_{i,j}\) denote the number of strings \(T\) of length \(i\) whose longest common subsequence with \(S\) has length \(j\). However, this cannot be transferred because the main issue is that we don't know which characters correspond to this longest common subsequence.

Consider the naive process of computing the longest common subsequence. Let \(g_{i,j}\) denote the length of the longest common subsequence between the first \(i\) characters of \(T\) and the first \(j\) characters of \(S\). Then we have

\[ g_{i,j} = \max\{g_{i-1,j},g_{i,j-1},g_{i-1,j-1}+[T_i=S_j]\}. \]

We find that for a fixed \(i\), we only need to record the values of the one-dimensional array \(g_i\) to accurately maintain the state of the longest common subsequence between \(S\) and the first \(i\) characters of \(T\). Since the length of \(S\) is only \(15\), we find this idea is feasible.

So we redefine the state \(f_{i,x}\) as the number of strings \(T\) of length \(i\) such that the DP array (i.e., \(g_i\)) with \(S\) has state \(x\). This DP seems to have many states, but we observe that \(g_{i,j}-g_{i,j-1}\in\{0,1\}\), so we can maintain the difference array of \(g_i\), and the number of states is \(2^{|S|}\).

Now think about how to transition. It is easy to see that if we know the array \(g_i\) and \(T_{i+1}\), we can compute \(g_{i+1}\) via the naive LCS transition (i.e., the DP equation above). Thus, the naive LCS becomes the inner DP that helps \(f\) transition.

Therefore, we enumerate \(T_{i+1}\), compute the state \(x'\) after transitioning from \(x\), and then add \(f_{i,x}\) to \(f_{i+1,x'}\) to complete the state transition of the outer DP. Finally, we record \(\textit{ans}_i\) as the answer for LCS length \(i\), enumerate each state \(S\), and add \(f_{m,S}\) to \(\textit{ans}_{\operatorname{popcount}(S)}\).

Reference Code 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
#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];  // State transition table: nx[mask][ch] represents the new state after appending ch to mask state
int f[1001]
     [1 << 16];  // DP array: f[i][mask] represents the number of schemes with length i and state mask
int ans[MAXN +
        1];  // Final answer: ans[i] represents the number of schemes containing i matching positions when length is m

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

int calc(int mask, int x) {  // Given current mask state and current character
                             // x (0~3), returns the new state
  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) {  // Count how many 1s are in mask (i.e., how many positions are matched)
  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;
}

Example 2

[ZJOI2019] Mahjong

Assume there are \(n\) types of mahjong tiles, each type has \(4\) tiles. Define a meld as three consecutive tiles \(i,i+1,i+2\) (a sequence) or three identical tiles \(i,i,i\) (a triplet). Define a pair as two identical tiles \(i,i\). A sequence of mahjong tiles is said to be winning if and only if it (viewed as a multiset) can be partitioned into four melds and one pair, or seven distinct pairs. Given \(13\) tiles, ask the expected number of additional draws needed to satisfy the condition that there exists a winning subsequence.

Solution

First, for any hand of tiles, we only need to consider the count of each tile type, not their order. Therefore, for any prefix of any hand, we can transform it into a sequence of length \(n\), where each position is an integer from \(0\) to \(4\). Initially, the count of the \(i\)-th tile is \(a_i\), which is equivalent to restricting the \(i\)-th number \(x_i\) in the sequence to be an integer in \([a_i,4]\). Note that the transformed sequence does not consider the order of draws, but the sequence of mahjong tiles in the problem does consider order.

Let \(X\) be the minimum number of draws to achieve a winning hand. Directly computing the expectation \(\mathbf E[X]\) is difficult, so consider the following transformation. Let \(h_i\) denote the number of sequences of \(i\) draws that are not winning. In the corresponding mahjong tile sequence, these \(i\) tiles must come before the remaining \((4n-13-i)\) tiles, but the order among these \(i\) tiles and the remaining \((4n-13-i)\) tiles is arbitrary. Therefore, the number of mahjong tile sequences where only the first \(i\) draws are not winning is

\[ h_i\cdot i!(4n-13-i)!. \]

Since the total number of mahjong tile sequences is \((4n-13)!\), the probability that only the first \(i\) draws are not winning is

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

Using the tail sum formula, the desired expectation is

\[ \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)!}. \]

Thus, the problem reduces to computing \(h_i\). We use DP over DP to solve this problem.

First, consider the inner DP, which is to use dynamic programming to determine whether a (transformed) sequence corresponds to a winning hand. The case of seven pairs is easier; we focus on the first winning form. For this, let \(g_{0/1,i,j,k}\) denote the maximum number of melds after processing the first \(i\) tile types, with \(j\) remaining groups of \((i−1,i)\) and \(k\) remaining tiles of type \(i\), and with/without a pair (i.e., \(0/1\)). If after DP on a sequence, the resulting \(g_{1,n}\) contains a number greater than or equal to \(4\), then the sequence is winning.

The state transition for this DP is relatively complex. We discuss it in two steps. First, consider the transition of \(g_{0/1,i}\). This is equivalent to adding \(x_i\) tiles of type \(i\) to the current hand without forming a new pair, and seeing how the meld count transitions. Clearly, if after adding \(x_i\) tiles of type \(i\), we want \(\ell\) sequences, \(j\) groups of \((i-1,i)\), and \(k\) single \(i\) tiles, then we should transition from \((g_{0/1,i-1})_{\ell,j}\) (this choice minimizes waste) and use the remaining tiles \((x_i-\ell-j-k)\) to form as many triplets as possible. Enumerating all possibilities gives the following transition equation:

\[ \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\}. \]

Second, consider the case where a pair needs to be formed. When adding \(x_i\) tiles of type \(i\), there are three types of transitions:

  • Transition from \(g_{0,i-1}\) with \(x_i\) tiles to \(g_{0,i}\);
  • Transition from \(g_{1,i-1}\) with \(x_i\) tiles to \(g_{1,i}\);
  • If \(x_i\ge 2\), transition from \(g_{0,i-1}\) with \(x_i-2\) tiles to \(g_{1,i}\).

This gives all transitions from \(g_{i-1}\) to \(g_i\) by adding \(x_i\) tiles.

Having solved the state transition of the inner DP, we can construct the winning automaton. The transitions of the automaton are the above inner DP transitions. We also need to consider how to represent each state of the automaton. Each state corresponds to a possible value of \(g_i\). It has three dimensions \((0/1,j,k)\). Since the \((i-1,i)\) groups and \(i\) tiles retained in the \(j\) and \(k\) dimensions are for forming sequences in the future, and three identical sequences can always be recombined into three triplets, we only need to consider the need to form at most \(2\) identical sequences. Thus, for each tile type, we only need to retain at most \(2\), i.e., \(j,k\in\{0,1,2\}\). Therefore, \(g_i\) can be represented as a \(2\times 3\times 3\) array. Additionally, to handle the seven-pairs winning hand, we need to add a counter to each state to indicate the current maximum number of pairs that can be formed.

The range of each element in the array \(g_i\) could be \(\{-\infty\}\cup\mathbf N\), but since a meld count greater than or equal to \(4\) is winning, we can limit each element to at most \(4\). Since any winning sequence remains winning after adding any tiles, we can use DFA minimization ideas to compress all winning states into a single state. Therefore, for non-winning states, we only need to consider values in \(\{-\infty\}\cup\{0,1,2,3\}\). In implementation, \(-\infty\) can be represented by \(-1\).

Even so, the number of possible states is still quite large, totaling \(1+7\times 5^{18}\). Enumerating them is not practical. In fact, most of these possibilities never actually appear in a winning automaton. To avoid considering states that don't actually exist, we can use a BFS approach, starting from the initial state and expanding step by step until reaching winning states. The resulting automaton has \(N = 2092\) states.

Finally, consider performing DP on the winning automaton (i.e., the outer DP). Let \(f_{i,j,k}\) denote the number of sequences after processing the \(i\)-th tile type, having drawn \(j\) tiles in total, and being at state \(k\) in the winning automaton. During transition, enumerate the number of draws \(0\leq t\leq 4-a_i\), where \(a_i\) is the number of tiles of type \(i\) used in the initial \(13\) tiles. Multiply the previous sequence count by the number of ways to choose \(t\) tiles from the remaining \(4−a_i\) tiles, i.e., \(\dbinom{4-a_i}{t}\), and sum them up. Formally:

\[ f_{i+1,j+t,k'} = \sum_{t=0}^{4-a_i}\dbinom{4-a_i}{t}f_{i,j,k}. \]

Here, \(k'=\delta(k,a_i+t)\) denotes the state after adding \(a_i+t\) tiles to automaton state \(k\). After completing the outer DP, we can compute the number of sequences that are still not winning after drawing \(i\) tiles:

\[ h_i=\sum_{j=1}^{N} f_{n,i,j}. \]

Substituting into the expression given earlier yields the desired expectation.

Reference Code 

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
```cpp
#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 {  // Hu (Winning) Automation
 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) {
      // Update the current matrix with the values from existing matrix x,
      // simulating the state after adding t copies of the current tile.
      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)
      // Represent enumerating the number of tiles used for forming melds,
      // for keeping (i-1, i), for keeping i, and for directly forming melds, respectively.
      // Limit to at most 4 pairs during update.
    }
  };

  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) {      // Add x new tiles
      if (Check()) return Hu();  // If already winning, return winning state
      node res;
      res.F[0].Upd(F[0], x), res.F[1].Upd(F[1], x);  // Update
      res.t = t;
      if (x > 1) res.F[1].Upd(F[0], x - 2), ++res.t;
      if (res.Check()) res = Hu();  // Check if winning
      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) {
    // Build state transition graph:
    // For current state A[x], try adding 0~4 copies of the same tile, get 5 new states.
    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();  // State 1 is initial valid state, 2 is winning state
    mp[A[1]] = 1, mp[A[2]] = 2;
    tot = 2;
    Expand(1);
    for (int i = 3; i <= tot; i++)
      Expand(i);  // Enumerate all reachable states, build state graph (since state 2
                  // is winning, it doesn't need expansion)
  }

  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++)
              // Enumerate how many of this tile can be added (cannot exceed 4, already have a[i]),
              // perform combinatorial transition.
              (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();  // Build state graph
  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++)  // Node 2 is the winning node, so exclude it from statistics
      (ans += 1LL * Hfu.f[n][i][j] * fact[i] % mod * fact[m - i] % mod) %= mod;
    // f[i][j][k]: represents the number of ways after processing the i-th tile,
    // having drawn j tiles in total, and being at node k on the Hu automation.
    // Multiply each way by i! (for arranging these i tiles) × (m - i)!
    // (arranging the remaining tiles), sum up and take modulo.
  }

  cout << (1LL * ans * inv[m] + 1) % mod;
  // Multiply the final answer by inv[m] to divide by m!, i.e., remove the total permutations.
  // Finally +1 because dp is defined as not having won after drawing i tiles.
  return 0;
}
```

Exercises

本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用