PROBLEM LINK:
Author: Jatin Yadav
Tester: Radoslav Dimitrov
Editorialist: William Lin
DIFFICULTY:
Hard
PREREQUISITES:
Geometry, Linear Programming (thanks to @anand20)
PROBLEM:
There are N convex polygons. Moving a polygon to another point has a cost of the Manhattan distance. Find the minimum cost to move all the polygons such that there exists a permutation P where polygon P_i is inside P_{i+1}.
QUICK EXPLANATION:
First, sort the polygons by their size. Then, we can formulate this as a linear programming problem: Our variables will be the amount we move each polygon in the four directions and our constraints will make sure polygon i is inside polygon i+1.
EXPLANATION:
Observation 1. We should sort the polygons by their size, since it’s impossible to fit a big polygon into a smaller one.
Let r_i, l_i, u_i, and d_i be the distance that we move polygon i to the right, left, up, and down, respectively.
Observation 2. The condition that polygon i has to be contained in polygon i+1 can be rewritten as a combination of linear inequalities.
Note that a convex polygon is the intersection of half-planes, with the sides of the polygon being on the border of the half-planes.
If we can make sure that for each half-plane, polygon i is in that half-plane, then polygon i will be inside polygon i+1.
Note that we only have to check the vertices of polygon i to ensure that the entire polygon is within the half-plane.
Let the first polygon be P and the second polygon be Q. To check that P_k is within the half-plane from side Q_j Q_{j+1}, we can use cross products. See here for more details.
Let \overrightarrow{v}=\overrightarrow{Q_j P_k} and \overrightarrow{w}=\overrightarrow{Q_j Q_{j+1}}. We want to make sure that the z component of \overrightarrow{v} \times \overrightarrow{w} \le 0. We can substitute variables into the inequality to obtain (v_x+r_i-r_{i+1}-l_i+l_{i+1})w_y-(v_y+u_i-u_{i+1}-d_i+d_{i+1})w_x \le 0. After expanding, we get w_yr_i-w_yr_{i+1}-w_yl_i+w_yl_{i+1}-w_xu_i+w_xu_{i+1}+w_xd_i-w_xd_{i+1} \le v_yw_x-v_xw_y.
To satisfy the conditions for this problem, we should create these inequalities for all possible i, j, and k.
Lastly, we should minimize \sum_{i=1}^N r_i+l_i+u_i+d_i and enforce the variables to be non-negative.
We have our linear inequalities and a linear function to minimize (which is a linear program), so we can finish our solution with a linear programming solver.
SOLUTIONS:
Setter's Solution
#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define sd(x) scanf("%d", &(x))
#define pii pair<int, int>
#define F first
#define S second
#define all(c) ((c).begin()), ((c).end())
#define sz(x) ((int)(x).size())
#define ld long double
template<class T,class U>
ostream& operator<<(ostream& os,const pair<T,U>& p){
os<<"("<<p.first<<", "<<p.second<<")";
return os;
}
template<class T>
ostream& operator <<(ostream& os,const vector<T>& v){
os<<"{";
for(int i = 0;i < (int)v.size(); i++){
if(i)os<<", ";
os<<v[i];
}
os<<"}";
return os;
}
#ifdef LOCAL
#define cerr cout
#else
#endif
#define TRACE
#ifdef TRACE
#define trace(...) __f(#__VA_ARGS__, __VA_ARGS__)
template <typename Arg1>
void __f(const char* name, Arg1&& arg1){
cerr << name << " : " << arg1 << std::endl;
}
template <typename Arg1, typename... Args>
void __f(const char* names, Arg1&& arg1, Args&&... args){
const char* comma = strchr(names + 1, ',');cerr.write(names, comma - names) << " : " << arg1<<" | ";__f(comma+1, args...);
}
#else
#define trace(...)
#endif
typedef double T;
typedef vector<T> VT;
typedef vector<VT> VVT;
typedef vector<int> VI;
#define EPS 1e-9
struct LPSolver {
int m, n;
VI B, N;
VVT D;
LPSolver(const VVT &A, const VT &b, const VT &c) : m(b.size()), n(c.size()), N(n + 1), B(m), D(m + 2, VT(n + 2)) {
for (int i = 0; i < m; i++) for (int j = 0; j < n; j++) D[i][j] = A[i][j];
for (int i = 0; i < m; i++) {
B[i] = n + i;
D[i][n] = -1;
D[i][n + 1] = b[i];
}
for (int j = 0; j < n; j++) {
N[j] = j;
D[m][j] = -c[j];
}
N[n] = -1;
D[m + 1][n] = 1;
}
void Pivot(int r, int s) {
for (int i = 0; i < m + 2; i++) if (i != r) {
for (int j = 0; j < n + 2; j++) if (j != s) {
D[i][j] -= D[r][j] * D[i][s] / D[r][s];
}
}
for (int j = 0; j < n + 2; j++) if (j != s) D[r][j] /= D[r][s];
for (int i = 0; i < m + 2; i++) if (i != r) D[i][s] /= -D[r][s];
D[r][s] = 1.0 / D[r][s];
swap(B[r], N[s]);
}
int Simplex(int phase) {
int x = phase == 1 ? m + 1 : m;
while (1) {
int s = -1;
for (int j = 0; j <= n; j++) {
if (phase == 2 && N[j] == -1) continue;
if (s == -1 || D[x][j] < D[x][s] || D[x][j] == D[x][s] && N[j] < N[s]) s = j;
}
if (D[x][s] > -EPS) return 1;
int r = -1;
for (int i = 0; i < m; i++) {
if (D[i][s] < EPS) continue;
if (r == -1 || D[i][n + 1] / D[i][s] < D[r][n + 1] / D[r][s] || D[i][n + 1] / D[i][s] == D[r][n + 1] / D[r][s] && B[i] < B[r]) r = i;
}
if (r == -1) return 0;
Pivot(r, s);
}
}
T Solve(VT& x) {
int r = 0;
for (int i = 1; i < m; i++) if (D[i][n + 1] < D[r][n + 1]) r = i;
if (D[r][n + 1] < -EPS) {
Pivot(r, n);
if (!Simplex(1) || D[m + 1][n + 1] < -EPS) return -numeric_limits<T>::infinity();
for (int i = 0; i < m; i++) if (B[i] == -1) {
int s = -1;
for (int j = 0; j <= n; j++) {
if (s == -1 || D[i][j] < D[i][s] || D[i][j] == D[i][s] && N[j] < N[s]) s = j;
}
Pivot(i, s);
}
}
if (!Simplex(2)) return numeric_limits<T>::infinity();
x = VT(n);
for (int i = 0; i < m; i++) if (B[i] < n) x[B[i]] = D[i][n + 1];
return D[m][n + 1];
}
};
ll area(const vector<pii> & vec){
int k = sz(vec);
ll ar = 0;
for(int i = 0; i < k; i++){
pii P = vec[i], Q = vec[(i + 1) % k];
ar += P.F * (ll) Q.S - P.S * (ll) Q.F;
}
return abs(ar);
}
int main(){
int n; sd(n);
vector<vector<pii>> vec2(n), vec(n);
vector<ll> ar(n);
for(int i = 0; i < n; i++){
int k; sd(k);
vec2[i].resize(k);
for(int j = 0;j < k; j++) sd(vec2[i][j].F), sd(vec2[i][j].S);
ar[i] = area(vec2[i]);
}
VVT A;
VT b;
VT c(6 * n);
fill(c.begin() + 4 * n, c.end(), -1);
for(int i = 0; i < 2 * n; i++){
VT F(6 * n);
F[i] = -1;
F[i + 2 * n] = 1;
F[i + 4 * n] = -1;
A.push_back(F);
b.push_back(0);
F[i] = 1; F[i + 2 * n] = -1;
A.push_back(F);
b.push_back(0);
}
vector<int> perm(n);
iota(all(perm), 0);
sort(all(perm), [&](int i, int j){
return ar[i] < ar[j];
});
for(int i = 1; i < n; i++)if(!(area(vec2[perm[i]]) >= area(vec2[perm[i - 1]]))){
assert(0);
}
for(int i = 0; i < n; i++)
vec[i] = vec2[perm[i]];
for(int i = 1; i < n; i++){
int k = sz(vec[i]);
for(int j = 0; j < k; j++){
int u = vec[i][j].F, v = vec[i][j].S;
int p = vec[i][(j+1) % k].F - vec[i][j].F, q = vec[i][(j + 1) % k].S - vec[i][j].S;
int r = i - 1;
for(pii P : vec[r]){
VT lft(6 * n);
int x = P.F, y = P.S;
lft[i] = q; lft[i + 2 * n] = -q;
lft[r] = -q; lft[r + 2 * n] = q;
lft[i + n] = -p; lft[i + 3 * n] = p;
lft[r + n] = p; lft[r + 3 * n] = -p;
double val = (x - u) * (ll) q - (y - v) * (ll) p;
for(auto & it : lft) it *= -1; val *= -1;
A.push_back(lft);
b.push_back(val);
}
}
}
LPSolver lps(A, b, c);
VT y;
double z = lps.Solve(y);
if(isnan(z) || -z > 1e7){
printf("-1\n");
} else printf("%.10lf\n", -z);
}
Tester's Solution
#include <bits/stdc++.h>
#define endl '\n'
#define SZ(x) ((int)x.size())
#define ALL(V) V.begin(), V.end()
#define L_B lower_bound
#define U_B upper_bound
#define pb push_back
using namespace std;
template<class T, class T1> int chkmin(T &x, const T1 &y) { return x > y ? x = y, 1 : 0; }
template<class T, class T1> int chkmax(T &x, const T1 &y) { return x < y ? x = y, 1 : 0; }
const int MAXN = (1 << 8);
// https://github.com/alxli/Algorithm-Anthology
template<class Matrix>
int simplex_solve(const Matrix &a, const std::vector<double> &b,
const std::vector<double> &c, std::vector<double> *x,
const bool MAXIMIZE = true, const double EPS = 1e-10) {
int m = a.size(), n = c.size();
Matrix t(m + 2, std::vector<double>(n + 2));
t[1][1] = 0;
for (int j = 1; j <= n; j++) {
t[1][j + 1] = MAXIMIZE ? c[j - 1] : -c[j - 1];
}
for (int i = 1; i <= m; i++) {
for (int j = 1; j <= n; j++) {
t[i + 1][j + 1] = -a[i - 1][j - 1];
}
t[i + 1][1] = b[i - 1];
}
for (int j = 1; j <= n; j++) {
t[0][j + 1] = j;
}
for (int i = n + 1; i <= m + n; i++) {
t[i - n + 1][0] = i;
}
double p1 = 0, p2 = 0;
bool done = true;
do {
double mn = std::numeric_limits<double>::max(), xmax = 0, v;
for (int j = 2; j <= n + 1; j++) {
if (t[1][j] > 0 && t[1][j] > xmax) {
p2 = j;
xmax = t[1][j];
}
}
for (int i = 2; i <= m + 1; i++) {
v = fabs(t[i][1] / t[i][p2]);
if (t[i][p2] < 0 && mn > v) {
mn = v;
p1 = i;
}
}
std::swap(t[p1][0], t[0][p2]);
for (int i = 1; i <= m + 1; i++) {
if (i != p1) {
for (int j = 1; j <= n + 1; j++) {
if (j != p2) {
t[i][j] -= t[p1][j]*t[i][p2] / t[p1][p2];
}
}
}
}
t[p1][p2] = 1.0 / t[p1][p2];
for (int j = 1; j <= n + 1; j++) {
if (j != p2) {
t[p1][j] *= fabs(t[p1][p2]);
}
}
for (int i = 1; i <= m + 1; i++) {
if (i != p1) {
t[i][p2] *= t[p1][p2];
}
}
for (int i = 2; i <= m + 1; i++) {
if (t[i][1] < 0) {
return -1;
}
}
done = true;
for (int j = 2; j <= n + 1; j++) {
if (t[1][j] > 0) {
done = false;
}
}
} while (!done);
x->clear();
for (int j = 1; j <= n; j++) {
for (int i = 2; i <= m + 1; i++) {
if (fabs(t[i][0] - j) < EPS) {
x->push_back(t[i][1]);
}
}
}
return 0;
}
typedef double T;
typedef vector<T> VT;
typedef vector<VT> VVT;
typedef vector<int> VI;
#define EPS 1e-10
struct LPSolver {
int m, n;
VI B, N;
VVT D;
LPSolver(const VVT &A, const VT &b, const VT &c) : m(b.size()), n(c.size()), N(n + 1), B(m), D(m + 2, VT(n + 2)) {
for (int i = 0; i < m; i++) for (int j = 0; j < n; j++) D[i][j] = A[i][j];
for (int i = 0; i < m; i++) {
B[i] = n + i;
D[i][n] = -1;
D[i][n + 1] = b[i];
}
for (int j = 0; j < n; j++) {
N[j] = j;
D[m][j] = -c[j];
}
N[n] = -1;
D[m + 1][n] = 1;
}
void Pivot(int r, int s) {
for (int i = 0; i < m + 2; i++) if (i != r) {
for (int j = 0; j < n + 2; j++) if (j != s) {
D[i][j] -= D[r][j] * D[i][s] / D[r][s];
}
}
for (int j = 0; j < n + 2; j++) if (j != s) D[r][j] /= D[r][s];
for (int i = 0; i < m + 2; i++) if (i != r) D[i][s] /= -D[r][s];
D[r][s] = 1.0 / D[r][s];
swap(B[r], N[s]);
}
int Simplex(int phase) {
int x = phase == 1 ? m + 1 : m;
while (1) {
int s = -1;
for (int j = 0; j <= n; j++) {
if (phase == 2 && N[j] == -1) continue;
if (s == -1 || D[x][j] < D[x][s] || D[x][j] == D[x][s] && N[j] < N[s]) s = j;
}
if (D[x][s] > -EPS) return 1;
int r = -1;
for (int i = 0; i < m; i++) {
if (D[i][s] < EPS) continue;
if (r == -1 || D[i][n + 1] / D[i][s] < D[r][n + 1] / D[r][s] || D[i][n + 1] / D[i][s] == D[r][n + 1] / D[r][s] && B[i] < B[r]) r = i;
}
if (r == -1) return 0;
Pivot(r, s);
}
}
T Solve(VT& x) {
int r = 0;
for (int i = 1; i < m; i++) if (D[i][n + 1] < D[r][n + 1]) r = i;
if (D[r][n + 1] < -EPS) {
Pivot(r, n);
if (!Simplex(1) || D[m + 1][n + 1] < -EPS) return -numeric_limits<T>::infinity();
for (int i = 0; i < m; i++) if (B[i] == -1) {
int s = -1;
for (int j = 0; j <= n; j++) {
if (s == -1 || D[i][j] < D[i][s] || D[i][j] == D[i][s] && N[j] < N[s]) s = j;
}
Pivot(i, s);
}
}
if (!Simplex(2)) return numeric_limits<T>::infinity();
x = VT(n);
for (int i = 0; i < m; i++) if (B[i] < n) x[B[i]] = D[i][n + 1];
return D[m][n + 1];
}
};
struct PT {
int x, y;
PT() { x = y = 0; }
PT(int _x, int _y) { x = _x; y = _y; }
PT operator-(const PT &other) {
return PT(x - other.x, y - other.y);
}
PT operator+(const PT &other) {
return PT(x + other.x, y + other.y);
}
int operator&(const PT &other) {
return x * other.y - y * other.x;
}
};
int n;
vector<PT> vec[MAXN];
int area[MAXN];
void read() {
cin >> n;
for(int i = 0; i < n; i++) {
int k;
cin >> k;
while(k--) {
int x, y;
cin >> x >> y;
vec[i].pb(PT(x, y));
}
int area = 0;
for(int j = 0, z = 1; j < SZ(vec[i]); j++, z = (z + 1) % SZ(vec[i])) {
area += vec[i][j] & vec[i][z];
}
area = max(area, -area);
::area[i] = area;
}
}
vector<int> X,Y;
vector<vector<double> > A;
vector<double> b,c;
double z;
int m;
void pivot(int x,int y){
swap(X[y],Y[x]);
b[x]/=A[x][y];
for(int i = 0; i < m; i++) if(i!=y)A[x][i]/=A[x][y];
A[x][y]=1/A[x][y];
for(int i = 0; i < n; i++) if(i!=x&&abs(A[i][y])>EPS){
b[i]-=A[i][y]*b[x];
for(int j = 0; j < m; j++) if(j!=y)A[i][j]-=A[i][y]*A[x][j];
A[i][y]=-A[i][y]*A[x][y];
}
z+=c[y]*b[x];
for(int i = 0; i < m; i++) if(i!=y)c[i]-=c[y]*A[x][i];
c[y]=-c[y]*A[x][y];
}
pair<double,vector<double> > simplex( // maximize c^T x s.t. Ax<=b, x>=0
vector<vector<double> > _A, vector<double> _b, vector<double> _c){
// returns pair (maximum value, solution vector)
A=_A;b=_b;c=_c;
n=b.size();m=c.size();z=0.;
X=vector<int>(m);Y=vector<int>(n);
for(int i = 0; i < m; i++)X[i]=i;
for(int i = 0; i < n; i++)Y[i]=i+m;
while(1){
int x=-1,y=-1;
double mn=-EPS;
for(int i = 0; i < n; i++)if(b[i]<mn)mn=b[i],x=i;
if(x<0)break;
for(int i = 0; i < m; i++)if(A[x][i]<-EPS){y=i;break;}
if(y < 0) return {-1, vector<double>()};
pivot(x,y);
}
while(1){
double mx=EPS;
int x=-1,y=-1;
for(int i = 0; i < m; i++)if(c[i]>mx)mx=c[i],y=i;
if(y<0)break;
double mn=1e200;
for(int i = 0; i < n; i++)if(A[i][y]>EPS&&b[i]/A[i][y]<mn)mn=b[i]/A[i][y],x=i;
if(x < 0) return {-1, vector<double>()};
pivot(x,y);
}
vector<double> r(m);
for(int i = 0; i < n; i++)if(Y[i]<m)r[Y[i]]=b[i];
return make_pair(z,r);
}
int perm[MAXN];
bool cmp(int i, int j) { return area[i] < area[j]; }
void solve() {
for(int i = 0; i < n; i++) {
perm[i] = i;
}
sort(perm, perm + n, cmp);
// Standard idea with 4*n variables for translation on x and y
int cnt_vars = 4 * n;
vector<vector<double> > A;
vector<double> b, c(cnt_vars, -1), emp(cnt_vars, 0), x;
for(int i = 0; i + 1 < n; i++) {
for(int j = 0, z = 1; j < SZ(vec[perm[i + 1]]); j++, z = (z + 1) % SZ(vec[perm[i + 1]])) {
PT ab = vec[perm[i + 1]][z] - vec[perm[i + 1]][j], o = vec[perm[i + 1]][j];
for(PT &p: vec[perm[i]]) {
// ab.x * (p.y - o.y + d1y[i] - d2y[i] - d1y[i + 1] + d2y[i + 1])
// - ab.y * (p.x - o.x + d1x[i] - d2x[i] - d1x[i + 1] + d2x[i + 1]) >= 0
vector<double> row = emp;
row[perm[i]] = ab.y;
row[perm[i] + n] = -ab.y;
row[perm[i + 1]] = -ab.y;
row[perm[i + 1] + n] = ab.y;
row[perm[i] + 2 * n] = -ab.x;
row[perm[i] + 3 * n] = ab.x;
row[perm[i + 1] + 2 * n] = ab.x;
row[perm[i + 1] + 3 * n] = -ab.x;
A.pb(row);
b.pb(-ab.y * (p.x - o.x) + ab.x * (p.y - o.y));
}
}
}
/*
int ret = simplex_solve(A, b, c, &x, false);
if(ret == -1) {
cout << -1 << endl;
return;
}
double ans = 0;
for(int i = 0; i < cnt_vars; i++) {
ans += c[i] * x[i];
}
*/
LPSolver lps(A, b, c);
VT y;
double ans = lps.Solve(y);
if(isnan(ans) || -ans > 1e7){
cout << -1 << endl;
return;
}
//ans = simplex(A, b, c).first;
cout << setprecision(9) << fixed << -ans << endl;
}
int main() {
ios_base::sync_with_stdio(false);
cin.tie(NULL);
read();
solve();
return 0;
}
Editorialist's Solution
#include <bits/stdc++.h>
using namespace std;
//https://github.com/ivangalban/TeamReference/blob/master/Numeric/simplex.cpp
const double eps = 1e-9, oo = numeric_limits<double>::infinity();
typedef vector<double> vec;
typedef vector<vec> mat;
double simplexMethodPD(mat &A, vec &b, vec &c)
{
int n = c.size(), m = b.size();
mat T(m + 1, vec(n + m + 1));
vector<int> base(n + m), row(m);
for(int j = 0; j < m; ++j)
{
for (int i = 0; i < n; ++i)
T[j][i] = A[j][i];
T[j][n + j] = 1;
base[row[j] = n + j] = 1;
T[j][n + m] = b[j];
}
for (int i = 0; i < n; ++i)
T[m][i] = c[i];
while (1)
{
int p = 0, q = 0;
for (int i = 0; i < n + m; ++i)
if (T[m][i] <= T[m][p])
p = i;
for (int j = 0; j < m; ++j)
if (T[j][n + m] <= T[q][n + m])
q = j;
double t = min(T[m][p], T[q][n + m]);
if (t >= -eps)
{
vec x(n);
for (int i = 0; i < m; ++i)
if (row[i] < n) x[row[i]] = T[i][n + m];
// x is the solution
return -T[m][n + m]; // optimal
}
if (t < T[q][n + m])
{
// tight on c -> primal update
for (int j = 0; j < m; ++j)
if (T[j][p] >= eps)
if (T[j][p] * (T[q][n + m] - t) >=
T[q][p] * (T[j][n + m] - t))
q = j;
if (T[q][p] <= eps)
return oo; // primal infeasible
}
else
{
// tight on b -> dual update
for (int i = 0; i < n + m + 1; ++i)
T[q][i] = -T[q][i];
for (int i = 0; i < n + m; ++i)
if (T[q][i] >= eps)
if (T[q][i] * (T[m][p] - t) >=
T[q][p] * (T[m][i] - t))
p = i;
if (T[q][p] <= eps)
return -oo; // dual infeasible
}
for (int i = 0; i < m + n + 1; ++i)
if (i != p) T[q][i] /= T[q][p];
T[q][p] = 1; // pivot(q, p)
base[p] = 1;
base[row[q]] = 0;
row[q] = p;
for (int j = 0; j < m + 1; ++j)
if (j != q)
{
double alpha = T[j][p];
for (int i = 0; i < n + m + 1; ++i)
T[j][i] -= T[q][i] * alpha;
}
}
return oo;
}
//actual solution
struct polygon {
vector<array<int, 2>> a;
int area=0;
polygon() {
//input polygon
int m;
cin >> m;
for(int i=0, x, y; i<m; ++i) {
cin >> x >> y;
a.push_back({x, y});
}
//find 2*area
for(int i=0; i<m; ++i)
area+=a[i][0]*a[(i+1)%m][1]-a[(i+1)%m][0]*a[i][1];
}
bool operator<(const polygon &o) const {
return area<o.area;
}
};
int n;
vector<polygon> p;
int main() {
ios::sync_with_stdio(0);
cin.tie(0);
//input
cin >> n;
for(int i=0; i<n; ++i)
p.push_back(polygon());
//sort polygons by size
sort(p.begin(), p.end());
//find the constraints
mat A;
vec b;
for(int i=0; i+1<n; ++i) {
//constraints between p[i] and p[i+1]
int s=p[i+1].a.size();
for(int j=0; j<s; ++j) {
//process half-plane of side j of p[i+1]
for(int k=0; k<p[i].a.size(); ++k) {
//make sure point k is in the half-plane
int ux=p[i+1].a[(j+1)%s][0]-p[i+1].a[j][0], uy=p[i+1].a[(j+1)%s][1]-p[i+1].a[j][1];
int vx=p[i].a[k][0]-p[i+1].a[j][0], vy=p[i].a[k][1]-p[i+1].a[j][1];
//add constraint for the point and half-plane
vec a(4*n);
a[4*i]=uy;
a[4*(i+1)]=-uy;
a[4*i+1]=-uy;
a[4*(i+1)+1]=uy;
a[4*i+2]=-ux;
a[4*(i+1)+2]=ux;
a[4*i+3]=ux;
a[4*(i+1)+3]=-ux;
A.push_back(a);
b.push_back(vy*ux-vx*uy);
}
}
}
//cost vector
vec c(4*n, 1);
//solve!
double ans=simplexMethodPD(A, b, c);
if(isinf(ans))
cout << -1;
else
cout << fixed << setprecision(9) << ans;
}
Please give me suggestions if anything is unclear so that I can improve. Thanks 
