Nice Milk - tobiasの日記

問題文
UVa
http://uva.onlinejudge.org/external/101/10117.html

高さhの牛乳に凸多角形のパンをk回浸す。
牛乳に浸されたパンの面積の最大値を求めよ。

パンの辺との距離がhかつ辺に平行な直線でパンを切り取る。
これをmin(n,k)回繰り返して残るパンの面積の最小値を求める。
(元のパンの面積)-(残ったパンの面積の最小値)を出力する。

↓のコードはUVaでは通ったけどPOJではTLE食らってる。

#include<cmath>
#include<algorithm>
#include<iostream>
#include<vector>
#include<climits>
#include<cfloat>
#include<cstdio>
#define curr(P, i) P[(i) % P.size()]
#define next(P, i) P[(i+1) % P.size()]
#define Line pair<point,point>
     
using namespace std;
     
double EPS = 1e-10;
double PI=3.1415926535;
     
double add(double a, double b){
  if(abs(a+b) < EPS * (abs(a)+abs(b)))return 0;
  return a+b;
}
     
struct point{
  double x, y;
  point(){}
  point(double x,double y) : x(x) , y(y){}
     
  point operator + (point p){
    return point(add(x,p.x), add(y,p.y));
  }
     
  point operator - (point p){
    return point(add(x,-p.x), add(y,-p.y));
  }
     
  point operator * (double d){
    return point(x*d,y*d);
  }
};
     
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 sqrt(a.x*a.x+a.y*a.y);
}
     
int ccw(point a, point b, point c) {
  b = b-a; c = c-a;
  if (cross(b, c) > 0)   return +1;
  if (cross(b, c) < 0)   return -1;
  if (dot(b, c) < 0)     return +2;      
  if (norm(b) < norm(c)) return -2;     
  return 0;
}
     
point intersection_l(point a1, point a2, point b1, point b2) {
  return a1 + (a2 - a1) * (cross(b2 - b1,b1 - a1) / cross(b2 - b1,a2 - a1));
}
   
vector<point> convex_cut(vector<point>pol,point a,point b){
  vector<point>q;
  for (int i=0; i<pol.size();i++) {
    point A=curr(pol,i),B=next(pol, i);
    if (ccw(a,b,A)!=-1)q.push_back(A);
    if (ccw(a,b,A)*ccw(a,b,B)<0)
      q.push_back(intersection_l(A,B,a,b));
  }
  return q;
}
   
double areaPol(vector<point> t){
  double area=0;
  for(int i=0;i<t.size();i++)
    area+=(curr(t,i).x-next(t,i).x)*(curr(t,i).y+next(t,i).y);
       
  return abs(area)/2.0;
}
    
point symmetry(point p1, point p2, point Q){
  double xa,ya,xb,yb,t1,t2;
  point R;
   
  xa=Q.x-p1.x,ya=Q.y-p1.y;
  xb=p2.x-p1.x,yb=p2.y-p1.y;
     
  t1=xa*xb+ya*yb,t2=xb*xb+yb*yb;
     
  R.x=2*(p1.x+xb*t1/t2)-Q.x;
  R.y=2*(p1.y+yb*t1/t2)-Q.y;
  return R;
}
   
Line parallel_line(point a,point b,double h){
    double t = atan((a.y-b.y)/(a.x-b.x));
        
    point a_(a.x + h*cos(t+PI/2), a.y + h*sin(t+PI/2));
    point b_(b.x + h*cos(t+PI/2), b.y + h*sin(t+PI/2));
   
    if(ccw(a,b,a_)==-1)a_=symmetry(a,b,a_),b_=symmetry(a,b,b_);
    
    return make_pair(a_,b_);
}
   
int n;
double h;
vector<Line>L;
     
double func(int k,int j,vector<point> pol){
    
    if(k==0 || j==n)return areaPol(pol);
        
    double res=DBL_MAX;
    for(int i=j;i<n;i++)
        res=min(res,func(k-1,i+1,convex_cut(pol,L[i].first,L[i].second)));
        
    return res;
}
     
int main(void){
     
  int k;
     
  while(cin >> n >> k >> h,n|k){
    vector<point>pol(n);
    L.clear();
      
    for(int i=0;i<n;i++)
      cin >> pol[i].x >> pol[i].y;
    
    for(int i=0;i<n;i++)
        L.push_back(parallel_line(curr(pol,i),next(pol,i),h));
      
    printf("%.2f\n",areaPol(pol)-func(min(n,k),0,pol));
  }
     
  return 0;
}