[OI笔记]数位DP合集 & 对数位DP的一点理解
s_r_f
2020-04-28 10:28:25
# $About$
几天前和$pcf$一起口胡了一波典型的数位$DP$题$,$
感觉如果不写代码的话就要变成口胡选手了$,$就来写写$blog$ $+$ 刷一波题。
# 前置芝士：$DP$
本文仅针对对$DP$有一定了解$,$并有一定基础的选手$.$
本文不会过多提及数位$DP$本身$,$主要还是讲应用$.$
# 定义
数位$DP,$就是一类和数位有关的$DP,$一般是计数$/$求和 类型的题$.$
数位$DP$有两种理解方式$:$
### 第一种理解
数位$DP$的状态只和 **已经确定好的数的信息(比如说值/位数/前几个数等)** 和 **没有确定好的数字的位数** 有关$.$
在状态比较多的题中$,$推荐使用记忆化搜索来实现$DP$的过程$.$
数位$DP$分为两个部分$:$ $DP$ 和 对给你的数查询答案 $.$
其中 $,$ 对给你的数查询答案 这一部分的本质就是 $,$ **$O($位数 × 进制 $)$次关于 $DP$ 状态的查询** $,$ 下文中除非特殊情况不会再赘述 $.$
---
### 第二种理解
在状态中计入$0/1$表示是否达到上界$,$就是另一种$DP$方式$,$但是对于每组数据需要重新$DP$一遍$.$
---
**提示:**
下文中复杂度分析中可能存在字母$T,N,M,$
在无特殊说明的情况下$T$为询问组数$,$ $N$为进制 $,$ $M$为最大位数$.$
# $Problems$
## $Part.I$ 基础数位$DP$
---
[洛谷 P4317 花神的数论题](https://www.luogu.com.cn/problem/P4317)
比较基础的数位$DP.$
考虑记$dp(n,m)$表示有$n$个二进制位$,$其中有$m$个$1$的二进制数个数$.$
$($ 实际上这个$dp(n,m)$ 就是 $C_{n}^{m},$但因为组合数学和本文内容没什么关系$,$故省略$.$ $)$
为什么这个题里面没有考虑有没有前导$0$的情况呢$?$因为前导$0$对答案没有影响$.$
然后就可以直接回答$O(log)$次关于$DP$状态的查询了$.$
复杂度$O(M^2+TNM),$其中$T=1,M=60,N=2$
$DP$部分代码$:$
```cpp
LL dp[64][64]; bool vis[64][64];
inline LL DP(int n,int m){
if (n < m || m < 0) return 0; if (n == 0) return 1;
if (vis[n][m]) return dp[n][m]; vis[n][m] = 1;
return dp[n][m] = DP(n-1,m) + DP(n-1,m-1);
}
```
---
[SP1433 KPSUM - The Sum](https://www.luogu.com.cn/problem/SP1433)
这个题有两种思路$:$
一种是枚举交错和然后求出每种交错和的数的个数$,$并且还要考虑奇偶性$,$比较繁琐$.$
另一种是直接把这些数作为一个整体来$dp$交错和$.$
记$dpsum(presum,prelen,n,1/0)$表示已经确定的前缀的交错和为$presum$ $,$已经确定的前缀**排除掉前导零**的长度为$prelen,$还有$n$位没有确定$,$目前的位数中是$/$否有非零数 的交错和$.$
为了$dp$这个数值$,$我们还需要知道$dplen(prelen,n,1/0)$表示已经已经确定的前缀**排除掉前导零**的长度为 $prelen$ $,$ 还有$n$位没有确定$,$目前的位数中是 $/$ 否有非零数 的数字的长度总和$.$
$dp$转移的时候考虑下位置的奇偶性即可$.$
注意到$len/prelen$本身并不重要$,$只和奇偶性有关$,$所以$dp$过程中的$len/prelen$可以对 $2$ 取模 $,$ 可以进一步压缩状态$.$
复杂度$O(N^2M+TNM),$其中$T$的范围没给$,$ $N=10,M=15$
代码见 [题解链接](https://www.luogu.com.cn/blog/s-r-f/solution-sp1433)
---
[「BalticOI 2013」非回文数 Palindrome-Free Numbers ](https://loj.ac/problem/2683)
这题相比上一题的交错和而言$,$思维难度有一点提高$,$~~而代码难度却大幅降低。~~
考虑如果存在回文$,$那么一定存在长度$\leq3$的回文$,$所以我们只要记录上一个数字和上上个数字就可以了$.$
然后同样的$,$很明显前导$0$会影响答案$,$所以记状态的时候再记一维$0/1$即可$.$
复杂度$O(N^3M+TNM),$其中$T=1,N=10,M=18$
[LOJ评测链接](https://loj.ac/submission/797375)
为啥我代码是LOJ上最优解啊(雾,不过70ms一共67个点,有没有人挑战一下卡到每个点只跑1ms啊
一个需要高精度$(M=1000)$的双倍经验$:$ [P3413 SAC#1 - 萌数](https://www.luogu.com.cn/problem/P3413)
----------
[LOJ#6274. 数字 ](https://loj.ac/problem/6274)
考虑数位$dp,$记$dp(n,lx,rx,ly,ry)$ $($ 后4维都是$0/1$ $)$表示目前在决策第$n$位$,$ 目前$x/y$的值是否正好$=$上$/$下界
然后注意一下在转移$x$ $and$ $y=0$ 的位数的时候$,$要把得到的答案取$max.$
代码见[评测链接](https://loj.ac/problem/6274)
--------------------------------
## $Part.II$ 对状态进行剪枝或合并
[CF55D Beautiful numbers](https://www.luogu.com.cn/problem/CF55D)
双倍经验$:$ [SP8177 JZPEXT - Beautiful numbers EXTREME](https://www.luogu.com.cn/problem/SP8177)但是似乎出锅了不能$remote$ $judge?$
首先不难设出状态$:$ $dp(s,r,n)$ 表示目前已经确定下来的前缀$mod$ $2520$ 为 $r,$ 目前的 $1..9$ 是否存在状压到 $s$ 里 $,$ 还有 $n$ 位没有确定$.$
但如果直接设状态状态总数为 $2\times10^7,$再乘上每次转移是 $O(10)$
显然会 $T$ 飞
由于我们只需要考虑集合里面数的$lcm$大小$,$所以有些状态是可以合并的$,$就写个预处理把它合并一下就好了$.$
$P.S:$ 我代码里面没有合并干净$,$是$73$种状态$,$但是也能过$.$
如果真正合并干净了那就是$48$种状态$.$
复杂度$O(V1*V2*NM+TNM),$其中$V1$为$lcm_{i=1}^{10}i=2520,$ $V2$为合并过后的状态总数$.$
[SP8177题解链接](https://www.luogu.com.cn/blog/s-r-f/solution-sp8177)
$DP$部分代码$:$
```cpp
LL pw[19];
int trans[1<<9],tid[1<<9]; bool used[1<<9];
int cnt,v[73+5],nxt[73+5][10];
inline void init(){
int i,j,k,id,l = 1<<9,s;
for (pw[0] = i = 1; i <= 18; ++i) pw[i] = pw[i-1] * 10;
for (i = 0; i < l; ++i){
trans[i] = i;
for (j = 9; j >= 1; --j) if (trans[i] >> (j-1) & 1)
for (k = j + j; k <= 9; k += j) if (trans[i] >> (k-1) & 1){
trans[i] &= ((1<<9)-1)^(1<<(j-1)); break;
}
used[trans[i]] = 1;
}
for (i = 0; i < l; ++i) if (used[i]) ++cnt,v[cnt] = i,tid[i] = cnt;//cnt==73
for (i = 0; i < l; ++i) if (id = tid[i]){
for (v[id] = 1,j = 1; j <= 9; ++j) if (i>>(j-1)&1) v[id] = lcm(v[id],j);
for (nxt[id][0] = id,j = 1; j <= 9; ++j) nxt[id][j] = tid[trans[i|(1<<j-1)]];
}
}
inline bool check(int id,int r){ return r % v[id] == 0; }
LL dp[74][2520][20]; bool vis[74][2520][20];
inline LL DP(int s,int r,int n){
if (!n) return check(s,r) ? 1 : 0;
if (vis[s][r][n]) return dp[s][r][n]; vis[s][r][n] = 1;
LL tot = 0;
for (int i = 0; i <= 9; ++i) tot += DP(nxt[s][i],(r*10+i)%2520,n-1);
return dp[s][r][n] = tot;
}
```
---------------------------
[[AHOI2009]同类分布](https://www.luogu.com.cn/problem/P4127)
首先枚举余数$mod$ $,$然后对于每个余数$DP.$
设$dp(s,r,n)$表示目前已经确定的前缀的数字的总和为$s,$目前已经确定的前缀对$mod$取模之后的结果$.$
然后就可以过了$.$
但是如果直接$dp$的话$,$一个点在$luogu$机器上开$O2$要跑$1.2s,$所以可以剪枝来优化常数$.$
剪枝之后不开$O2$能跑在 $1s$ 以内 $.$
因为$mod$的值域为$O(NM),$而对每个$mod$ $dp$的状态数总数是$O(N^3M^2)$的
所以复杂度$O(N^5M^3+TNM),$其中$N=10,M=18,T=1$
当然因为加了剪枝以及实际上的$mod$只到$M*9=162,$所以效率还是可以的$.$
$DP$ 部分代码 $:$
```cpp
const int S = 162+1;
LL dp[S][S][19]; bool vis[S][S][19]; int M;
inline void Clear(int mod){
register int i,j,k;
M = mod;
for (i = 0; i <= mod; ++i) for (j = 0; j < mod; ++j) for (k = 0; k <= 18; ++k) vis[i][j][k] = dp[i][j][k] = 0;
}
inline LL DP(int s,int r,int n){
if (s < 0 || s > n * 9) return 0; if (!n) return (s == 0 && r == 0) ? 1 : 0;
if (vis[s][r][n]) return dp[s][r][n]; vis[s][r][n] = 1;
LL tot = 0; for (int i = 0; i <= 9; ++i) tot += DP(s-i,(r*10+i)%M,n-1);
return dp[s][r][n] = tot;
}
inline LL solve(LL n){
int mod,i,j,x,s,r; LL tot = 0; bool flag = 0;
for (mod = 1; mod <= 162; ++mod){
Clear(mod); flag = 0;
s = r = 0;
for (i = 18; i >= 0; --i) if ((x = n/pw[i]%10) || flag){
for (j = 0; j < x; ++j) tot += DP(mod-s-j,(r*10+j)%M,i);
s += x,r = (r*10+x) % M; flag = 1;
}
}
return tot;
}
```
---------------------
[ [SDOI2013]淘金](https://www.luogu.com.cn/problem/P3303)
首先写个爆搜$,$发现对于$1\leq i \leq10^{12}$的$i$ $,$ $f(i)$ 只有$8282$ 种 $.$ 所以考虑对于这 $8282$ 种数求出有多少数字满足$f(i) = $这个数字$.$
数位$dp,$ 记 $dp(x,n,1/0)$ 表示目前还有$n$位数字没有决定$,$**除去前导零之后**剩下的数字乘积为$x$ $($ 所有的$x$应当事先搜出来放到数组里 $),$ 目前有$/$无非零数的方案数$.$
求出来之后$,$令$a_i$为第$i$个数字被用到的次数$,$
把$a_i$从大到小排序$,$用堆贪心就可以了$.$
复杂度$O((8282*NM+TNM) \times hash()).$
代码见[我的题解](https://www.luogu.com.cn/blog/s-r-f/solution-p3303)
--------------------------------
## $Part.III$ 利用其它方式优化计算
---------------------------
#### $1.$算贡献
[CF908G New Year and Original Order](https://www.luogu.com.cn/problem/CF908G)
这题也有两种思路$,$不过交错和那道题两种思路不管优劣都能过$,$但是在这个题目里较劣复杂度的做法就过不去了$.$
一种思路是对于$O(700*10)$个询问直接整体$DP$求答案$,$
记$dp[i][j]$表示目前已经安排到了数字$i,$已经放了$j$个数位的方案数和答案$($ 可以用一个$struct$来存 $)$
$.$每次$DP$的复杂度都是$O(700^2*10)$的$,$总复杂度$O(700^3*10^2)$过不去$,$并且不能在只$dp$一次的情况下处理多组询问$.$
这种思路的复杂度是$O(TM^3N^2)$
另一种思路是$,$我们**考虑贡献**$.$
记$sum(n) = \sum\limits_{i=0}^{n-1}10^i$
那么对于一个数字答案显然为 $\sum\limits_{i=0}^{9}$ $sum($ 满足 $k\leq i$ 的数字的个数 $).$
所以对这个 $($ 满足 $k\leq i$ 的数字的个数 $)$ $dp,$就可以降低复杂度$.$
记$dp(n,m,k)$表示目前还有$n$位还没确定$,$目前已经确定下来的前缀中有$m$位数字$\geq$ $k$ $,$直接$dp$即可$.$
处理一组询问的复杂度仍然是$O(10*700)$即$O(NM).$
这种做法的复杂度是 $O(M^2N+TNM),$ 可以通过本题$.$
两种思路的代码见[我的题解](https://www.luogu.com.cn/blog/s-r-f/solution-cf908g)
------------------------
#### $2.$结合位运算/反演等其它计数工具
[P3791 普通数学题](https://www.luogu.com.cn/problem/P3791)
考虑到这个题是$xor,$所以容易想到是进制为$2$的数位$dp.$
不难发现$\sum\limits_{i=1}^{n}d(i)=\sum\limits_{i=1}^n \lfloor \frac{n}{i} \rfloor,$可以$O(\sqrt n)$计算$.$
然后就$O(M^2)$枚举两维的限制$,$再计算即可$.$
怎么计算呢$?$
设两维分别有$a/b$位没有确定$.$
可以发现$,$前面有一些位置的部分是被钦定的$,$而后一部分$($ 后$max(a,b)$位 $)$是没有被确定的$,$且可以任意取$.$
然后答案就是两个$\sum\limits_{i=1}^{n}d(i)$的差$ \times 2^{min(a,b)}$ $.$
不难发现虽然询问有$O(M^2)$个$,$但是本质不同的只有$O(M)$个 $,$ 因为关于$\sum d(i)$的询问只和**未确定数位较多的那个数字的已确定的前缀**有关 $.$
所以复杂度为$O(\sqrt{n} \times M)=O(\sqrt{N^M} \times M),$其中$M=33.$
**注意这个复杂度里的$n$是值域$,$其实际值为$N^M$**
所以这个题就可以拿来加强了$,$比如$d(i$ $xor$ $j$ $xor$ $x$ $)-> $ $φ(i$ $xor$ $j$ $xor$ $x$ $)$或者$d(i$ $xor$ $j$ $xor$ $x$ $)$ $->$ $μ(i$ $xor$ $j$ $xor$ $x$ $)$ /cy
代码$:$
```cpp
const int P = 998244353,MAX = 33;
map<LL,int>FF;
inline int F(LL n){
if (n <= 0) return 0; if (FF.count(n)) return FF[n];
LL sum = 0,l = 1,r;
while (l <= n) r = n/(n/l),sum = (sum + (r-l+1) * (n/l)) % P,l = r+1;
return FF[n] = sum;
}
inline LL get(int l,int r){ return (1ll<<r+1) - (1ll<<l); }
inline int calc(LL x,LL v1,LL v2,int a,int b){
if (a < b) swap(a,b); static LL s;
s = (x^v1^v2)&get(a,MAX);
return (1ll<<b)%P * (F(s+(1ll<<a)-1)-F(s-1)+P) % P;
}
LL n,m,x,ans;
int main(){
int i,j; LL v1,v2;
cin >> n >> m >> x; ++n,++m;
for (i = MAX,v1 = 0; i >= 0; --i) if (n>>i&1){
for (j = MAX,v2 = 0; j >= 0; --j) if (m>>j&1) ans += calc(x,v1,v2,i,j),v2 += 1ll<<j;
v1 += 1ll<<i;
}
cout << (ans%P+P)%P << '\n';
return 0;
}
```
----------------
[CF582D Number of Binominal Coefficients](https://www.luogu.com.cn/problem/CF582D)
数位$DP.$
考虑$lucas$定理的过程$,$设$a = k,b = n-k,$则在$p$进制下$a+b$有$>=k$次进位$,$且$a+b<=A.$
设$dp[n][m][1/0][1/0]$表示当前考虑到第$n$位$,$ 还需要进位$m$次$,$这一位的数值是$/$否紧贴着上界$,$这一位是$/$否往上一位进位$.$
复杂度$O(3500^2),$其中$3500$为$A$在$p$进制下最大位数$.$
这个题目为什么不能用$O(3500*p)$次对$dp$状态的询问呢$?$
因为$p$的值域太大$,$并且这个题没有多组询问$.((($
$DP$部分代码$:$
```cpp
int dp[N][N][2][2];
bool vis[N][N][2][2];
inline LL s(int l,int r){ return l > r ? 0 : (l+r) * (r-l+1ll) / 2 % P; }
inline LL calc(int n){
if (n < 0) return 0; if (n < p) return 1ll * (n+1) * (n+2) / 2 % P;
n = min(n,p*2-2);
return (1ll * (n-p+2) * p % P + s(-p+n+2,p-1)) % P;
}
inline LL calc(int l,int r){ l = max(l,0),r = min(r,p*2-2); return (l > r) ? (0) : (calc(r) - calc(l-1) + P) % P; }
inline LL DP(int n,int m,int s1,int s2){
m = max(m,0);
if (!n) return (!s2 && !m) ? 1 : 0;
if (vis[n][m][s1][s2]) return dp[n][m][s1][s2]; vis[n][m][s1][s2] = 1;
int ans = 0;
if (s1){
if (!s2){
upd(ans,calc(0,p-2) * DP(n-1,m-1,1,1) % P);
upd(ans,calc(0,p-1) * DP(n-1,m,1,0) % P);
}
else{
upd(ans,calc(p-1,p+p-2) * DP(n-1,m-1,1,1) % P);
upd(ans,calc(p,p+p-1) * DP(n-1,m,1,0) % P);
}
}
else{
if (!s2){
upd(ans,calc(0,A[n]-1) * DP(n-1,m,1,0) % P);
upd(ans,calc(A[n],A[n]) * DP(n-1,m,0,0) % P);
upd(ans,calc(0,A[n]-2) * DP(n-1,m-1,1,1) % P);
upd(ans,calc(A[n]-1,A[n]-1) * DP(n-1,m-1,0,1) % P);
}
else{
upd(ans,calc(p,p+A[n]-1) * DP(n-1,m,1,0) % P);
upd(ans,calc(p+A[n],p+A[n]) * DP(n-1,m,0,0) % P);
upd(ans,calc(p-1,p+A[n]-2) * DP(n-1,m-1,1,1) % P);
upd(ans,calc(p+A[n]-1,p+A[n]-1) * DP(n-1,m-1,0,1) % P);
}
}
dp[n][m][s1][s2] = ans;
return ans;
}
```
---
[CF585F Digits of Number Pi](https://www.luogu.com.cn/problem/CF585F)
可以利用$AC$自动机$,$也可以利用$SAM$来表示状态$.$
然后数位$dp$就行了.
dp部分代码$:$
```cpp
pair<int,int> trans[V][D][10];
struct SuffixAutomation{
int cnt,ed; int ch[V][10],fa[V],len[V];
inline void init(){ cnt=ed=1; }
inline void extend(int c){
int p = ed,np = ++cnt; len[np] = len[p] + 1;
while (p && !ch[p][c]) ch[p][c] = np,p = fa[p];
if (!p) fa[np] = 1;
else{
int q = ch[p][c];
if (len[q] == len[p] + 1) fa[np] = q;
else{
int nq = ++cnt;
fa[nq] = fa[q],len[nq] = len[p]+1,memcpy(ch[nq],ch[q],40);
fa[q] = fa[np] = nq;
while (ch[p][c] == q) ch[p][c] = nq,p = fa[p];
}
}
ed = np;
}
inline pair<int,int> getnxt(int now,int l,int c){
while (fa[now] && !ch[now][c]){
now = fa[now]; if (l > len[now]) l = len[now];
}
if (ch[now][c]) now = ch[now][c],++l;
if (l >= (m>>1)) now = l = -1;
return make_pair(now,l);
}
inline void get(){
int i,j,k;
for (i = 1; i <= cnt; ++i) for (j = 0; j <= (m>>1); ++j) for (k = 0; k <= 9; ++k)
trans[i][j][k] = getnxt(i,j,k);
}
}SAM;
int f[2][2][M]; bool visf[2][2][M];
inline int F(int tx,int ty,int i){
if (i > m) return 1;
if (visf[tx][ty][i]) return f[tx][ty][i]; visf[tx][ty][i] = 1;
int r = 0,c;
for (c = 0; c <= 9; ++c){
if (!tx && c < x[i]) continue;
if (!ty && c > y[i]) continue;
upd(r,F(tx|(c>x[i]),ty|(c<y[i]),i+1));
}
return f[tx][ty][i] = r;
}
int g[2][2][M][V][D]; bool visg[2][2][M][V][D];
inline int G(int tx,int ty,int i,int now,int l){
if (i > m) return 0;
if (visg[tx][ty][i][now][l]) return g[tx][ty][i][now][l]; visg[tx][ty][i][now][l] = 1;
int r = 0,c; pair<int,int>t;
for (c = 0; c <= 9; ++c){
if (!tx && c < x[i]) continue;
if (!ty && c > y[i]) continue;
t = trans[now][l][c];
if (t.first != -1) upd(r,G(tx|(c>x[i]),ty|(c<y[i]),i+1,t.first,t.second));
else upd(r,F(tx|(c>x[i]),ty|(c<y[i]),i+1));
}
return g[tx][ty][i][now][l] = r;
}
```
---
4、更加复杂的例子
[牛客挑战赛40C-小V和字符串](https://ac.nowcoder.com/acm/contest/5555/C)
首先考虑$F$是怎么计算的$.$
不难发现我们有一种按位置扫描的暴力$:$记录当前的答案和当前两个串$A,B$的$1$的个数之差$.$
然后按这个数位$DP,$再记个两个$0/1$表示字典序即可$.$
$DP$部分代码$:$
```cpp
int visc[N][2][2][N<<1]; int dpcnt[N][2][2][N<<1];
inline int Cnt(int pos,int bs,int ab,int sta){
if (pos == n+1) return sta == 1002 ? 1 : 0;
if (visc[pos][bs][ab][sta]) return dpcnt[pos][bs][ab][sta]; visc[pos][bs][ab][sta] = 1;
int tot = 0;
for (int va = 0; va <= 1; ++va) for (int vb = 0; vb <= 1; ++vb){
if (va>vb&&!ab) continue;
if (vb>S[pos]-'0' && !bs) continue;
upd(tot,Cnt(pos+1,bs||(vb<(S[pos]-'0')),ab||(va<vb),sta + va - vb));
}
return dpcnt[pos][bs][ab][sta] = tot;
}
int viss[N][2][2][N<<1]; int dps[N][2][2][N<<1];
inline int Sum(int pos,int bs,int ab,int sta){
if (pos == n+1) return 0;
if (viss[pos][bs][ab][sta]) return dps[pos][bs][ab][sta]; viss[pos][bs][ab][sta] = 1;
int tot = 0,vv;
for (int va = 0; va <= 1; ++va) for (int vb = 0; vb <= 1; ++vb){
if (va>vb&&!ab) continue;
if (vb>S[pos]-'0' && !bs) continue;
upd(tot,Sum(pos+1,bs|(vb<S[pos]-'0'),ab|(va<vb),sta + va - vb));
vv = 0;
if (va == vb) vv = 0;
else if (sta == 1002){
if (va == vb) vv = 0; else vv = P-pos;
}
else if (sta < 1002){
if (va) vv = pos; else vv = P-pos;
}
else if (sta > 1002){
if (va) vv = P-pos; else vv = pos;
}
upd(tot,(LL)vv*Cnt(pos+1,bs|(vb<S[pos]-'0'),ab|(va<vb),sta + va - vb)%P);
}
return dps[pos][bs][ab][sta] = tot;
}
```
### 一些例题
[P1836 数页码](https://www.luogu.com.cn/problem/P1836)
[P4999 烦人的数学作业](https://www.luogu.com.cn/problem/P4999)
[[CQOI2016]手机号码](https://www.luogu.com.cn/problem/P4124)
[[SCOI2009] windy 数](https://www.luogu.com.cn/problem/P2657)