[BZOJ 1185][HNOI 2007]最小矩形覆盖(凸包+旋转卡壳)

题目链接

http://www.lydsy.com/JudgeOnline/problem.php?id=1185

思路

这个题真的很神奇啊。。。虽然有人说被卡精度了，不过反正我是没遇到这样的问题(三态函数大法好)
实际上这个题与我之前写的另一题（最大土地面积）的做法比较相近，也是先维护一个凸包，然后搞旋转卡壳。个人认为此题的旋转卡壳略难一些。
首先我们发现，一个最小覆盖矩形一定是与点集的凸包上的部分点相交，手模一下可以发现，这个矩形的一条边一定交了凸包上的相邻两个点，而其他三条边各至少交了凸包上的一个点，因此我们可以枚举那个交了相邻两个点的矩形的边（以下简称底边），然后再用旋转卡壳找出矩形的另外三条边上的那三个点。显然这三个点各是点集中的最左点、最右点、最上点(把枚举的底边看成水平并且是在凸包的最下面的话)，然后可以发现，卡壳过程中最上点与底边构成的三角形的面积是个单峰函数，而卡壳过程中最左点和最右点各自与底边上任一点的连线在底边上的投影长度也是个单峰函数。既然是单峰函数，那么就能拿旋转卡壳搞。得到底边以及上述的三个点后，求出这个矩形的面积和四个顶点也是比较容易的事情了，更新答案即可。注意最终答案的四个顶点的输出顺序！不过由于之前我在做凸包时，栈里凸包的所有点都是按照逆时针顺序排列的，因此最终输出答案时并不需要做太多处理。
最终此题的复杂度是$O(nlogn)$的。
下图是我对我的程序中的部分变量的解释

然后再贡献一个题目样例的图形表示，题目样例还是卡得很紧的，有凸包上连线垂直的情况，但是题目样例本身就是一个凸包，这一点需要注意

代码

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <cmath>

#define MAXN 51000
#define EPS 1e-10

using namespace std;

int n;
double minsqr; //minsqr=最小矩形面积

int dcmp(double x)
{
    if(fabs(x)<EPS) return 0;
    if(x>-EPS) return 1;
    return -1;
}

struct Point
{
    double x,y;
    Point(){}
    Point(double _x,double _y):x(_x),y(_y){}
    void read()
    {
        scanf("%lf%lf",&x,&y);
    }
}points[MAXN],stack[MAXN<<1],ans[10]; //ans保存最小覆盖矩形的四个顶点

Point operator-(Point a,Point b)
{
    return Point(a.x-b.x,a.y-b.y);
}

Point operator*(Point a,double b)
{
    return Point(a.x*b,a.y*b);
}

Point operator+(Point a,Point b)
{
    return Point(a.x+b.x,a.y+b.y);
}

double cross(Point a,Point b,Point c) //a->b 叉积 a->c
{
    return (b.x-a.x)*(c.y-a.y)-(c.x-a.x)*(b.y-a.y);
}

double dot(Point a,Point b,Point c) //a->b 点积 b->c
{
    return (b.x-a.x)*(c.x-b.x)+(b.y-a.y)*(c.y-b.y);
}

double dist(Point a,Point b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

int top=0;

bool cmp(Point a,Point b)
{
    if(!dcmp(a.x-b.x)) return a.y<b.y;
    return a.x<b.x;
}

bool operator<(Point a,Point b)
{
    if(!dcmp(a.y-b.y)) return a.x<b.x;
    return a.y<b.y;
}

void Graham()
{
    sort(points+1,points+n+1,cmp);
    for(int i=1;i<=n;i++) //下凸壳
    {
        while(top>=2&&dcmp(cross(stack[top-1],stack[top],points[i])<=0)) top--;
        stack[++top]=points[i];
    }
    int tmp=top;
    for(int i=n-1;i>=1;i--) //上凸壳
    {
        while(top>=tmp+1&&dcmp(cross(stack[top-1],stack[top],points[i])<=0)) top--;
        stack[++top]=points[i];
    }
    top--; //栈顶显然是points[1],重复了一次，将重复的丢掉
}

void rotcalip() //旋转卡壳
{
    int left=1,right=1,up=1; //左、右、上三个点
    for(int i=1;i<=top;i++)
    {
        while(dcmp(cross(stack[up+1],stack[i],stack[i+1])-cross(stack[up],stack[i],stack[i+1]))>=0) up=up%top+1;
        while(dcmp(dot(stack[right+1],stack[i+1],stack[i])-dot(stack[right],stack[i+1],stack[i]))>=0) right=right%top+1;
        if(i==1) left=right;
        while(dcmp(dot(stack[left+1],stack[i],stack[i+1])-dot(stack[left],stack[i],stack[i+1]))>=0) left=left%top+1;
        double L=dist(stack[i],stack[i+1]); //L=点i到i+1的距离
        double H=cross(stack[i+1],stack[up],stack[i])/L; //H是矩形的高
        double bottom=(fabs(dot(stack[i],stack[i+1],stack[left])/L)+fabs(dot(stack[i],stack[i+1],stack[right])/L));
        double nowsqr=bottom*H;
        if(nowsqr<minsqr) //这个解比当前答案更优
        {
            minsqr=nowsqr;
            ans[0]=stack[i]+(stack[i+1]-stack[i])*((fabs(dot(stack[i],stack[i+1],stack[right]))/L+L)/L);
            ans[1]=ans[0]+(stack[right]-ans[0])*(H/dist(stack[right],ans[0]));
            ans[2]=ans[1]+(stack[up]-ans[1])*(bottom/dist(stack[up],ans[1]));
            ans[3]=ans[2]+(stack[left]-ans[2])*(H/dist(stack[left],ans[2]));
        }
    }
}

int main()
{
    minsqr=1e12;
    scanf("%d\n",&n);
    for(int i=1;i<=n;i++)
        points[i].read();
    Graham();
    rotcalip();
    printf("%.5lf\n",minsqr);
    int tmp=0;
    for(int i=0;i<=3;i++)
    {
        if(ans[i]<ans[tmp]) tmp=i;
    }
    for(int i=0;i<=3;i++)
    {
        printf("%.5lf %.5lf\n",ans[(i+tmp)%4].x,ans[(i+tmp)%4].y);
    }
    return 0;
}