Smallest Enclosing Square
National Olympiad — National round (2013 season)
Submission name: land.cpp
Problem summary
Find the side length of the smallest square, with any orientation, containing a given convex quadrilateral.
Input
Four pairs of coordinates in cyclic boundary order.
Output
The minimum side length formatted to three decimals.
Constraints
Integer coordinates in [0,1000]; the quadrilateral is strictly convex and no three vertices are collinear. Full credit requires absolute error at most 0.001.
Worked example
An independently prepared example for this explanation.
Input
10 6 5 0 3 3 3 10
Output
7.938
Solution analysis
- For an orientation θ, project all vertices onto u=(cosθ,sinθ) and v=(−sinθ,cosθ). The enclosing square needs side max(width(θ),height(θ)). Orientations repeat every π/2.
- Projection extrema change only when an axis becomes perpendicular to the difference of two vertices. Collect those angles modulo π/2, together with 0 and π/2, and split the interval there. Within each open interval, the four supporting vertices are fixed, so width and height are each A cosθ+B sinθ.
- Each positive width or height function is concave on its support interval because its second derivative is its negative. Thus the minimum of their maximum occurs at an interval endpoint or where width equals height. Solve the latter trigonometric equation analytically and evaluate all candidates. For the sample this gives side 7.938; restricting to axis-aligned squares would miss the optimum.
Why it works
Projection ranges are necessary and sufficient for containment in an oriented square. The breakpoints include every change of supporting vertices. Between them, away from a width-height intersection, the larger range is one fixed positive concave function and has its minimum at an endpoint. Therefore endpoints and intersections contain a global minimizer, all of which the algorithm evaluates.
Watch out. A square may be rotated. Minimize the larger projection range, not the bounding-box area. Include angle boundaries and use long double for trigonometric evaluation.
Commented C++17 solution
An independently written complete program with standard headers.
#include<iostream>
#include<vector>
#include<array>
#include<algorithm>
#include<cmath>
#include<iomanip>
// Algorithm: For an orientation θ, project all vertices onto u=(cosθ,sinθ) and v=(−sinθ,cosθ).
using namespace std;
struct Point{
long double x;
long double y;
};
int main(){
array<Point,4>points;
for(Point&point:points){
cin>>point.x>>point.y;
}
const long double pi=acosl(-1.0L);
const long double quarter=pi/2;
vector<long double>angles={0,quarter};
for(int i=0;i<4;i++){
for(int j=i+1;j<4;j++){
long double angle=fmodl(atan2l(points[j].y-points[i].y,points[j].x-points[i].x),quarter);
if(angle<0){
angle+=quarter;
}
angles.push_back(angle);
}
}
sort(angles.begin(),angles.end());
auto ranges=[&](long double angle){
long double cosine=cosl(angle);
long double sine=sinl(angle);
long double minX=1e30L;
long double maxX=-1e30L;
long double minY=1e30L;
long double maxY=-1e30L;
for(Point point:points){
long double x=point.x*cosine+point.y*sine;
long double y=-point.x*sine+point.y*cosine;
minX=min(minX,x);
maxX=max(maxX,x);
minY=min(minY,y);
maxY=max(maxY,y);
}
return max(maxX-minX,maxY-minY);
};
long double answer=1e30L;
for(long double angle:angles){
answer=min(answer,ranges(angle));
}
for(size_t interval=1;interval<angles.size();interval++){
long double middle=(angles[interval-1]+angles[interval])/2;
long double cosine=cosl(middle);
long double sine=sinl(middle);
int minX=0;
int maxX=0;
int minY=0;
int maxY=0;
auto xProjection=[&](int i){
return points[i].x*cosine+points[i].y*sine;
};
auto yProjection=[&](int i){
return -points[i].x*sine+points[i].y*cosine;
};
for(int i=1;i<4;i++){
if(xProjection(i)<xProjection(minX)){
minX=i;
}
if(xProjection(i)>xProjection(maxX)){
maxX=i;
}
if(yProjection(i)<yProjection(minY)){
minY=i;
}
if(yProjection(i)>yProjection(maxY)){
maxY=i;
}
}
long double a=points[maxX].x-points[minX].x-points[maxY].y+points[minY].y;
long double b=points[maxX].y-points[minX].y+points[maxY].x-points[minY].x;
if(fabsl(a)+fabsl(b)<1e-18L){
continue;
}
long double crossing=atan2l(-a,b);
for(int turn=-1;turn<=1;turn++){
long double angle=crossing+turn*pi;
if(angle>=angles[interval-1]&&angle<=angles[interval]){
answer=min(answer,ranges(angle));
}
}
}
cout<<fixed<<setprecision(3)<<answer<<'\n';
}