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 花神的数论题

比较基础的数位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$部分代码$:

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

这个题有两种思路:

一种是枚举交错和然后求出每种交错和的数的个数,并且还要考虑奇偶性,比较繁琐.

另一种是直接把这些数作为一个整体来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

代码见 题解链接

「BalticOI 2013」非回文数 Palindrome-Free Numbers

这题相比上一题的交错和而言,思维难度有一点提高,而代码难度却大幅降低。

考虑如果存在回文,那么一定存在长度\leq3的回文,所以我们只要记录上一个数字和上上个数字就可以了.

然后同样的,很明显前导0会影响答案,所以记状态的时候再记一维0/1即可.

复杂度O(N^3M+TNM),其中T=1,N=10,M=18

LOJ评测链接

为啥我代码是LOJ上最优解啊(雾,不过70ms一共67个点,有没有人挑战一下卡到每个点只跑1ms啊

一个需要高精度(M=1000)的双倍经验: P3413 SAC#1 - 萌数

LOJ#6274. 数字

考虑数位dp,记dp(n,lx,rx,ly,ry) ( 后4维都是0/1 )表示目前在决策第n位, 目前x/y的值是否正好=上/下界

然后注意一下在转移x and y=0 的位数的时候,要把得到的答案取max.

代码见评测链接

Part.II 对状态进行剪枝或合并

CF55D Beautiful numbers

双倍经验: SP8177 JZPEXT - Beautiful numbers EXTREME但是似乎出锅了不能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题解链接

DP$部分代码$:

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]同类分布

首先枚举余数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$ 部分代码 $:

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]淘金

首先写个爆搜,发现对于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()).

代码见我的题解

Part.III 利用其它方式优化计算

1.算贡献

CF908G New Year and Original Order

这题也有两种思路,不过交错和那道题两种思路不管优劣都能过,但是在这个题目里较劣复杂度的做法就过不去了.

一种思路是对于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), 可以通过本题.

两种思路的代码见我的题解

2.结合位运算/反演等其它计数工具

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

代码:

  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

数位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$部分代码$:

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

可以利用AC自动机,也可以利用SAM来表示状态.

然后数位dp就行了.

dp部分代码:

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和字符串

首先考虑F是怎么计算的.

不难发现我们有一种按位置扫描的暴力:记录当前的答案和当前两个串A,B的1的个数之差.

然后按这个数位DP,再记个两个0/1表示字典序即可.

DP$部分代码$:

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 数页码

P4999 烦人的数学作业

[CQOI2016]手机号码

[SCOI2009] windy 数