LiveArchive 6503 : Golf Field

問題文 : https://icpcarchive.ecs.baylor.edu/external/65/6503.pdf

問題概要 :
２次元平面上にn個の点がある. ( 4 <= n <= 30,000 )
この中から異なる４点を選び凸包を作る.
得られる凸包の面積の最大値を求めよ.

解法 :
凸包 + キャリパー法 + 三分探索

以下の問題の完全上位互換
ncastar.hatenablog.com

解法は基本的に同じだが、
上記の問題よりも頂点数が多くなっているため愚直に２重ループを回せない
そのため２重ループの代わりにキャリパー法を使う

コード :

#include<bits/stdc++.h>
 
#define REP(i,s,n) for(int i=s;i<n;++i)
#define rep(i,n) REP(i,0,n)
#define EPS (1e-10)
#define equals(a,b) (fabs((a)-(b)) < EPS)
#define COUNTER_CLOCKWISE 1
#define CLOCKWISE -1 
#define ONLINE_BACK 2
#define ONLINE_FRONT -2
#define ON_SEGMENT 0
 
using namespace std;
 
// Geo --->
 
class Point{
public:
  double x,y;
 
  Point(double x = 0,double y = 0): x(x),y(y){}
 
  Point operator + (Point p){return Point(x+p.x,y+p.y);}
  Point operator - (Point p){return Point(x-p.x,y-p.y);}
  Point operator * (double a){return Point(a*x,a*y);}
  Point operator / (double a){return Point(x/a,y/a);}
  Point operator * (Point a){ return Point(x * a.x - y * a.y, x * a.y + y * a.x); }
 
  bool operator < (const Point& p) const{ return !equals(x,p.x)?x<p.x:(!equals(y,p.y)&&y<p.y); }
 
  bool operator == (const Point& p)const{ return fabs(x-p.x) < EPS && fabs(y-p.y) < EPS; }
 
};
 
struct Segment{
  Point p1,p2;
  Segment(Point p1 = Point(),Point p2 = Point()):p1(p1),p2(p2){}
  bool operator <  (const Segment& s) const { return ( p2 == s.p2 ) ? p1 < s.p1 : p2 < s.p2; }
  bool operator == (const Segment& s) const { return ( s.p1 == p1 && s.p2 == p2 ) || ( s.p1 == p2 && s.p2 == p1 ); }
 
};
 
typedef Point Vector;
typedef Segment Line;
typedef vector<Point> Polygon;
 
ostream& operator << (ostream& os,const Point& a){ return os << "(" << a.x << "," << a.y << ")"; }
 
ostream& operator << (ostream& os,const Segment& a){ return os << "( " << a.p1 << " , " << a.p2 << " )"; }
 
double dot(Point a,Point b){ return a.x*b.x + a.y*b.y; }
 
double cross(Point a,Point b){ return a.x*b.y - a.y*b.x; }
 
double norm(Point a){ return a.x*a.x+a.y*a.y; }
 
double abs(Point a){ return sqrt(norm(a)); }
 
//rad は角度をラジアンで持たせること
Point rotate(Point a,double rad){ return Point(cos(rad)*a.x - sin(rad)*a.y,sin(rad)*a.x + cos(rad)*a.y); }
 
// 度をラジアンに変換
double toRad(double agl){ return agl*M_PI/180.0; }
 
// a => prev, b => cur, c=> next
// prev から cur へ行って next へ行く際の角度を求める
double getArg(Point a,Point b,Point c){
  double arg1 = atan2(b.y-a.y,b.x-a.x);
  double arg2 = atan2(c.y-b.y,c.x-b.x);
  double arg = fabs( arg1 - arg2 );
  while( arg > M_PI ) arg -= 2.0 * M_PI;
  return fabs(arg);
}
 
int ccw(Point p0,Point p1,Point p2){
  Point a = p1-p0;
  Point b = p2-p0;
  if(cross(a,b) > EPS)return COUNTER_CLOCKWISE;
  if(cross(a,b) < -EPS)return CLOCKWISE;
  if(dot(a,b) < -EPS)return ONLINE_BACK;
  if(norm(a) < norm(b))return ONLINE_FRONT;
  return ON_SEGMENT;
}
 
bool intersectLL(Line l, Line m) {
  return abs(cross(l.p2-l.p1, m.p2-m.p1)) > EPS || // non-parallel
         abs(cross(l.p2-l.p1, m.p1-l.p1)) < EPS;   // same line
}
bool intersectLS(Line l, Line s) {
  return cross(l.p2-l.p1, s.p1-l.p1)*       // s[0] is left of l
         cross(l.p2-l.p1, s.p2-l.p1) < EPS; // s[1] is right of l
}
bool intersectLP(Line l,Point p) {
  return abs(cross(l.p2-p, l.p1-p)) < EPS;
}
bool intersectSS(Line s, Line t) {
  return ccw(s.p1,s.p2,t.p1)*ccw(s.p1,s.p2,t.p2) <= 0 &&
         ccw(t.p1,t.p2,s.p1)*ccw(t.p1,t.p2,s.p2) <= 0;
}
bool intersectSP(Line s, Point p) {
  return abs(s.p1-p)+abs(s.p2-p)-abs(s.p2-s.p1) < EPS; // triangle inequality
}
 
