In mathematics you don't understand things, you just get used to them.

动态 书籍 杂项 关于 友链

标签 dp 下的文章

重新学了一遍记忆化搜索实现的数位 dp，谈一下理解，毕竟我是普及组小朋友。

这个东西的实际意义大概就是通过将「有上界限制的分支使用搜索实现，一般的分支使用 dp 转移」的手法，将带有上界限制的数位计数问题一般化为所有数位都可以任意取的计数问题。

拜托，存数位请一定不要将 num[1] 置为最高位！比如 $1000$、$100$，两个数在指针相等时所指的数位不一样。果然还是要看 DJ 的！

link。

首先考虑暴力，枚举规划前缀 $[1, i]$ 和前缀 mex $x$，则我们需要 $x$ 个数来填了 $[0, x)$，还剩下 $i-x$ 个数随便填 $[0, x) \cup (x, n]$，如果直接组合数算可能会出问题，考虑 dp。

定义 $f[i][x][j]$ 表示规划前缀 $[1, i]$，当前的 mex 为 $x$，还有 $j$ 个数当前不对 mex 的取值造成影响（也就是说他们都大于 $x$，这 $j$ 个数是指在 $a$ 数组中的，所以我们不必关心这 $j$ 个数具体是什么）。转移就分两种情况：

  • $a[i+1] \neq x$：$(i+1, x, j) \gets (i+1, x, j)+(i, x, j)*(x+j)$ 和 $(i+1, x, j+1) \gets (i, x, j)$，第一个就是考虑当加入的 $a[i+1]$ 属于那 $j$ 个数或者属于 $[0, x)$，第二个很简单；
  • $a[i+1] = x$：设当前 mex 变成了 $y$，则有 $(i+1, y, j-y+x+1) \gets (i+1, y, j-y+x+1)+(i, x, j) \times \binom{j}{y-x-1} \times (y-x-1)!$，意义明显，注意后面那个是排列数而不是组合数。

然后这个是 $O(n^2k^2)$ 的，把刷表改成填表后前缀和优化即可。

using modint = modint998244353;
int n, K, a[2100];
modint dp[2][2100][2100], sum[2][2100][2100], fct[2100], ifct[2100];
inline int L(int i) { return max(0, a[i]-K); }
inline int R(int i) { return min(i, a[i]+K); }
signed main() {
    ios::sync_with_stdio(0);
    cin.tie(0);
    fct[0] = 1;
    for (int i=1; i<2100; ++i) {
        fct[i] = fct[i-1]*i;
    }
    ifct[2099] = fct[2099].inv();
    for (int i=2098; i>=0; --i) {
        ifct[i] = ifct[i+1]*(i+1);
    }
    cin >> n >> K;
    for (int i=1; i<=n; ++i) {
        cin >> a[i];
    }
    dp[0][0][0] = sum[0][0][0] = 1;
    for (int i=1,cur=1; i<=n; ++i,cur^=1) {
        for (int j=0; j<=i; ++j) {
            for (int x=L(i); x<=R(i) && x<=j; ++x) {
                dp[cur][x][j] = dp[cur^1][x][j]*j;
                if (j) {
                    dp[cur][x][j] += dp[cur^1][x][j-1];
                }
                if (j && x) {
                    dp[cur][x][j] += sum[cur^1][min(x-1, R(i-1))][j-1]*ifct[j-x];
                }
                sum[cur][x][j] = dp[cur][x][j]*fct[j-x];
                if (x) {
                    sum[cur][x][j] += sum[cur][x-1][j];
                }
            }
        }
        for (int j=0; j<=i-1; ++j) {
            for (int x=L(i-1); x<=R(i-1) && x<=j; ++x) {
                dp[cur^1][x][j] = sum[cur^1][x][j] = 0;
            }
        }
    }
    modint ans = 0;
    for (int i=0; i<=n; ++i) {
        for (int j=L(n); j<=R(n) && j<=i; ++j) {
            ans += dp[n&1][j][i]*fct[n-j]*ifct[n-i];
        }
    }
    cout << ans.val() << "\n";
}

