POLCHAIN - Editorial

PROBLEM LINK:

Practice
Div-1 Contest

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 :slight_smile:

8 Likes

Can you provide some good resources on linear programming solver ?
How it works ? When can we use that ? What is time complexity of solving this problem ?

1 Like

Simplex Algorithm : http://fourier.eng.hmc.edu/e176/lectures/NM/node32.html

4 Likes

how to solve for subtask 3 where all polygons are axis aligned rectangles

1 Like

Use DP and check every possible placement of all rectangles.
To optimize it further we can just solve for optimal cost because of x-axis and cost because of y-axis individually ( because since rectangles are axis aligned, the cost due to x-axis will be independent of cost due to y-axis).

Let’s define a rectangle with 4 points
mnx,mxx,mny,mxy.
( mn =minimum , mx = maximum)

  1. sort the rectangles according to area and then see if rectangle with smaller area can lie inside bigger rectangle or not. Print -1 if it’s not possible.
  2. Now,Let’s solve for x coordinate, ( we will use same method to solve for y coordinate)
    Let’s formulate our DP
    dp[i][j] denotes minimum cost when ith rectangle’s mnx ( that is minimum x coordinate in given rectangle) lies exactly at place j and all other smaller rectangles lies inside it( note that j can be negative as well so we need to handle that by shifting graph).

dp[i][j]= min{dp[i-1][k]+abs(mnx-j)}
where k is such that i-1th rectangle will lie inside ith rectangle if we place mnx of ith rectangle at jth position.

Start dp with smallest rectangle. Find cost of placing it’s mnx at jth position and Mark cost in dp[0][j].
Then go for next rectangle,
for each j, check all possible places allowed for previous rectangle and store the minimum value at dp[i][j] ( includes cost of two i+1 rectangles)

Feel free to ask if you don’t understand anything.
Time complexity = O(n*k1*k1 + n*k2*k2)
Where,
k1=width of the grid and
k2=height of the grid
We can obviously use segtree and optimize it further to
O(n*k1*log(k1) + n*k2*log(k2) )
But it’s not required in given problem.
Link to implementation : https://www.codechef.com/viewsolution/29596030

3 Likes

Thanks a lot for sharing.

I think a lot of people treat linear programming solvers as a black box (like me). You only need to know when to apply it, and the link that @anand20 posted is good.

5 Likes

Hi,

I guess the solution is not taking care of regular polygons.Is the problem only dealing with convex polygon?,if that’s the case why its not mentioned in the problem. Problem says polygon which can be concave or any other polygons

The problem mentions about convex polygons…read carefully :slightly_smiling_face:

2 Likes