Point projection(Line l,Point p) {
  double t = dot(p-l.p1, l.p1-l.p2) / norm(l.p1-l.p2);
  return l.p1 + (l.p1-l.p2)*t;
}
Point reflection(Line l,Point p) {
  return p + (projection(l, p) - p) * 2;
}
double distanceLP(Line l, Point p) {
  return abs(p - projection(l, p));
}
double distanceLL(Line l, Line m) {
  return intersectLL(l, m) ? 0 : distanceLP(l, m.p1);
}
 
double distanceLS(Line l, Line s) {
  if (intersectLS(l, s)) return 0;
  return min(distanceLP(l, s.p1), distanceLP(l, s.p2));
}
double distanceSP(Line s, Point p) {
  Point r = projection(s, p);
  if (intersectSP(s, r)) return abs(r - p);
  return min(abs(s.p1 - p), abs(s.p2 - p));
}
 
double distanceSS(Line s, Line t) {
  if (intersectSS(s, t)) return 0;
  return min(min(distanceSP(s, t.p1), distanceSP(s, t.p2)),
             min(distanceSP(t, s.p1), distanceSP(t, s.p2)));
}
 
Point crosspoint(Line l,Line m){
  double A = cross(l.p2-l.p1,m.p2-m.p1);
  double B = cross(l.p2-l.p1,l.p2-m.p1);
  if(abs(A) < EPS && abs(B) < EPS){
    vector<Point> vec;
    vec.push_back(l.p1),vec.push_back(l.p2),vec.push_back(m.p1),vec.push_back(m.p2);
    sort(vec.begin(),vec.end());
    assert(vec[1] == vec[2]); //同じセグメントかもよ
    return vec[1];
    //return m.p1;
  }
  if(abs(A) < EPS)assert(false);
  return m.p1 + (m.p2-m.p1)*(B/A);
}
 
Polygon andrewScan(Polygon s)
{
  Polygon u,l;
  if(s.size() < 3)return s;
 
  sort(s.begin(),s.end());
 
  u.push_back(s[0]);
  u.push_back(s[1]);
  l.push_back(s[s.size()-1]);
  l.push_back(s[s.size()-2]);
 
  for(int i=2;i<s.size();i++)
    {
      for(int n=u.size();n >= 2 && ccw(u[n-2],u[n-1],s[i]) != CLOCKWISE; n--)
    u.pop_back();
      u.push_back(s[i]);
    }
 
  for(int i=s.size()-3; i>=0 ; i--)
    {
      for(int n=l.size(); n >= 2 && ccw(l[n-2],l[n-1],s[i]) != CLOCKWISE; n--)
    l.pop_back();
      l.push_back(s[i]);
    }
 
  reverse(l.begin(),l.end());
 
  for(int i = u.size()-2; i >= 1; i--) l.push_back(u[i]);
 
  return l;
}
 
double heron(Point pa,Point pb,Point pc) {
  double a = abs(pa-pb);
  double b = abs(pb-pc);
  double c = abs(pc-pa);
  double s = ( a + b + c ) / 2.0;
  return sqrt( s * ( s - a ) * ( s - b ) * ( s - c ) );
}
 
double LT(double a,double b) { return !equals(a,b) && a < b; }
double LTE(double a,double b) { return equals(a,b) || a < b; }
// Geo <---
 
int n;
 
 
double ternarySearch(int a,int b,vector<Point> &ps){
  int nV = (int)ps.size();
  int L = 1, R = ( a > b ) ? ( ( a - b - 1 + nV ) % nV ) : ( ( ( a + nV ) - b - 1  ) % nV );
  int sp = b;
  int pL = L, pR = R;
 
  rep(_,120){
    if( LT(heron(ps[((L+R*2)/3+sp)%nV],ps[a],ps[b]),heron(ps[((L*2+R)/3+sp)%nV],ps[a],ps[b])) ) {
      R = ( L + R * 2 ) / 3;
    } else {
      L = ( L * 2 + R ) / 3;
    }
    if( R == pR && L == pL ) break;
    pR = R, pL = L;
  }
  return heron(ps[(int)(((L+R)*0.5)+sp)%nV],ps[a],ps[b]);
}
 
bool cmp_x(const Point &p,const Point &q) {
  if( p.x != q.x ) return p.x < q.x;
  return p.y < q.y;
}
 
void compute(vector<Point> &ps) {
  ps = andrewScan(ps);
  double maxi = 0.0;
  n = ps.size();
  if( n == 2 ) assert(false);
  int i = 0, j = 0;
  rep(k,n) {
    if( !cmp_x(ps[i],ps[k]) ) i = k;
    if(  cmp_x(ps[j],ps[k]) ) j = k;
  }
 
  int si = i, sj = j;
  while( i != sj || j != si ) {
    maxi = max(maxi,ternarySearch(i,j,ps)+ternarySearch(j,i,ps));
    if( cross(( ps[(i+1)%n] - ps[i] ),( ps[(j+1)%n]-ps[j])) < 0 ) {
      i = ( i + 1 ) % n;
    } else {
      j = ( j + 1 ) % n;
    }
  }
  printf("%.1f\n",maxi);
}
 
int main() {
  int T;
  while( cin >> T ) {
    while( T-- ) {
      cin >> n;
      vector<Point> ps(n);
      rep(i,n) cin >> ps[i].x >> ps[i].y;
      compute(ps);
    }
  }
  return 0;
}