link。

一个 dp（拜谢 ly）和切入点都略有不同的做法，并不需要观察啥性质。

原问题针对子序列进行规划，自然地想到转而对前缀进行规划。接下来我们考虑一个前缀 $[1, i]$ 以及一个 $j \in [1, i]$ 对答案的贡献，可以写出 $\displaystyle \textit{cont}(j \text{ to } [1, i]) = [\max_{i < k} a_k > a_j] \times \text{the number of LISs containing } j \text{ indexed in } [1, i]$。

这个可以做两个 dp 解决，第一个好解决的静态 dp，即结束在 $j$ 的 LIS 方案数；第二个 dp 有些烦：$j$ 在动，我们考虑的前缀 $[1, i]$ 也在移动。

到这里其实换一下着手点第二个 dp 就变成静态的了，即考虑位置 $j$，直接算 $(j, i)$ 的贡献即可，其中 $i$ 是「满足 $a_i > a_j$ 的最大的 $i$」。于是我们的第二个 dp 就可以被描述为从 $j$ 开始，结束在 $i$ 之前（不包括，因为要保证 $\max_{i < k} a_k > a_j$）的 LIS 方案数。答案即 $\displaystyle \sum_{i=1}^n dp_i \times dp'_i$。

第二个 dp 的具体实现，还需要一个辅助的 dp，其定义是第二个 dp 的定义中去掉贡献范围上界（即 $i$），转移画画图就能理解了。

用下 Fenwick Tree 啥的就能 $O(n \log_2 n)$ 了。

using modint = modint1000000007;
int n, a[200100], pos[200100], id[200100];
modint dp[200100], dp2[200100], dp3[200100], bit[200100], bit2[200100];
void cng(int x, modint w) {
    for (; x<=n; x+=x&-x) {
        bit[x] += w;
    }
}
modint qry(int x) {
    modint res = 0;
    for (; x; x-=x&-x) {
        res += bit[x];
    }
    return res;
}
void cng2(int x, modint w) {
    for (; x<=n; x+=x&-x) {
        bit2[x] += w;
    }
}
modint qry2(int x) {
    modint res = 0;
    for (; x; x-=x&-x) {
        res += bit2[x];
    }
    return res;
}
void solve() {
    cin >> n;
    bastr<int> dsc;
    for (int i=1; i<=n; ++i) {
        cin >> a[i];
        dsc += a[i];
    }
    sort(dsc.begin(), dsc.end());
    dsc.erase(unique(dsc.begin(), dsc.end()), dsc.end());
    for (int i=1; i<=n; ++i) {
        a[i] = lower_bound(dsc.begin(), dsc.end(), a[i])-dsc.begin()+1;
    }
    iota(id+1, id+n+1, 1);
    sort(id+1, id+n+1, [&](int x, int y) {
        return a[x] < a[y];
    });
    for (int i=1,now=n; i<=n; ++i) {
        while (now >= 1 && a[now] <= a[id[i]]) {
            now--;
        }
        pos[id[i]] = now;
    }
    for (int i=1; i<=n; ++i) {
        bit[i] = 0;
    }
    for (int i=1; i<=n; ++i) {
        dp[i] = qry(a[i]-1)+1;
        cng(a[i], dp[i]);
    }
    for (int i=1; i<=n; ++i) {
        bit[i] = bit2[i] = 0;
    }
    modint ans = 0;
    for (int i=n; i>=1; --i) {
        dp2[i] = qry(cs(dsc)-a[i])+1;
        cng(cs(dsc)-a[i]+1, dp2[i]);
        if (pos[i] > i) {
            dp3[i] = qry(cs(dsc)-a[pos[i]]+1)+qry2(a[pos[i]]-1)-qry2(a[i]);
            cng2(a[i], dp3[i]);
        }
        else {
            dp3[i] = 0;
        }
        ans += dp[i]*dp3[i];
    }
    cout << ans.val() << "\n";
}

