poj1259 The picnic【最大空凸包详解 计算几何+动态规划】

题目大意：

给定平面上$n$个点，求面积最大的空凸包。$1\le n\le 100$

解题思路：

感谢http://blog.youkuaiyun.com/nyroro/article/details/45268767这篇博客的讲解。

算法一

先枚举每个点作为凸包的最下面的点O，不考虑它下方的点，其余点以O为原点极角排序。

再枚举凸包最后一个点$i$与O相连,设$dp[i][j](i>j)$表示组成凸包的最后一个三角形是$Oij$的最大凸包面积，那么容易得到dp方程：

$dp[i][j]=max(S_\Delta Oij+dp[j][k]),(\Delta Oij内无点且边Oi上无点且k在\vec {ij}右侧保证凸性)$；

注意到如果边$Oi$上有点，它就只能作为凸包的边界，所以不能更新$dp[i][j]$，但还是要用$max(S_\Delta Oij+dp[j][k])$更新最终答案，否则用$dp[i][j]$更新答案。

但这样单次dp的复杂度是$O(n^3)$的，考虑优化这个过程，所以有了算法二。

算法二

还是枚举凸包最后一个点$i$，dp意义也相同。可以发现，合法的$j$可用以下方法得到：

1.令$j_1=i-1$,可以保证$\Delta Oij$合法；如果$j_1$在$Oi$上，那么从大到小枚举第一个不在$Oi$上的点作为$j_1$，但这样$Oi$上就有点了，处理方法同上。

2.令$j_2$为最大的在$\vec{ij_1}$右侧的点，$j_3$为最大的在$\vec{ij_2}$右侧的点……可以证明这样可以找出所有合法的$j$。

3.注意到找出的$j$序列是满足凸性的，且$dp[i][j_n]$的第一个可以转移来的点就是$j_{n+1}$,所以可以用$g[i][j]$表示$max(dp[i][k],1\le k\le j)$，那么就有$dp[i][j_n]=S_\Delta Oij+g[i][j_{n+1}]$
注：1~j中不合法的点dp值为0，不影响$g$的更新，所以$g$数组一定满足凸性；但如果$Oi$上有点还是不更新$dp$和$g$，但还是要更新最终答案。

这样枚举原点是$O(n)$的，枚举凸包最后一点$i$是$O(n)$的，枚举合法的$j$是$O(n)$的，转移是$O(1)$的，更新$g$数组是$O(n)$的，所以总复杂度是$O(n^3)$的。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<cmath>
using namespace std;

int getint()
{
    int i=0,f=1;char c;
    for(c=getchar();(c!='-')&&(c<'0'||c>'9');c=getchar());
    if(c=='-')f=-1,c=getchar();
    for(;c>='0'&&c<='9';c=getchar())i=(i<<3)+(i<<1)+c-'0';
    return i*f;
}

const int N=105;
struct point
{
    double x,y;
    point(){}
    point(double _x,double _y):x(_x),y(_y){}
    inline friend point operator - (const point &a,const point &b)
    {return point(a.x-b.x,a.y-b.y);}
    inline friend double operator * (const point &a,const point &b)
    {return a.x*b.y-a.y*b.x;}
    inline double dis(){return x*x+y*y;}
}a[N],p[N],O;
int T,n,m;
double dp[N][N],ans;

inline bool cmp(const point &a,const point &b)
{
    double res=(a-O)*(b-O);
    if(res)return res>0;
    return (a-O).dis()<(b-O).dis();
}

void solve()
{
    memset(dp,0,sizeof(dp));
    sort(p+1,p+m+1,cmp);
    for(int i=1;i<=m;i++)
    {
        int j=i-1;
        while(j&&!((p[i]-O)*(p[j]-O)))j--;
        bool bz=(j==i-1);
        while(j)
        {
            int k=j-1;
            while(k&&(p[i]-p[k])*(p[j]-p[k])>0)k--;
            double area=fabs((p[i]-O)*(p[j]-O))/2;
            if(k)area+=dp[j][k];
            if(bz)dp[i][j]=area;
            ans=max(ans,area),j=k;
        }
        if(bz)for(int j=1;j<i;j++)dp[i][j]=max(dp[i][j],dp[i][j-1]);
    }
}

int main()
{
    //freopen("lx.in","r",stdin);
    T=getint();
    while(T--)
    {
        n=getint();ans=0;
        for(int i=1;i<=n;i++)a[i].x=getint(),a[i].y=getint();
        for(int i=1;i<=n;i++)
        {
            O=a[i],m=0;
            for(int j=1;j<=n;j++)
                if(a[j].y>a[i].y||a[j].y==a[i].y&&a[j].x>a[i].x)p[++m]=a[j];
            solve();
        }
        printf("%0.1lf\n",ans);
    }
    return 0;
}