最长的 K 序列递增子序列
Longest K Sequential Increasing Subsequences
为什么我创建了一个重复的线程
我在阅读 Longest increasing subsequence with K exceptions allowed. I realised that the person who was asking the question hadn't really understood the problem, because he was referring to a link 后创建了这个线程，它解决了“允许一次更改的最长递增子数组”问题。所以他得到的答案实际上与LIS问题无关。
问题描述
假设给定一个数组A，长度为N。
找到允许 K 个例外的最长递增子序列。
例子
1)
N=9 , K=1
A=[3,9,4,5,8,6,1,3,7]
答案：7
解释：
最长的递增子序列是：3,4,5,8(or 6),1(exception),3,7 -> total=7
- N=11 , K=2
A=[5,6,4,7,3,9,2,5,1,8,7]
答案：8
到目前为止我做了什么...
如果K=1则只允许一个异常。如果使用在O(NlogN)中计算最长递增子序列的已知算法（click here to see this algorithm），那么我们可以计算从A[0]到A[的LIS N-1] 数组 A 的每个元素。我们将结果保存在一个新数组 L 中，大小为 N。查看示例 n.1，L 数组将是：
L=[1,2,2,3,4,4,4,4,5].
使用反向逻辑，我们计算数组R，其中每个元素包含当前从N-1到0的最长递减序列
LIS有一个例外就是sol=max(sol,L[i]+R[i+1]),
其中 sol 被初始化为 sol=L[N-1]。
所以我们从 0 开始计算 LIS，直到索引 i（异常），然后停止并开始一个新的 LIS，直到 N-1.
A=[3,9,4,5,8,6,1,3,7]
L=[1,2,2,3,4,4,4,4,5]
R=[5,4,4,3,3,3,3,2,1]
Sol = 7
-> 分步说明：
init: sol = L[N]= 5
i=0 : sol = max(sol,1+4) = 5
i=1 : sol = max(sol,2+4) = 6
i=2 : sol = max(sol,2+3) = 6
i=3 : sol = max(sol,3+3) = 6
i=4 : sol = max(sol,4+3) = 7
i=4 : sol = max(sol,4+3) = 7
i=4 : sol = max(sol,4+2) = 7
i=5 : sol = max(sol,4+1) = 7
复杂度：
O( NlogN + NlogN + N ) = O(NlogN)
因为数组 R, L 需要 NlogN 时间来计算，我们还需要 Θ(N) 才能找到 sol.
k=1 问题的代码
#include <stdio.h>
#include <vector>
std::vector<int> ends;
int index_search(int value, int asc) {
int l = -1;
int r = ends.size() - 1;
while (r - l > 1) {
int m = (r + l) / 2;
if (asc && ends[m] >= value)
r = m;
else if (asc && ends[m] < value)
l = m;
else if (!asc && ends[m] <= value)
r = m;
else
l = m;
}
return r;
}
int main(void) {
int n, *S, *A, *B, i, length, idx, max;
scanf("%d",&n);
S = new int[n];
L = new int[n];
R = new int[n];
for (i=0; i<n; i++) {
scanf("%d",&S[i]);
}
ends.push_back(S[0]);
length = 1;
L[0] = length;
for (i=1; i<n; i++) {
if (S[i] < ends[0]) {
ends[0] = S[i];
}
else if (S[i] > ends[length-1]) {
length++;
ends.push_back(S[i]);
}
else {
idx = index_search(S[i],1);
ends[idx] = S[i];
}
L[i] = length;
}
ends.clear();
ends.push_back(S[n-1]);
length = 1;
R[n-1] = length;
for (i=n-2; i>=0; i--) {
if (S[i] > ends[0]) {
ends[0] = S[i];
}
else if (S[i] < ends[length-1]) {
length++;
ends.push_back(S[i]);
}
else {
idx = index_search(S[i],0);
ends[idx] = S[i];
}
R[i] = length;
}
max = A[n-1];
for (i=0; i<n-1; i++) {
max = std::max(max,(L[i]+R[i+1]));
}
printf("%d\n",max);
return 0;
}
泛化到 K 个异常
我已经提供了 K=1 的算法。我不知道如何更改上述算法以适用于 K 异常。如果有人能帮助我，我会很高兴。
此答案已从 my answer 修改为 Computer Science Stackexchange 上的类似问题。
最多 k 个异常的 LIS 问题允许使用拉格朗日松弛的 O(n log² n) 算法。当 k 大于 log n 时，这在 O(nk log n) DP 上渐近改进，我们也将简要解释。
令DP[a][b]表示最长递增子序列的长度，最多有b个例外（前一个整数大于下一个整数的位置）结束于元素 b 一个。这个DP在算法中没有涉及，但是定义它使得算法的证明更容易。
为方便起见，我们假设所有元素都是不同的，并且数组中的最后一个元素是它的最大值。请注意，这并不限制我们，因为我们只需将 m / 2n 添加到每个数字的第 m 次出现，并将无穷大添加到数组并从答案中减去 1。设 V 为 1 <= V[i] <= n 是第 i 个元素的值的排列。
为了解决复杂度为 O(nk log n) 的问题，我们保持不变量，即 DP[a][b] 已针对 b < j 计算。从 0 到 k 循环 j，在第 j 次迭代计算所有 a 的 DP[a][j]。为此，从 1 到 n 循环 i。我们维护 x < i 上的 DP[x][j-1] 的最大值和一个前缀最大数据结构，在索引 i 处，对于 x < i，在位置 V[x] 处具有 DP[x][j]，并且 0在其他位置。
我们有 DP[i][j] = 1 + max(DP[i'][j], DP[x][j-1]) 我们遍历 i', x < i, V [i'] < V[i]。 DP[x][j-1] 的前缀最大值为我们提供了第二类项的最大值，查询前缀最大值数据结构中的前缀 [0, V[i]] 为我们提供了第一类项的最大值类型。然后更新prefix maximum和prefix maximum数据结构。
这是该算法的 C++ 实现。请注意，此实现不假定数组的最后一个元素是其最大值，或者数组不包含重复项。
#include <iostream>
#include <vector>
#include <algorithm>
using namespace std;
// Fenwick tree for prefix maximum queries
class Fenwick {
private:
vector<int> val;
public:
Fenwick(int n) : val(n+1, 0) {}
// Sets value at position i to maximum of its current value and
void inc(int i, int v) {
for (++i; i < val.size(); i += i & -i) val[i] = max(val[i], v);
}
// Calculates prefix maximum up to index i
int get(int i) {
int res = 0;
for (++i; i > 0; i -= i & -i) res = max(res, val[i]);
return res;
}
};
// Binary searches index of v from sorted vector
int bins(const vector<int>& vec, int v) {
int low = 0;
int high = (int)vec.size() - 1;
while(low != high) {
int mid = (low + high) / 2;
if (vec[mid] < v) low = mid + 1;
else high = mid;
}
return low;
}
// Compresses the range of values to [0, m), and returns m
int compress(vector<int>& vec) {
vector<int> ord = vec;
sort(ord.begin(), ord.end());
ord.erase(unique(ord.begin(), ord.end()), ord.end());
for (int& v : vec) v = bins(ord, v);
return ord.size();
}
// Returns length of longest strictly increasing subsequence with at most k exceptions
int lisExc(int k, vector<int> vec) {
int n = vec.size();
int m = compress(vec);
vector<int> dp(n, 0);
for (int j = 0;; ++j) {
Fenwick fenw(m+1); // longest subsequence with at most j exceptions ending at this value
int max_exc = 0; // longest subsequence with at most j-1 exceptions ending before this
for (int i = 0; i < n; ++i) {
int off = 1 + max(max_exc, fenw.get(vec[i]));
max_exc = max(max_exc, dp[i]);
dp[i] = off;
fenw.inc(vec[i]+1, off);
}
if (j == k) return fenw.get(m);
}
}
int main() {
int n, k;
cin >> n >> k;
vector<int> vec(n);
for (int i = 0; i < n; ++i) cin >> vec[i];
int res = lisExc(k, vec);
cout << res << '\n';
}
现在我们将 return 到 O(n log² n) 算法。 Select 一些整数 0 <= r <= n。定义DP'[a][r] = max(DP[a][b] - rb)，其中最大值接管b，MAXB[a][r]为最大b使得DP'[a][ r] = DP[a][b] - rb，而 MINB[a][r] 类似于最小这样的 b。我们将证明 DP[a][k] = DP'[a][r] + rk 当且仅当 MINB[a][r] <= k <= MAXB[a][r]。此外，我们将证明对于任何 k 存在一个 r，该不等式成立。
注意 MINB[a][r] >= MINB[a][r'] 和 MAXB[a][r] >= MAXB[a][r'] 如果 r < r'，因此如果我们假设两个声明的结果，我们可以对 r 进行二进制搜索，尝试 O(log n) 值。因此，如果我们可以在 O(n log n) 时间内计算 DP'、MINB 和 MAXB，那么我们的复杂度为 O(n log² n)。
为此，我们需要一个存储元组 P[i] = (v_i, low_i, high_i) 的线段树，并支持以下操作：
给定一个范围[a,b]，找出该范围内的最大值（最大值v_i，a <= i <= b），以及最小低点和最大高点与范围内的那个值配对。
设置元组P[i]的值
这很容易实现，假设对线段树有一定的了解，每个操作的复杂度为 O(log n) 次。具体可以参考下面算法的实现。
我们现在将展示如何在 O(n log n) 中计算 DP'、MINB 和 MAXB。修复 r。构建最初包含 n+1 个空值（-INF、INF、-INF）的线段树。我们认为对于小于当前位置 i 的 j，P[V[j]] = (DP'[j], MINB[j], MAXB[j])。如果 r > 0，则将 DP'[0] = 0、MINB[0] = 0 和 MAXB[0] 设置为 0，否则设置为 INF 和 P[0] = (DP'[0], MINB[0], MAXB[ 0]).
循环 i 从 1 到 n。有两种类型的以 i 结尾的子序列：一种是前一个元素大于 V[i]，另一种是小于 V[i]。要考虑第二种情况，请查询 [0, V[i]] 范围内的线段树。让结果为 (v_1, low_1, high_1)。设置 off1 = (v_1 + 1, low_1, high_1)。对于第一种，查询[V[i], n]范围内的线段树。让结果为 (v_2, low_2, high_2)。设置 off2 = (v_2 + 1 - r, low_2 + 1, high_2 + 1)，其中我们因创建异常而受到 r 的惩罚。
然后我们把off1和off2合并成off。如果 off1.v > off2.v 设置 off = off1，如果 off2.v > off1.v 设置 off = off2。否则，设置 off = (off1.v, min(off1.low, off2.low), max(off1.high, off2.high))。然后设置DP'[i] = off.v, MINB[i] = off.low, MAXB[i] = off.high and P[i] = off.
由于我们在每个 i 处进行两次线段树查询，因此总共需要 O(n log n) 时间。很容易通过归纳证明我们计算出正确的值 DP'、MINB 和 MAXB。
所以简而言之，算法是：
预处理，修改值，使它们形成排列，最后的值是最大值。
二进制搜索正确的 r，初始边界 0 <= r <= n
用空值初始化线段树，设置DP'[0]、MINB[0]和MAXB[0]。
从 i = 1 到 n 的循环，在第 i 步
- 查询线段树的范围[0, V[i]]和[V[i], n],
- 根据这些查询计算 DP'[i]、MINB[i] 和 MAXB[i]，并且
- 将线段树中V[i]位置的值设置为元组(DP'[i], MINB[i], MAXB[i]).
如果 MINB[n][r] <= k <= MAXB[n][r], return DP'[n][r] + kr - 1.
否则，如果MAXB[n][r] < k，则正确的r小于当前的r。如果 MINB[n][r] > k，则正确的 r 大于当前的 r。将 r 和 return 上的边界更新为第 1 步。
这是此算法的 C++ 实现。它还会找到最优子序列。
#include <iostream>
#include <vector>
#include <algorithm>
using namespace std;
using ll = long long;
const int INF = 2 * (int)1e9;
pair<ll, pair<int, int>> combine(pair<ll, pair<int, int>> le, pair<ll, pair<int, int>> ri) {
if (le.first < ri.first) swap(le, ri);
if (ri.first == le.first) {
le.second.first = min(le.second.first, ri.second.first);
le.second.second = max(le.second.second, ri.second.second);
}
return le;
}
// Specialised range maximum segment tree
class SegTree {
private:
vector<pair<ll, pair<int, int>>> seg;
int h = 1;
pair<ll, pair<int, int>> recGet(int a, int b, int i, int le, int ri) const {
if (ri <= a || b <= le) return {-INF, {INF, -INF}};
else if (a <= le && ri <= b) return seg[i];
else return combine(recGet(a, b, 2*i, le, (le+ri)/2), recGet(a, b, 2*i+1, (le+ri)/2, ri));
}
public:
SegTree(int n) {
while(h < n) h *= 2;
seg.resize(2*h, {-INF, {INF, -INF}});
}
void set(int i, pair<ll, pair<int, int>> off) {
seg[i+h] = combine(seg[i+h], off);
for (i += h; i > 1; i /= 2) seg[i/2] = combine(seg[i], seg[i^1]);
}
pair<ll, pair<int, int>> get(int a, int b) const {
return recGet(a, b+1, 1, 0, h);
}
};
// Binary searches index of v from sorted vector
int bins(const vector<int>& vec, int v) {
int low = 0;
int high = (int)vec.size() - 1;
while(low != high) {
int mid = (low + high) / 2;
if (vec[mid] < v) low = mid + 1;
else high = mid;
}
return low;
}
// Finds longest strictly increasing subsequence with at most k exceptions in O(n log^2 n)
vector<int> lisExc(int k, vector<int> vec) {
// Compress values
vector<int> ord = vec;
sort(ord.begin(), ord.end());
ord.erase(unique(ord.begin(), ord.end()), ord.end());
for (auto& v : vec) v = bins(ord, v) + 1;
// Binary search lambda
int n = vec.size();
int m = ord.size() + 1;
int lambda_0 = 0;
int lambda_1 = n;
while(true) {
int lambda = (lambda_0 + lambda_1) / 2;
SegTree seg(m);
if (lambda > 0) seg.set(0, {0, {0, 0}});
else seg.set(0, {0, {0, INF}});
// Calculate DP
vector<pair<ll, pair<int, int>>> dp(n);
for (int i = 0; i < n; ++i) {
auto off0 = seg.get(0, vec[i]-1); // previous < this
off0.first += 1;
auto off1 = seg.get(vec[i], m-1); // previous >= this
off1.first += 1 - lambda;
off1.second.first += 1;
off1.second.second += 1;
dp[i] = combine(off0, off1);
seg.set(vec[i], dp[i]);
}
// Is min_b <= k <= max_b?
auto off = seg.get(0, m-1);
if (off.second.second < k) {
lambda_1 = lambda - 1;
} else if (off.second.first > k) {
lambda_0 = lambda + 1;
} else {
// Construct solution
ll r = off.first + 1;
int v = m;
int b = k;
vector<int> res;
for (int i = n-1; i >= 0; --i) {
if (vec[i] < v) {
if (r == dp[i].first + 1 && dp[i].second.first <= b && b <= dp[i].second.second) {
res.push_back(i);
r -= 1;
v = vec[i];
}
} else {
if (r == dp[i].first + 1 - lambda && dp[i].second.first <= b-1 && b-1 <= dp[i].second.second) {
res.push_back(i);
r -= 1 - lambda;
v = vec[i];
--b;
}
}
}
reverse(res.begin(), res.end());
return res;
}
}
}
int main() {
int n, k;
cin >> n >> k;
vector<int> vec(n);
for (int i = 0; i < n; ++i) cin >> vec[i];
vector<int> ans = lisExc(k, vec);
for (auto i : ans) cout << i+1 << ' ';
cout << '\n';
}
我们现在来证明这两个说法。我们想证明
DP'[a][r] = DP[a][b] - rb 当且仅当 MINB[a][r] <= b <= MAXB[a][r ]
对于所有a,k存在一个整数r,0 <= r <= n,使得MINB[a][r] <= k <= MAXB[a][r]
这两者都是由问题的凹性得出的。凹性意味着 DP[a][k+2] - DP[a][k+1] <= DP[a][k+1] - DP[a][k] 对于所有 a, k。这很直观：我们被允许的例外越多，允许的例外越少对我们有帮助。
修复a和r。设 f(b) = DP[a][b] - rb，d(b) = f(b+1) - f(b)。根据问题的凹性，我们有 d(k+1) <= d(k)。假设对于所有 i，x < y 且 f(x) = f(y) >= f(i)。因此 d(x) <= 0，因此对于 [x, y) 中的 i，d(i) <= 0。但是 f(y) = f(x) + d(x) + d(x + 1) + ... + d(y - 1)，因此对于 [x, y) 中的 i，d(i) = 0。因此 f(y) = f(x) = f(i) for i in [x, y]。这证明了第一个说法。
为了证明第二个，设r = DP[a][k+1] - DP[a][k] 并如前定义f, d。然后 d(k) = 0，因此对于 i < k，d(i) >= 0，对于 i > k，d(i) <= 0，因此 f(k) 是所需的最大值。
证明凹性更难。有关证明，请参阅 my answer 处的 cs.stackexchange。
为什么我创建了一个重复的线程
我在阅读 Longest increasing subsequence with K exceptions allowed. I realised that the person who was asking the question hadn't really understood the problem, because he was referring to a link 后创建了这个线程，它解决了“允许一次更改的最长递增子数组”问题。所以他得到的答案实际上与LIS问题无关。
问题描述
假设给定一个数组A，长度为N。 找到允许 K 个例外的最长递增子序列。
例子
1)
N=9 , K=1
A=[3,9,4,5,8,6,1,3,7]
答案：7
解释：
最长的递增子序列是：3,4,5,8(or 6),1(exception),3,7 -> total=7
- N=11 , K=2
A=[5,6,4,7,3,9,2,5,1,8,7]
答案：8
到目前为止我做了什么...
如果K=1则只允许一个异常。如果使用在O(NlogN)中计算最长递增子序列的已知算法（click here to see this algorithm），那么我们可以计算从A[0]到A[的LIS N-1] 数组 A 的每个元素。我们将结果保存在一个新数组 L 中，大小为 N。查看示例 n.1，L 数组将是： L=[1,2,2,3,4,4,4,4,5].
使用反向逻辑，我们计算数组R，其中每个元素包含当前从N-1到0的最长递减序列
LIS有一个例外就是sol=max(sol,L[i]+R[i+1]), 其中 sol 被初始化为 sol=L[N-1]。 所以我们从 0 开始计算 LIS，直到索引 i（异常），然后停止并开始一个新的 LIS，直到 N-1.
A=[3,9,4,5,8,6,1,3,7]
L=[1,2,2,3,4,4,4,4,5]
R=[5,4,4,3,3,3,3,2,1]
Sol = 7
-> 分步说明：
init: sol = L[N]= 5
i=0 : sol = max(sol,1+4) = 5
i=1 : sol = max(sol,2+4) = 6
i=2 : sol = max(sol,2+3) = 6
i=3 : sol = max(sol,3+3) = 6
i=4 : sol = max(sol,4+3) = 7
i=4 : sol = max(sol,4+3) = 7
i=4 : sol = max(sol,4+2) = 7
i=5 : sol = max(sol,4+1) = 7
复杂度： O( NlogN + NlogN + N ) = O(NlogN)
因为数组 R, L 需要 NlogN 时间来计算，我们还需要 Θ(N) 才能找到 sol.
k=1 问题的代码
#include <stdio.h>
#include <vector>
std::vector<int> ends;
int index_search(int value, int asc) {
int l = -1;
int r = ends.size() - 1;
while (r - l > 1) {
int m = (r + l) / 2;
if (asc && ends[m] >= value)
r = m;
else if (asc && ends[m] < value)
l = m;
else if (!asc && ends[m] <= value)
r = m;
else
l = m;
}
return r;
}
int main(void) {
int n, *S, *A, *B, i, length, idx, max;
scanf("%d",&n);
S = new int[n];
L = new int[n];
R = new int[n];
for (i=0; i<n; i++) {
scanf("%d",&S[i]);
}
ends.push_back(S[0]);
length = 1;
L[0] = length;
for (i=1; i<n; i++) {
if (S[i] < ends[0]) {
ends[0] = S[i];
}
else if (S[i] > ends[length-1]) {
length++;
ends.push_back(S[i]);
}
else {
idx = index_search(S[i],1);
ends[idx] = S[i];
}
L[i] = length;
}
ends.clear();
ends.push_back(S[n-1]);
length = 1;
R[n-1] = length;
for (i=n-2; i>=0; i--) {
if (S[i] > ends[0]) {
ends[0] = S[i];
}
else if (S[i] < ends[length-1]) {
length++;
ends.push_back(S[i]);
}
else {
idx = index_search(S[i],0);
ends[idx] = S[i];
}
R[i] = length;
}
max = A[n-1];
for (i=0; i<n-1; i++) {
max = std::max(max,(L[i]+R[i+1]));
}
printf("%d\n",max);
return 0;
}
泛化到 K 个异常
我已经提供了 K=1 的算法。我不知道如何更改上述算法以适用于 K 异常。如果有人能帮助我，我会很高兴。
此答案已从 my answer 修改为 Computer Science Stackexchange 上的类似问题。
最多 k 个异常的 LIS 问题允许使用拉格朗日松弛的 O(n log² n) 算法。当 k 大于 log n 时，这在 O(nk log n) DP 上渐近改进，我们也将简要解释。
令DP[a][b]表示最长递增子序列的长度，最多有b个例外（前一个整数大于下一个整数的位置）结束于元素 b 一个。这个DP在算法中没有涉及，但是定义它使得算法的证明更容易。
为方便起见，我们假设所有元素都是不同的，并且数组中的最后一个元素是它的最大值。请注意，这并不限制我们，因为我们只需将 m / 2n 添加到每个数字的第 m 次出现，并将无穷大添加到数组并从答案中减去 1。设 V 为 1 <= V[i] <= n 是第 i 个元素的值的排列。
为了解决复杂度为 O(nk log n) 的问题，我们保持不变量，即 DP[a][b] 已针对 b < j 计算。从 0 到 k 循环 j，在第 j 次迭代计算所有 a 的 DP[a][j]。为此，从 1 到 n 循环 i。我们维护 x < i 上的 DP[x][j-1] 的最大值和一个前缀最大数据结构，在索引 i 处，对于 x < i，在位置 V[x] 处具有 DP[x][j]，并且 0在其他位置。
我们有 DP[i][j] = 1 + max(DP[i'][j], DP[x][j-1]) 我们遍历 i', x < i, V [i'] < V[i]。 DP[x][j-1] 的前缀最大值为我们提供了第二类项的最大值，查询前缀最大值数据结构中的前缀 [0, V[i]] 为我们提供了第一类项的最大值类型。然后更新prefix maximum和prefix maximum数据结构。
这是该算法的 C++ 实现。请注意，此实现不假定数组的最后一个元素是其最大值，或者数组不包含重复项。
#include <iostream>
#include <vector>
#include <algorithm>
using namespace std;
// Fenwick tree for prefix maximum queries
class Fenwick {
private:
vector<int> val;
public:
Fenwick(int n) : val(n+1, 0) {}
// Sets value at position i to maximum of its current value and
void inc(int i, int v) {
for (++i; i < val.size(); i += i & -i) val[i] = max(val[i], v);
}
// Calculates prefix maximum up to index i
int get(int i) {
int res = 0;
for (++i; i > 0; i -= i & -i) res = max(res, val[i]);
return res;
}
};
// Binary searches index of v from sorted vector
int bins(const vector<int>& vec, int v) {
int low = 0;
int high = (int)vec.size() - 1;
while(low != high) {
int mid = (low + high) / 2;
if (vec[mid] < v) low = mid + 1;
else high = mid;
}
return low;
}
// Compresses the range of values to [0, m), and returns m
int compress(vector<int>& vec) {
vector<int> ord = vec;
sort(ord.begin(), ord.end());
ord.erase(unique(ord.begin(), ord.end()), ord.end());
for (int& v : vec) v = bins(ord, v);
return ord.size();
}
// Returns length of longest strictly increasing subsequence with at most k exceptions
int lisExc(int k, vector<int> vec) {
int n = vec.size();
int m = compress(vec);
vector<int> dp(n, 0);
for (int j = 0;; ++j) {
Fenwick fenw(m+1); // longest subsequence with at most j exceptions ending at this value
int max_exc = 0; // longest subsequence with at most j-1 exceptions ending before this
for (int i = 0; i < n; ++i) {
int off = 1 + max(max_exc, fenw.get(vec[i]));
max_exc = max(max_exc, dp[i]);
dp[i] = off;
fenw.inc(vec[i]+1, off);
}
if (j == k) return fenw.get(m);
}
}
int main() {
int n, k;
cin >> n >> k;
vector<int> vec(n);
for (int i = 0; i < n; ++i) cin >> vec[i];
int res = lisExc(k, vec);
cout << res << '\n';
}
现在我们将 return 到 O(n log² n) 算法。 Select 一些整数 0 <= r <= n。定义DP'[a][r] = max(DP[a][b] - rb)，其中最大值接管b，MAXB[a][r]为最大b使得DP'[a][ r] = DP[a][b] - rb，而 MINB[a][r] 类似于最小这样的 b。我们将证明 DP[a][k] = DP'[a][r] + rk 当且仅当 MINB[a][r] <= k <= MAXB[a][r]。此外，我们将证明对于任何 k 存在一个 r，该不等式成立。
注意 MINB[a][r] >= MINB[a][r'] 和 MAXB[a][r] >= MAXB[a][r'] 如果 r < r'，因此如果我们假设两个声明的结果，我们可以对 r 进行二进制搜索，尝试 O(log n) 值。因此，如果我们可以在 O(n log n) 时间内计算 DP'、MINB 和 MAXB，那么我们的复杂度为 O(n log² n)。
为此，我们需要一个存储元组 P[i] = (v_i, low_i, high_i) 的线段树，并支持以下操作：
给定一个范围[a,b]，找出该范围内的最大值（最大值v_i，a <= i <= b），以及最小低点和最大高点与范围内的那个值配对。
设置元组P[i]的值
这很容易实现，假设对线段树有一定的了解，每个操作的复杂度为 O(log n) 次。具体可以参考下面算法的实现。
我们现在将展示如何在 O(n log n) 中计算 DP'、MINB 和 MAXB。修复 r。构建最初包含 n+1 个空值（-INF、INF、-INF）的线段树。我们认为对于小于当前位置 i 的 j，P[V[j]] = (DP'[j], MINB[j], MAXB[j])。如果 r > 0，则将 DP'[0] = 0、MINB[0] = 0 和 MAXB[0] 设置为 0，否则设置为 INF 和 P[0] = (DP'[0], MINB[0], MAXB[ 0]).
循环 i 从 1 到 n。有两种类型的以 i 结尾的子序列：一种是前一个元素大于 V[i]，另一种是小于 V[i]。要考虑第二种情况，请查询 [0, V[i]] 范围内的线段树。让结果为 (v_1, low_1, high_1)。设置 off1 = (v_1 + 1, low_1, high_1)。对于第一种，查询[V[i], n]范围内的线段树。让结果为 (v_2, low_2, high_2)。设置 off2 = (v_2 + 1 - r, low_2 + 1, high_2 + 1)，其中我们因创建异常而受到 r 的惩罚。
然后我们把off1和off2合并成off。如果 off1.v > off2.v 设置 off = off1，如果 off2.v > off1.v 设置 off = off2。否则，设置 off = (off1.v, min(off1.low, off2.low), max(off1.high, off2.high))。然后设置DP'[i] = off.v, MINB[i] = off.low, MAXB[i] = off.high and P[i] = off.
由于我们在每个 i 处进行两次线段树查询，因此总共需要 O(n log n) 时间。很容易通过归纳证明我们计算出正确的值 DP'、MINB 和 MAXB。
所以简而言之，算法是：
预处理，修改值，使它们形成排列，最后的值是最大值。
二进制搜索正确的 r，初始边界 0 <= r <= n
用空值初始化线段树，设置DP'[0]、MINB[0]和MAXB[0]。
从 i = 1 到 n 的循环，在第 i 步
- 查询线段树的范围[0, V[i]]和[V[i], n],
- 根据这些查询计算 DP'[i]、MINB[i] 和 MAXB[i]，并且
- 将线段树中V[i]位置的值设置为元组(DP'[i], MINB[i], MAXB[i]).
如果 MINB[n][r] <= k <= MAXB[n][r], return DP'[n][r] + kr - 1.
否则，如果MAXB[n][r] < k，则正确的r小于当前的r。如果 MINB[n][r] > k，则正确的 r 大于当前的 r。将 r 和 return 上的边界更新为第 1 步。
这是此算法的 C++ 实现。它还会找到最优子序列。
#include <iostream>
#include <vector>
#include <algorithm>
using namespace std;
using ll = long long;
const int INF = 2 * (int)1e9;
pair<ll, pair<int, int>> combine(pair<ll, pair<int, int>> le, pair<ll, pair<int, int>> ri) {
if (le.first < ri.first) swap(le, ri);
if (ri.first == le.first) {
le.second.first = min(le.second.first, ri.second.first);
le.second.second = max(le.second.second, ri.second.second);
}
return le;
}
// Specialised range maximum segment tree
class SegTree {
private:
vector<pair<ll, pair<int, int>>> seg;
int h = 1;
pair<ll, pair<int, int>> recGet(int a, int b, int i, int le, int ri) const {
if (ri <= a || b <= le) return {-INF, {INF, -INF}};
else if (a <= le && ri <= b) return seg[i];
else return combine(recGet(a, b, 2*i, le, (le+ri)/2), recGet(a, b, 2*i+1, (le+ri)/2, ri));
}
public:
SegTree(int n) {
while(h < n) h *= 2;
seg.resize(2*h, {-INF, {INF, -INF}});
}
void set(int i, pair<ll, pair<int, int>> off) {
seg[i+h] = combine(seg[i+h], off);
for (i += h; i > 1; i /= 2) seg[i/2] = combine(seg[i], seg[i^1]);
}
pair<ll, pair<int, int>> get(int a, int b) const {
return recGet(a, b+1, 1, 0, h);
}
};
// Binary searches index of v from sorted vector
int bins(const vector<int>& vec, int v) {
int low = 0;
int high = (int)vec.size() - 1;
while(low != high) {
int mid = (low + high) / 2;
if (vec[mid] < v) low = mid + 1;
else high = mid;
}
return low;
}
// Finds longest strictly increasing subsequence with at most k exceptions in O(n log^2 n)
vector<int> lisExc(int k, vector<int> vec) {
// Compress values
vector<int> ord = vec;
sort(ord.begin(), ord.end());
ord.erase(unique(ord.begin(), ord.end()), ord.end());
for (auto& v : vec) v = bins(ord, v) + 1;
// Binary search lambda
int n = vec.size();
int m = ord.size() + 1;
int lambda_0 = 0;
int lambda_1 = n;
while(true) {
int lambda = (lambda_0 + lambda_1) / 2;
SegTree seg(m);
if (lambda > 0) seg.set(0, {0, {0, 0}});
else seg.set(0, {0, {0, INF}});
// Calculate DP
vector<pair<ll, pair<int, int>>> dp(n);
for (int i = 0; i < n; ++i) {
auto off0 = seg.get(0, vec[i]-1); // previous < this
off0.first += 1;
auto off1 = seg.get(vec[i], m-1); // previous >= this
off1.first += 1 - lambda;
off1.second.first += 1;
off1.second.second += 1;
dp[i] = combine(off0, off1);
seg.set(vec[i], dp[i]);
}
// Is min_b <= k <= max_b?
auto off = seg.get(0, m-1);
if (off.second.second < k) {
lambda_1 = lambda - 1;
} else if (off.second.first > k) {
lambda_0 = lambda + 1;
} else {
// Construct solution
ll r = off.first + 1;
int v = m;
int b = k;
vector<int> res;
for (int i = n-1; i >= 0; --i) {
if (vec[i] < v) {
if (r == dp[i].first + 1 && dp[i].second.first <= b && b <= dp[i].second.second) {
res.push_back(i);
r -= 1;
v = vec[i];
}
} else {
if (r == dp[i].first + 1 - lambda && dp[i].second.first <= b-1 && b-1 <= dp[i].second.second) {
res.push_back(i);
r -= 1 - lambda;
v = vec[i];
--b;
}
}
}
reverse(res.begin(), res.end());
return res;
}
}
}
int main() {
int n, k;
cin >> n >> k;
vector<int> vec(n);
for (int i = 0; i < n; ++i) cin >> vec[i];
vector<int> ans = lisExc(k, vec);
for (auto i : ans) cout << i+1 << ' ';
cout << '\n';
}
我们现在来证明这两个说法。我们想证明
DP'[a][r] = DP[a][b] - rb 当且仅当 MINB[a][r] <= b <= MAXB[a][r ]
对于所有a,k存在一个整数r,0 <= r <= n,使得MINB[a][r] <= k <= MAXB[a][r]
这两者都是由问题的凹性得出的。凹性意味着 DP[a][k+2] - DP[a][k+1] <= DP[a][k+1] - DP[a][k] 对于所有 a, k。这很直观：我们被允许的例外越多，允许的例外越少对我们有帮助。
修复a和r。设 f(b) = DP[a][b] - rb，d(b) = f(b+1) - f(b)。根据问题的凹性，我们有 d(k+1) <= d(k)。假设对于所有 i，x < y 且 f(x) = f(y) >= f(i)。因此 d(x) <= 0，因此对于 [x, y) 中的 i，d(i) <= 0。但是 f(y) = f(x) + d(x) + d(x + 1) + ... + d(y - 1)，因此对于 [x, y) 中的 i，d(i) = 0。因此 f(y) = f(x) = f(i) for i in [x, y]。这证明了第一个说法。
为了证明第二个，设r = DP[a][k+1] - DP[a][k] 并如前定义f, d。然后 d(k) = 0，因此对于 i < k，d(i) >= 0，对于 i > k，d(i) <= 0，因此 f(k) 是所需的最大值。
证明凹性更难。有关证明，请参阅 my answer 处的 cs.stackexchange。