link。

容易发现，如果将 $x$ 写作 $\displaystyle \prod_{i = 1}^k p_i^{\alpha_i}$ 的形式，$\displaystyle J(x) = 1+\sum p_i^{\alpha_i}+\sum\sum p_i^{\alpha_i}p_j^{\alpha_j}+\dots = \sum_{T \in 2^S} \sum_{i \in T} p_i^{\alpha_i}$，其中 $S = \{1, \dots, k\}$。写到这里可以发现 Joker function 是个 multiplicative function，所以 Joker function 又可以写作 $\displaystyle J(x) = \prod_{i = 1}^k J(p_i^{\alpha_i}) = \prod_{i = 1}^k \left(p_i^{\alpha_i}+1\right)$。

考虑 dp，设 $f[i][j]$ 表示用前 $i$ 个因数通过上述形式凑出 $j$ 的方案数，转移很平凡，具体看代码。这个 dp 可以用数据结构来存。这个 dp 复杂度的保证在于 $\displaystyle \max_{i \in [1, 10^{12}]} {\sigma_0(i)} = 6720$，这就叫人比较无语。至于判断一个数是否是 $p^k+1$ 的形式，可以根号分治来做，记阈值为 $T = 10^6$，对于 $\leqslant T$ 的元素，可以线性筛出所有质数，标记所有质数的次幂，数量级粗略是 $O(T\log T)$，完全跑不满；对于 $> T$ 的元素，因为 $(T+1)^2 > 10^{12}$，所以只需要判断是否为质数即可，Miller Rabin 或者其他 nb 的判素数方法即可。写不来 MR 所以直接蒯了个 mrsrz 的。

namespace Miller_Rabin{
    const int P[]={2,3,7,97,19260817};
    inline LL mul(LL a,LL b,const LL&md){
        LL tmp=a*b-(LL)((long double)a*b/md+.5)*md;
        return tmp+(tmp>>63&md);
    } 
    inline LL pow(LL a,LL b,const LL&md){
        LL ret=1;
        for(;b;b>>=1,a=mul(a,a,md))if(b&1)ret=mul(ret,a,md);
        return ret;
    }
    bool test(LL a,LL p){
        LL x=p-1;
        int d=0;
        while(!(x&1))++d,x>>=1;
        LL t=pow(a,x,p);
        while(d--){
            const LL ls=t;
            t=mul(t,t,p);
            if(t==1&&ls!=1&&ls!=p-1)return 0;
        }
        return t==1;
    }
    bool check(LL n){
        if(n<2)return 0;
        if(n==2||n==3||n==7||n==97||n==P[4])return 1;
        for(int i=0;i<5;++i)if(n%P[i]==0)return 0;
        for(int i=0;i<5;++i)
        if(!test(P[i]%n,n))return 0;
        return 1;
    }
}
using Miller_Rabin::check;
LL A;
bitset<1000100> tag;
int prm[1000100], cnt;
map<LL, int> mp; // a in it iff a = p^k
bastr<LL> offset[2000100];
void sieve(int n) {
    for (int i=2; i<=n; ++i) {
        if (!tag.test(i)) {
            prm[++cnt] = i;
        }
        for (int j=1; j<=cnt && i<=n/prm[j]; ++j) {
            tag.set(i*prm[j]);
            if (i%prm[j] == 0) {
                break;
            }
        }
    }
    for (int i=1; i<=cnt; ++i) {
        for (LL j=prm[i]; j<=1e12; j*=prm[i]) {
            mp[j] = i;
        }
    }
}
void executer() {
    cin >> A;
    map<LL, LL> dp[2];
    bastr<int> eff;
    for (LL i=1; i<=A/i; ++i) {
        if (A%i == 0) {
            if (mp.count(i-1)) {
                offset[mp[i-1]] += i;
                eff += mp[i-1];
            }
            if (i*i != A) {
                if (mp.count(A/i-1)) {
                    offset[mp[A/i-1]] += A/i;
                    eff += mp[A/i-1];
                }
                else if (check(A/i-1)) {
                    offset[++cnt] += A/i;
                    eff += cnt;
                }
            }
        }
    }
    sort(eff.begin(), eff.end());
    eff.erase(unique(eff.begin(), eff.end()), eff.end());
    dp[0][1] = 1;
    int cur = 1;
    for (auto u : eff) {
        dp[cur] = dp[cur^1];
        for (auto v : dp[cur^1]) {
            for (auto p : offset[u]) {
                if (A/p >= v.first && A%(v.first*p) == 0) {
                    dp[cur][v.first*p] += v.second;
                }
            }
        }
        cur ^= 1;
    }
    cout << dp[cs(eff)&1][A] << "\n";
}

