#include "aliens.h"
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
namespace NWilberSolver {
template <typename TLookupType>
void ImplementSMAWK(TLookupType lookup,
const std::vector<size_t>& rows,
const std::vector<size_t>& columns,
std::vector<size_t>& result) {
if (rows.empty()) {
return;
}
/// REDUCE
std::vector<size_t> survivingColumns;
for (size_t column: columns) {
while (true) {
if (survivingColumns.empty()) {
break;
}
size_t row = rows[survivingColumns.size() - 1];
const auto& firstValue = lookup(row, column);
const auto& secondValue = lookup(row, survivingColumns.back());
if (firstValue >= secondValue) {
break;
}
survivingColumns.pop_back();
}
if (survivingColumns.size() < rows.size()) {
survivingColumns.push_back(column);
}
}
std::vector<size_t> oddRows;
for (size_t i = 1; i < rows.size(); i += 2) {
oddRows.push_back(rows[i]);
}
ImplementSMAWK(lookup, oddRows, columns, result);
std::unordered_map<size_t, size_t> columnIndices;
for (size_t i = 0; i != survivingColumns.size(); ++i) {
columnIndices[survivingColumns[i]] = i;
}
/// INTERPOLATE
size_t start = 0;
for (size_t r = 0; r < rows.size(); r += 2) {
size_t row = rows[r];
size_t stop = survivingColumns.size() - 1;
size_t argmin = survivingColumns[start];
auto minimum = lookup(row, argmin);
if (r + 1 < rows.size()) {
stop = columnIndices[result[rows[r + 1]]];
}
for (size_t column = start + 1; column <= stop; ++column) {
auto val = lookup(row, survivingColumns[column]);
if (column == start || val < minimum) {
argmin = survivingColumns[column];
minimum = val;
}
}
result[row] = argmin;
start = stop;
}
}
template<typename TLookupType>
std::vector<size_t> SMAWK(int n, int m, TLookupType lookup, bool transpose = false) {
std::vector<size_t> rows(n);
std::vector<size_t> columns(m);
std::iota(rows.begin(), rows.end(), 0);
std::iota(columns.begin(), columns.end(), 0);
if (transpose) {
std::swap(rows, columns);
}
std::vector<size_t> result(rows.size());
auto lookupTranspose = [&](size_t i, size_t j) {
return transpose ? lookup(j, i) : lookup(i, j);
};
ImplementSMAWK(lookupTranspose, rows, columns, result);
return result;
}
using ll = long long;
struct TPairHash {
public:
template <typename T, typename U>
std::size_t operator()(const std::pair<T, U> &x) const {
return std::hash<T>()(x.first) ^ std::hash<U>()(x.second);
}
};
template <typename TPenaltyFunct, typename T = ll>
std::vector<std::pair<size_t, T>> Solve(size_t n, TPenaltyFunct penalty, T init) {
size_t r = 0;
std::vector<std::pair<size_t, T>> D(n + 1);
D[0] = {0, init};
std::unordered_map<std::pair<size_t, size_t>, T, TPairHash> cache;
auto G = [&D, &penalty, &cache](size_t i, size_t j) {
auto it = cache.find(std::make_pair(i, j));
if (it != cache.end()) {
return it->second;
}
auto result = D[i].second + penalty(i, j);
cache[std::make_pair(i, j)] = result;
return result;
};
for (size_t c = 0; c < n;) {
size_t p = std::min(2 * c - r + 1, n);
size_t gRows = c - r + 1;
size_t gCols = p - c;
auto GLookup = [&G, &r, &c, &p](size_t i, size_t j) {
return G(r + i, c + 1 + j);
};
std::vector<size_t> F = SMAWK(gRows, gCols, GLookup, true);
for (size_t it = c + 1; it <= p; ++it) {
D[it] = std::make_pair(r + F[it - c - 1],
GLookup(F[it - c - 1], it - c - 1));
}
size_t bRows = p - c - 1;
size_t bCols = p - c - 1;
auto BLookup = [&G, &r, &c, &p](size_t i, size_t j) {
return G(c + 1 + i, c + 2 + j);
};
std::vector<size_t> H = SMAWK(bRows, bCols, BLookup, true);
size_t j = p + 1;
for (size_t it = c + 2; it <= p; ++it) {
if (BLookup(H[it - c - 2], it - c - 2) <
GLookup(F[it - c - 1], it - c - 1)) {
j = it;
break;
}
}
if (j == p + 1) {
c = p;
} else {
D[j] = std::make_pair(c + 1 + H[j - c - 2],
BLookup(H[j - c - 2], j - c - 2));
r = c + 1;
c = j;
}
}
return D;
}
}
namespace Solver {
const int N = 1e5 + 5;
std::pair<int, int> input[N];
int l[N], r[N];
int n, m, k;
ll dp[N];
int trace[N];
void process() {
std::sort(input, input + n, [](const auto& x, const auto& y) {
return x.first < y.first || (x.first == y.first && x.second > y.second);
});
int rmax = -1;
int cnt = 0;
for (int i = 0; i != n; ++i) {
if (rmax < input[i].second) {
rmax = input[i].second;
++cnt;
l[cnt] = input[i].first;
r[cnt] = input[i].second;
}
}
n = cnt;
}
ll sqr(ll a, bool check = false) {
return check ? sqr(std::max(a, 0ll)) : a * 1ll * a;
}
std::pair<size_t, ll> check(ll lambda) {
auto penalty = [&](int i, int j) {
return lambda +
sqr(r[j] - l[i + 1]) -
sqr(r[i] - l[i + 1], true);
};
auto result = NWilberSolver::Solve(n, penalty, 0ll);
/*std::cerr << "lambda v2 = " << lambda << std::endl;
for (int i = 0; i <= n; ++i) {
std::cerr << "{" << result[i].first << ", " << result[i].second << "} ";
}
std::cerr << std::endl;*/
int layersCount = 0;
for (size_t i = n; i != 0; i = result[i].first) {
++layersCount;
}
return {layersCount, result[n].second};
}
ll solve() {
process();
ll l = 0;
ll r = m * 1ll * m;
while(l < r) {
ll mid = ((l + r) >> 1);
if (check(mid).first <= k) {
r = mid;
} else {
l = mid + 1;
}
}
return check(l).second - l * k;
}
};
ll take_photos(int n, int m, int k, vector<int> r, vector<int> c) {
Solver::n = n;
Solver::m = m;
Solver::k = k;
for (int i = 0; i < n; ++i) {
if (r[i] > c[i]) {
std::swap(r[i], c[i]);
}
Solver::input[i].first = r[i];
Solver::input[i].second = ++c[i];
}
return Solver::solve();
}
