首先，肯定只有凸包上的点会限制这个矩形，所以建立凸包。
然后可以知道，矩形上一定有一条边与凸包上的边重合，否则可以转一下使得它重合，答案会更小。
于是沿着凸包枚举这一条边，通过旋转卡壳找到离这条边最远的点以及这个矩形两端的点，这五个点构成的矩形就是一个可能的答案了。
各种判断用向量叉积和点积
注意一下输出\(-0.0000\)的情况
#include<bits/stdc++.h>
#define ld long double
#define eps 1e-8
//This code is written by Itst
using namespace std;
const int MAXN = 5e4 + 10;
struct vec{
ld x , y;
vec(ld _x = 0 , ld _y = 0){x = _x; y = _y;}
bool operator <(const vec a)const{
return x < a.x;
}
vec operator -(vec a){
return vec(x - a.x , y - a.y);
}
vec operator *(ld p){
return vec(p * x , p * y);
}
vec operator +(vec a){
return vec(x + a.x , y + a.y);
}
}now[MAXN] , temp[MAXN] , squ[4];
int cnt , N , top , st[MAXN] , ind[MAXN];
ld ans;
inline bool cmp(ld a , ld b){
return a - eps < b && a + eps > b;
}
inline ld cot(vec a , vec b){
return a.x * b.y - a.y * b.x;
}
inline ld dot(vec a , vec b){
return a.x * b.x + a.y * b.y;
}
inline ld calS(vec a , vec b , vec c){
return fabs(cot(c - a , b - a));
}
inline ld len(vec a){
return sqrt(a.x * a.x + a.y * a.y);
}
inline ld calS(vec a , vec b , vec c , vec d , vec e){
return calS(a , d , b) / len(b - a) * dot(c - e , b - a) / len(b - a);
}
inline vec rev(vec a){
return vec(-a.y , a.x);
}
void input(){
cin >> N;
for(int i = 1 ; i <= N ; ++i)
cin >> now[i].x >> now[i].y;
}
void init(){
ans = 1e18;
sort(now + 1 , now + N + 1);
for(int i = 1 ; i <= N ; ++i){
while(top >= 2 && cot(now[i] - now[st[top - 1]] , now[st[top]] - now[st[top - 1]]) > -eps)
--top;
st[++top] = i;
}
for(int i = 1 ; i <= top ; ++i)
ind[++cnt] = st[i];
top = 0;
for(int i = N ; i ; --i){
while(top >= 2 && cot(now[i] - now[st[top - 1]] , now[st[top]] - now[st[top - 1]]) > -eps)
--top;
st[++top] = i;
}
for(int i = 2 ; i < top ; ++i)
ind[++cnt] = st[i];
for(int i = 1 ; i <= cnt ; ++i)
temp[i] = now[ind[i]];
memcpy(now + 1 , temp + 1 , sizeof(vec) * cnt);
}
void work(){
int minX = 1 , maxX = 1 , minY = 1 , maxY = 1;
for(int i = 2 ; i <= N ; ++i){
if(now[minX].x > now[i].x)
minX = i;
if(now[maxX].x < now[i].x)
maxX = i;
if(now[minY].y > now[i].y)
minY = i;
if(now[maxY].y < now[i].y)
maxY = i;
}
ans = (now[maxY].y - now[minY].y) * (now[maxX].x - now[minX].x);
squ[0].x = squ[3].x = now[minX].x;
squ[1].x = squ[2].x = now[maxX].x;
squ[0].y = squ[1].y = now[minY].y;
squ[2].y = squ[3].y = now[maxY].y;
for(int i = 1 ; i <= cnt ; minY = minY % cnt + 1 , ++i){
while(calS(now[minY] , now[maxY] , now[minY % cnt + 1]) < calS(now[minY] , now[maxY % cnt + 1] , now[minY % cnt + 1]))
maxY = maxY % cnt + 1;
while(dot(now[minY % cnt + 1] - now[minY] , now[minX % cnt + 1] - now[minX]) < eps)
minX = minX % cnt + 1;
while(dot(now[minY % cnt + 1] - now[minY] , now[maxX % cnt + 1] - now[maxX]) > -eps)
maxX = maxX % cnt + 1;
ld t = calS(now[minY] , now[minY % cnt + 1] , now[maxX] , now[maxY] , now[minX]);
if(t < ans){
ans = t;
squ[0] = (now[minY % cnt + 1] - now[minY]) * (dot(now[minY] - now[maxX] , now[minY] - now[minY % cnt + 1]) / len(now[minY] - now[minY % cnt + 1]) / len(now[minY] - now[minY % cnt + 1])) + now[minY];
squ[1] = rev(now[minY % cnt + 1] - now[minY]) * (calS(now[minY] , now[minY % cnt + 1] , now[maxY]) / len(now[minY] - now[minY % cnt + 1]) / len(now[minY] - now[minY % cnt + 1])) + squ[0];
squ[2] = (now[minY] - now[minY % cnt + 1]) * (t / len(squ[1] - squ[0]) / len(now[minY % cnt + 1] - now[minY])) + squ[1];
squ[3] = squ[2] + (squ[0] - squ[1]);
}
}
}
void output(){
cout << fixed << setprecision(5) << (ans < 1e-5 ? 0 : ans) << endl;
int dir = 0;
for(int j = 1 ; j < 4 ; ++j)
if(squ[dir].y > squ[j].y || cmp(squ[dir].y , squ[j].y) && squ[j].x < squ[dir].x)
dir = j;
for(int i = 0 ; i < 4 ; ++i)
cout << fixed << setprecision(5) << (squ[(dir + i) % 4].x < 1e-5 ? 0 : squ[(dir + i) % 4].x) << ' ' << (squ[(dir + i) % 4].y < 1e-5 ? 0 : squ[(dir + i) % 4].y) << endl;
}
int main(){
#ifndef ONLINE_JUDGE
freopen("in" , "r" , stdin);
//freopen("out" , "w" , stdout);
#endif
input();
init();
work();
output();
return 0;
}