link。

如果做过 codeforces - 1144G 那这题最多 *2200。

序列中的最大值必然为其中一个拐点，不妨设 $a_p = a_\max$，先讨论另一个拐点 $i$ 在 $p$ 左侧的情况。于是问题转化为规划 $[1, i)$，$(i, p)$，$(p, n]$ 几个区间中的数给 $i$ 还是给 $p$。

  • 对于 $[1, i)$，令 $dp[i]$ 表示将 $[1, i]$ 分为两个 strictly increasing subsequence，其中不以 $i$ 结尾的 subsequence 的结尾元素的最小值，分讨转移即可。
  • 对于 $(i, p)$，同 codeforces - 1144G，which is 这题唯一的难点。
  • 对于 $(p, n]$，同情况 1。

不太会写代码，，，参考了下 editorial，，，

#include <bits/stdc++.h>
#define cm(x, y) x = min(x, y)
#define cm2(x, y) x = max(x, y)
using namespace std;
int n, a[500100], pos, dp[500100], dp2[500100], dp3[500100][2];
int solve() {
    pos = n;
    for (int i=0; i<n; ++i) {
        if (a[i] > a[pos]) {
            pos = i;
        }
    }
    dp[0] = -1;
    for (int i=1; i<=pos; ++i) {
        dp[i] = 1e9;
        if (a[i] > dp[i-1]) {
            cm(dp[i], a[i-1]);
        }
        if (a[i] > a[i-1]) {
            cm(dp[i], dp[i-1]);
        }
    }
    dp2[n-1] = -1;
    for (int i = n-2; i>=pos; --i) {
        dp2[i] = 1e9;
        if (a[i] > dp2[i+1]) {
            cm(dp2[i], a[i+1]);
        }
        if (a[i] > a[i+1]) {
            cm(dp2[i], dp2[i+1]);
        }
    }
    int res = 0;
    dp3[pos][0] = dp[pos];
    for (int i=pos+1; i<n; ++i) {
        dp3[i][0] = 1e9;
        dp3[i][1] = -1e9;
        if (a[i-1] > a[i]) {
            cm(dp3[i][0], dp3[i-1][0]);
        }
        if (dp3[i-1][1] > a[i]) {
            cm(dp3[i][0], a[i-1]);
        }
        if (a[i-1] < a[i]) {
            cm2(dp3[i][1], dp3[i-1][1]);
        }
        if (dp3[i-1][0] < a[i]) {
            cm2(dp3[i][1], a[i-1]);
        }
        res += dp3[i][1] > dp2[i];
    }
    return res;
}
signed main() {
    ios::sync_with_stdio(0);
    cin.tie(0);
    cin >> n;
    for (int i=0; i<n; ++i) {
        cin >> a[i];
    }
    int ret = solve();
    reverse(a, a+n);
    cout << ret+solve() << "\n";
}
  1. 1
  2. ...
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. ...
  11. 15