计算几何题集

这段时间刷了一群计算几何题，呼！

做计算几何时一定要保证自己做到以下两点：

　　1.下手之前保证自己已经考虑到每一个细节，每一步的作法都要非常清晰。

　　2.敲代码时必须非常冷静，必须保持头脑清醒，保证自己敲下来的代码必然是对的（计算几何的代码几乎无法调试）。

接下来就是这一道道题了：

POJ2653

题目意思是一堆木条从天而降，要找到哪些木条是在最上面的（像一个古老的游戏）。

做法不是很难（计算几何做法都不是很难），直接上O(n²)就行了，考察叉积的应用。

#include <iostream>
#include <cstdio>
using namespace std;
const double eps=1e-10;
struct point
{
    double x,y;
    point(){x=y=0;}
    point(double X,double Y){x=X,y=Y;}
};
point V(point start,point end){return point(end.x-start.x,end.y-start.y);}
double crossP(point a,point b){return a.x*b.y-a.y*b.x;}
struct line
{
    bool istop;
    int prev,next;
    point s,e;
}g[110000];int head;
bool have(line a,line b)
{
    return (crossP(V(a.s,a.e),V(a.s,b.s))*crossP(V(a.s,a.e),V(a.s,b.e))<-eps)
         &&(crossP(V(b.s,b.e),V(b.s,a.s))*crossP(V(b.s,b.e),V(b.s,a.e))<-eps);
}
int main(int argc, char *argv[])
{
    int n;
    while(scanf("%d",&n)&&n)
    {
        head=0;
        for(int i=1;i<=n;i++)
        {
            scanf("%lf%lf%lf%lf",&g[i].s.x,&g[i].s.y,&g[i].e.x,&g[i].e.y);
            g[i].istop=true;g[i].next=0;g[i].prev=head;g[head].next=i;
            for(int j=head;j;j=g[j].prev)
                if(have(g[i],g[j]))
                {
                    g[j].istop=false;
                    g[g[j].prev].next=g[j].next;
                    g[g[j].next].prev=g[j].prev;
                }
            head=i;
        }
        printf("Top sticks:");
        for(int i=1;i<n;i++)
            if(g[i].istop)printf(" %d,",i);
        printf(" %d.
",n);
    }
    return 0;
}

POJ2318

题意是将一个箱子用木板隔成好几格子，然后在箱子里撒一些点，求每个格子里点的个数。

直接枚举每个格子，用叉积判断点的是否在里面即可。

其实可以更快，比方说给点排个序，或者二分位置，可以把效率压下来。

#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
struct point 
{
    int x,y;
    point(){x=y=0;}
    point(int X,int Y){x=X,y=Y;}
}up[5100],down[5100];
inline point V(point s,point e){return point(e.x-s.x,e.y-s.y);}
inline int crossP(point a,point b){return a.x*b.y-a.y*b.x;}
int ans[5100];
int main(int argc, char *argv[])
{
    int n,m;point ul,lr;//upper-left lower-right
    while(scanf("%d",&n)&&n)
    {
        scanf("%d%d%d%d%d",&m,&ul.x,&ul.y,&lr.x,&lr.y);
        memset(ans,0,sizeof ans);
        up[0]=point(ul.x,ul.y);down[0]=point(ul.x,lr.y);
        up[n+1]=point(lr.x,ul.y);down[n+1]=point(lr.x,lr.y);
        for(int i=1;i<=n;i++)
        {
            scanf("%d%d",&up[i].x,&down[i].x);
            up[i].y=ul.y;down[i].y=lr.y;
        }
        while(m--)
        {
            point t;scanf("%d%d",&t.x,&t.y);
            for(int i=0;i<=n;i++)
                if(crossP(V(down[i],up[i]),V(down[i],t))<0&&crossP(V(down[i+1],up[i+1]),V(down[i+1],t))>0)
                {
                    ans[i]++;
                    break;
                }
        }
        for(int i=0;i<=n;i++)printf("%d: %d
",i,ans[i]);
        printf("
");
    }
    return 0;
}

poj1113

好了凸包君出现了。

答案是凸包的周长+半径为l的周长的长度。

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
using namespace std;
const double pi=3.14159265358979323846;
struct point{
    int x,y;
    point(){x=y=0;}
    point(int X,int Y){x=X,y=Y;}
}p[1100];
point V(point s,point e){return point(e.x-s.x,e.y-s.y);}
int crossP(point a,point b){return a.x*b.y-a.y*b.x;}
double dis(point a,point b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
bool cmp1(point a,point b){return a.x==b.x?a.y<b.y:a.x<b.x;}
bool cmp2(point a,point b){return a.x==b.x?a.y>b.y:a.x>b.x;}
point s[1100];int r;
int main(int argc, char *argv[])
{
    int n,l;scanf("%d%d",&n,&l);
    for(int i=1;i<=n;i++)scanf("%d%d",&p[i].x,&p[i].y);
    sort(p+1,p+n+1,cmp1);
    r=0;s[++r]=p[1];s[++r]=p[2];
    double ans=2*pi*l;
    for(int i=3;i<=n;i++)
    {
        while(r>1&&crossP(V(s[r-1],s[r]),V(s[r],p[i]))>0)r--;
        s[++r]=p[i];
    }
    for(int i=1;i<r;i++)ans+=dis(s[i],s[i+1]);
    sort(p+1,p+n+1,cmp2);
    r=0;s[++r]=p[1];s[++r]=p[2];
    for(int i=3;i<=n;i++)
    {
        while(r>1&&crossP(V(s[r-1],s[r]),V(s[r],p[i]))>0)r--;
        s[++r]=p[i];
    }
    for(int i=1;i<r;i++)ans+=dis(s[i],s[i+1]);
    printf("%d
",(int)(ans+0.5));
    return 0;
}

POJ1556

题意是求(0,5)到(10,5)的最短路，只要枚举所有的特殊点，互相连边，SPFA即可。

P.S.这是第一次被奇葩POJ卡… G++输出double不能用%lf要用%f，不然会WA（C++提交不会）。

#include <iostream>
#include <cstdio>
#include <cmath>
#include <queue>
#include <cstring>
using namespace std;
const double eps=1e-10;
int n,pnum;
//point
struct point
{
    double x,y;
    int ord;
    point(){x=y=0;}
    point(double X,double Y,int Ord){x=X,y=Y,ord=Ord;}
};
double dis(point a,point b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
//line
struct line
{
    point s,e;
    line(point S,point E){s=S,e=E;}
};
inline double gety(line l,double x)
{
    return (x-l.s.x)*(l.e.y-l.s.y)/(l.e.x-l.s.x)+l.s.y;
}
//wall
struct wall
{
    double x,y[5];
    point p[5];
}w[200];
int head[600];
//graph
struct edge{int to,next;double c;}g[80000];int gnum=0;
void addedge(int from,int to,double c)
{
    g[++gnum].to=to;g[gnum].c=c;g[gnum].next=head[from];head[from]=gnum;
}
double d[500];bool isin[500];
double SPFA()
{
    queue<int> q;
    for(int i=1;i<=pnum;i++)d[i]=1e10,isin[i]=false;
    d[1]=0;isin[1]=true;q.push(1);
    while(!q.empty())
    {    
        int now=q.front();
        for(int i=head[now];i;i=g[i].next)
            if(d[now]+g[i].c+eps<d[g[i].to])
            {
                d[g[i].to]=d[now]+g[i].c;
                if(!isin[g[i].to])
                {
                    q.push(g[i].to);
                    isin[g[i].to]=true;
                }
            }
        isin[now]=false;q.pop();
    }
    return d[2];
}
//main
bool check(line l)
{
    for(int i=1;i<=n;i++)
    {
        if(w[i].x+eps<l.s.x||w[i].x>l.e.x+eps)continue;
        double y=gety(l,w[i].x);
        if(y+eps<w[i].y[1])return false;
        else if(y-eps<w[i].y[2])continue;
        else if(y+eps<w[i].y[3])return false;
        else if(y-eps<w[i].y[4])continue ;
        else return false;
    }
    return true;
}
void addedge(point s,point e)
{
    if(check(line(s,e)))addedge(s.ord,e.ord,dis(s,e));
}
int main(int argc, char *argv[])
{
    //freopen("1.in","r",stdin);
    //freopen("1.out","w",stdout);
    while(scanf("%d",&n)!=EOF&&n!=-1)
    {
        memset(head,0,sizeof head);
        pnum=2;gnum=0;
        point start(0,5,1),end(10,5,2);
        for(int i=1;i<=n;i++)
        {    
            scanf("%lf",&w[i].x);
            for(int j=1;j<=4;j++)
            {
                scanf("%lf",&w[i].y[j]);
                w[i].p[j]=point(w[i].x,w[i].y[j],++pnum);
            }
        }
        addedge(start,end);
        for(int i=1;i<=n;i++)
            for(int j=1;j<=4;j++)
            {
                addedge(start,w[i].p[j]);
                addedge(w[i].p[j],end);
            }
        for(int i=1;i<=n;i++)
            for(int j=1;j<=4;j++)
                for(int i1=i+1;i1<=n;i1++)
                    for(int j1=1;j1<=4;j1++)
                    addedge(w[i].p[j],w[i1].p[j1]);
        printf("%.2f
",SPFA());//神奇的POJ告诉我如果这里用"%0.2lf
"用G++就过不了 C++能过 
    }
    return 0;
}

POJ1228

题意是判断给定的点所构成的凸包是否是稳定的。

稳定的定义是说，给出点集S,不存在一个点P使得

　　1.S的凸包≠S∪{P}的凸包

　　2.S中的点全部出现在S∪{P}的凸包上。

做这题需要一个结论：一个凸包是稳定的当且仅当这个凸包的每条边上都有数量大于等于3的点。

证明也很简单：如果有一条边上只有两个点，我们可以在凸包外侧离这条线极近的地方取一点P，P就能够使得1、2均成立。

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
using namespace std;
struct point{
    int x,y;
    friend bool operator==(point a,point b){return a.x==b.x&&a.y==b.y;}
    point(){x=y=0;}
    point(int X,int Y){x=X,y=Y;}
}p[1100];
point V(point s,point e){return point(e.x-s.x,e.y-s.y);}
int crossP(point a,point b){return a.x*b.y-a.y*b.x;}
bool cmp1(point a,point b){return a.x==b.x?a.y<b.y:a.x<b.x;}
bool cmp2(point a,point b){return a.x==b.x?a.y>b.y:a.x>b.x;}
bool inaline(point a,point b,point c)
{
    return !crossP(V(a,b),V(b,c));
}
point s[1100];int r;
int main(int argc, char *argv[])
{
    int T;scanf("%d",&T);
    while(T--)
    {
        int n;scanf("%d",&n);
        for(int i=1;i<=n;i++)scanf("%d%d",&p[i].x,&p[i].y);
        if(n<=5){printf("NO
");continue;}
        bool ans=true;
        sort(p+1,p+n+1,cmp1);
        r=0;s[++r]=p[1];s[++r]=p[2];
        for(int i=3;i<=n;i++)
        {
            while(r>1&&crossP(V(s[r-1],s[r]),V(s[r],p[i]))>0)r--;
            s[++r]=p[i];
        }
        if(r<=2||r>=n){printf("NO
");continue;}
        for(int i=3;i<=r-2;i++)
            if(!(inaline(s[i-1],s[i],s[i+1])||(inaline(s[i-2],s[i-1],s[i])&&inaline(s[i],s[i+1],s[i+2]))))ans=false;
        if(!ans){printf("NO
");continue;}    
        point em1=s[r-1],em2=s[r-2],sa1=s[2],sa2=s[3];
        sort(p+1,p+n+1,cmp2);
        r=0;s[++r]=p[1];s[++r]=p[2];
        for(int i=3;i<=n;i++)
        {
            while(r>1&&crossP(V(s[r-1],s[r]),V(s[r],p[i]))>0)r--;
            s[++r]=p[i];
        }
        if(r<=2||r>=n){printf("NO
");continue;}
        for(int i=3;i<=r-2;i++)
            if(!(inaline(s[i-1],s[i],s[i+1])||(inaline(s[i-2],s[i-1],s[i])&&inaline(s[i],s[i+1],s[i+2]))))ans=false;
        if(!ans){printf("NO
");continue;}
        if(!(inaline(s[r-1],s[r],sa1)||(inaline(s[r],sa1,sa2)&&inaline(s[r-2],s[r-1],s[r]))))ans=false;
        if(!(inaline(em1,s[1],s[2])||(inaline(em2,em1,s[1])&&inaline(s[1],s[2],s[3]))))ans=false;
        printf(ans?"YES
":"NO
");
    }
    return 0;
}

POJ1151

矩形面积并，经典线段树。

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
using namespace std;
const double eps=1e-10;
double sa[6100];
struct node{int l,r,havecover;double len;}t[6100];
void construct(int l,int r,int k)
{
    t[k].l=l;
    t[k].r=r;
    t[k].havecover=t[k].len=0;
    if(l<r)
    {
        int mid=(l+r)/2;
        construct(l,mid,2*k);
        construct(mid+1,r,2*k+1);
    }
}
void add(int l,int r,int k,int num)
{
    if(t[k].l==l&&t[k].r==r)
    {
        t[k].havecover+=num;
        if(t[k].havecover)t[k].len=sa[r+1]-sa[l];
        else t[k].len=t[k].l<t[k].r?t[2*k].len+t[2*k+1].len:0;
        return ;
    }
    int mid=(t[k].l+t[k].r)/2;
    if(l>mid)add(l,r,2*k+1,num);
    else if(r<=mid)add(l,r,2*k,num);
    else {
        add(l,mid,2*k,num);
        add(mid+1,r,2*k+1,num);
    }
    if(!t[k].havecover)t[k].len=t[2*k].len+t[2*k+1].len;
}
struct line{
    double x;
    int ys,ye,isleft;
    friend bool operator<(line a,line b){return a.x<b.x;}
}l[6100];int lsum;
struct isort
{
    int pos;
    double num;
    bool isstart;
    friend bool operator<(isort a,isort b){return a.num<b.num;}
}s[6100];int snum;
int main(int argc, char *argv[])
{
    int n,time=0;
    while(scanf("%d",&n)&&n)
    {
        lsum=snum=0;
        for(int i=1;i<=n;i++)
        {
            double x1,y1,x2,y2;scanf("%lf%lf%lf%lf",&x1,&y1,&x2,&y2);
            l[++lsum].x=x1;l[lsum].isleft=1;
            s[++snum].num=y1;s[snum].pos=lsum;s[snum].isstart=true;
            s[++snum].num=y2;s[snum].pos=lsum;s[snum].isstart=false;
            l[++lsum].x=x2;l[lsum].isleft=-1;
            s[++snum].num=y1;s[snum].pos=lsum;s[snum].isstart=true;
            s[++snum].num=y2;s[snum].pos=lsum;s[snum].isstart=false;
        }
        sort(s+1,s+snum+1);
        double last=s[1].num;int ord=1;
        sa[ord]=s[1].num;
        for(int i=1;i<=snum;i++)
        {
            if(s[i].num!=last)sa[++ord]=last=s[i].num;
            if(s[i].isstart)l[s[i].pos].ys=ord;
            else l[s[i].pos].ye=ord;
        }
        construct(1,ord-1,1);
        sort(l+1,l+lsum+1);
        double lastx=l[1].x,ans=0;
        cerr<<sa[4]-sa[1]<<endl;
        for(int i=1;i<=lsum;i++)
        {
            //cerr<<" == "<<t[1].len<<" "<<lastx<<" "<<l[i].x<<endl;
            //cerr<<l[i].ys<<" "<<l[i].ye<<endl;
            ans+=t[1].len*(l[i].x-lastx);
            add(l[i].ys,l[i].ye-1,1,l[i].isleft);
            lastx=l[i].x;
            //for(int j=1;j<=10;j++)cerr<<t[j].l<<" "<<t[j].r<<" "<<t[j].havecover<<" "<<t[j].len<<endl;
        }
        printf("Test case #%d
Total explored area: %.2f

",++time,ans);//POJ????%lf 
    }
    return 0;
}

POJ2187

求最远点对，旋转卡壳。

这里有我以前一直理解错的一个地方：旋转卡壳的单调性应该是由对踵点所构成三角形的面积决定的而不是由其距离决定的。

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
using namespace std;
inline int imax(int a,int b){return a>b?a:b;}
struct point
{
    int x,y;
    point(){x=y=0;}
    point(int x,int y){this->x=x;this->y=y;}
};
point V(point s,point e){return point(e.x-s.x,e.y-s.y);}
int crossP(point a,point b){return a.x*b.y-a.y*b.x;}
int dis(point a,point b){return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);}
point p[51000];
bool cmp1(point a,point b){return a.x==b.x?a.y<b.y:a.x<b.x;}
bool cmp2(point a,point b){return a.x==b.x?a.y>b.y:a.x>b.x;}

point c[51000];int cnum;
inline int prev(int num){return num-1?num-1:cnum;}
inline int next(int num){return (num+1<=cnum)?num+1:1;}
point s[51000];int r;
int main(int argc, char *argv[])
{
    freopen("1.in","r",stdin);
    freopen("1.out","w",stdout);
    int n;scanf("%d",&n);
    for(int i=1;i<=n;i++)scanf("%d%d",&p[i].x,&p[i].y);
    
    sort(p+1,p+n+1,cmp1);
    r=0;s[++r]=p[1];s[++r]=p[2];
    for(int i=3;i<=n;i++)
    {
        while(r>1&&crossP(V(s[r-1],s[r]),V(s[r],p[i]))>=0)r--;
        s[++r]=p[i];
    }
    for(int i=1;i<r;i++)c[++cnum]=s[i];
    sort(p+1,p+n+1,cmp2);
    r=0;s[++r]=p[1];s[++r]=p[2];
    for(int i=3;i<=n;i++)
    {
        while(r>1&&crossP(V(s[r-1],s[r]),V(s[r],p[i]))>=0)r--;
        s[++r]=p[i];
    }
    for(int i=1;i<r;i++)c[++cnum]=s[i];
    int ans=0,d=2;
    /*
    for(int i=1;i<=cnum;i++)
        for(d=i+1;d<=cnum;d++)
        ans=imax(ans,dis(c[i],c[d]));
    */
    for(int i=1;i<=cnum;i++)
    {
        while(abs(crossP(V(c[i],c[next(i)]),V(c[i],c[d])))<abs(crossP(V(c[i],c[next(i)]),V(c[i],c[next(d)]))))
            d=next(d);
        ans=imax(ans,dis(c[i],c[d]));ans=imax(ans,dis(c[next(i)],c[d]));
    }
    printf("%d
",ans);
    return 0;
} 

BZOJ1043: [HAOI2008]下落的圆盘

与圆有关。其实不论是并还是交，不论是周长还是面积，写法都是几乎完全一样的。

/**************************************************************
    Problem: 1043
    User: zhuohan123
    Language: C++
    Result: Accepted
    Time:280 ms
    Memory:1392 kb
****************************************************************/
 
#include <iostream>
#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;
const double pi=3.14159265358979323846,eps=1e-10;
inline double angout(double t){return t/pi*180;}
inline double theta(double y,double x)
{
    double ang=atan2(y,x);
    return ang>-eps?ang:2*pi+ang;
}
inline double tabs(double ang)
{
    if(ang<-eps)return ang+2*pi;
    if(ang>2*pi)return ang-2*pi;
    return ang;
}
struct point
{
    double x,y;
    point(){}
    point(double X,double Y){x=X,y=Y;}
};
inline double dis(point a,point b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
inline double sqrdis(point a,point b){return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);}
struct circle{point o;double r;}c[1100];
int n;
struct T
{
    double ang;int isstart;
    friend bool operator<(const T& a,const T& b)
    {
        if(abs(a.ang-b.ang)<eps)
        {
            if(a.isstart==b.isstart)return false;
            else if(a.isstart==0)return true;
            else if(a.isstart==1)return true;
            return false;
        }
        else return a.ang+eps<b.ang;
    }
     
}a[7000];int anum;
int main(int argc, char *argv[])
{
    //freopen("1.in","r",stdin);
    //freopen("1.out","w",stdout);
    scanf("%d",&n);
    for(int i=1;i<=n;i++)scanf("%lf%lf%lf",&c[i].r,&c[i].o.x,&c[i].o.y);
    double ans=0;
    for(int i=1;i<=n;i++)
    {
        bool isallcover=false;
        for(int j=i+1;j<=n;j++)
            if(dis(c[i].o,c[j].o)+c[i].r<c[j].r+eps){isallcover=true;break;}
        if(isallcover)continue;
        anum=0;
        for(int j=i+1;j<=n;j++)
        {
            double d=dis(c[i].o,c[j].o);
            if(d+eps>c[i].r+c[j].r)continue;
            if(d+c[j].r<c[i].r+eps)continue;
            double ang1=acos((c[i].r*c[i].r+sqrdis(c[i].o,c[j].o)-c[j].r*c[j].r)/(2*c[i].r*d)),
                   ang2=theta(c[j].o.y-c[i].o.y,c[j].o.x-c[i].o.x);
            double eang=tabs(ang1+ang2),sang=tabs(ang2-ang1);
            if(sang>eang+eps)
            {
                a[++anum].ang=0;a[anum].isstart=1;
                a[++anum].ang=eang;a[anum].isstart=-1;
                a[++anum].ang=sang;a[anum].isstart=1;
                a[++anum].ang=2*pi;a[anum].isstart=-1;
            }
            else
            {
                a[++anum].ang=sang;a[anum].isstart=1;
                a[++anum].ang=eang;a[anum].isstart=-1;
            }
        }
        a[++anum].ang=0;a[anum].isstart=0;
        sort(a+1,a+anum+1);
        a[++anum].ang=2*pi;a[anum].isstart=0;
        int k=0;
        for(int j=1;j<anum;j++)
        {
            k+=a[j].isstart;
            if(k==0)ans+=c[i].r*(a[j+1].ang-a[j].ang);
        }
    }
    printf("%.3lf
",ans);
    return 0;
}

POJ1981

这道题很有意思，求一个单位圆能够最多能够覆盖多少点。

做法很巧妙。分别以每个点为圆心，1为半径画圆，对于任意两个点，在这两个圆相交的区域中的任意点为圆心所画的单位圆必然同时包含这两个点。

所以，对于每一个圆，让它与其他所有的圆求交，看哪一段被其他圆覆盖的次数最多，更新答案。

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
using namespace std;
inline int imax(int a,int b){return a>b?a:b;}
const double pi=3.14159265358979323846,eps=1e-10;
struct point{
    //circle: R=1
    double x,y;
    point(){}
    point(double X,double Y){x=X,y=Y;}
}c[310];
inline double tabs(double ang){
    if(ang>2*pi+eps)return ang-2*pi;
    if(ang<-eps)return ang+2*pi;
    return ang;
}
inline double dis(point a,point b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
struct angle{
    double ang;int bra;
    angle(){ang=0,bra=0;}
    angle(double Ang,int Bra){ang=Ang,bra=Bra;}
    friend bool operator<(angle a,angle b)
    {
        if(abs(a.ang-b.ang)<eps)return a.bra>b.bra;
        return a.ang+eps<b.ang;
    }
}a[3100];int anum;
int main(int argc, char *argv[])
{
    int n;
    while(scanf("%d",&n)&&n)
    {
        for(int i=1;i<=n;i++)scanf("%lf%lf",&c[i].x,&c[i].y);
        int ans=0;
        for(int i=1;i<=n;i++)
        {
            anum=0;
            a[++anum]=angle(0,0);
            for(int j=1;j<=n;j++)
            {
                double d=dis(c[i],c[j]);
                if(d<eps||d>2+eps)continue ;
                double ang1=tabs(atan2(c[j].y-c[i].y,c[j].x-c[i].x)),
                       ang2=tabs(acos(d/2));
                   double angs=tabs(ang1-ang2),ange=tabs(ang1+ang2);
                   if(angs<ange+eps)
                   {
                       a[++anum]=angle(angs,1);
                       a[++anum]=angle(ange,-1);
                   }
                   else 
                   {
                       a[++anum]=angle(0,1);
                       a[++anum]=angle(ange,-1);
                       a[++anum]=angle(angs,1);
                       a[++anum]=angle(2*pi,-1);
                   }
            }
            sort(a+2,a+anum+1);
            a[++anum]=angle(2*pi,0);
            int k=0;
            for(int i=1;i<anum;i++)
            {
                k+=a[i].bra;
                ans=imax(k+1,ans);
            }
        }
        printf("%d
",ans);
    }
    return 0;
}

BZOJ2178: 圆的面积并

做法如题…

/**************************************************************
    Problem: 2178
    User: zhuohan123
    Language: C++
    Result: Accepted
    Time:780 ms
    Memory:1408 kb
****************************************************************/
 
#include <iostream>
#include <cstdio>
#include <cmath>
#include <algorithm>
//#define ONLINE_JUDGE
using namespace std;
const double pi=3.14159265358979323846,eps=1e-10;
double tabs(double ang)
{
    if(ang>2*pi+eps)return ang-2*pi;
    if(ang<-eps)return ang+2*pi;
    return ang;
}
struct point
{
    double x,y;
    point(){x=y=0;}
    point(double X,double Y){x=X,y=Y;}
};
inline double crossP(point a,point b){return a.x*b.y-a.y*b.x;}
inline double dis(point a,point b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
inline double sqrdis(point a,point b){return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);}
struct circle
{
    point o;double r;
    double area(double ang){return (ang-sin(ang))/2*r*r;}
    point getpo(double ang){return point(cos(ang)*r+o.x,sin(ang)*r+o.y);}
};
circle c[1100];
bool cantuse[1100];
struct angle
{
    double ang;int bra;
    angle(){ang=0,bra=0;}
    angle(double Ang,int Bra){ang=Ang,bra=Bra;}
    friend bool operator<(angle a,angle b)
    {
        if(abs(a.ang-b.ang)<eps)return a.bra>b.bra;
        return a.ang<b.ang+eps;
    }
}a[8100];int anum;
int main(int argc, char *argv[])
{
    #ifndef ONLINE_JUDGE
        freopen("1.in","r",stdin);
        freopen("1.out","w",stdout);
    #endif
    int n;scanf("%d",&n);
    for(int i=1;i<=n;i++)scanf("%lf%lf%lf",&c[i].o.x,&c[i].o.y,&c[i].r);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            if(i!=j&&!cantuse[j]&&(dis(c[i].o,c[j].o)+c[i].r<c[j].r+eps))
            {
                cantuse[i]=true;
                break;
            }
    double ans=0;
    for(int i=1;i<=n;i++)
        if(!cantuse[i])
        {
            anum=0;
            a[++anum]=angle(0,0);
            for(int j=1;j<=n;j++)
                if(!cantuse[j]&&i!=j)
                {
                    double d=dis(c[i].o,c[j].o);
                    if(d+eps>c[i].r+c[j].r)continue ;
                    double ang1=tabs(atan2(c[j].o.y-c[i].o.y,c[j].o.x-c[i].o.x)),
                           ang2=acos((c[i].r*c[i].r+sqrdis(c[i].o,c[j].o)-c[j].r*c[j].r)/(2*c[i].r*d));
                    double angs=tabs(ang1-ang2),ange=tabs(ang1+ang2);
                    if(angs>ange+eps)
                    {
                        a[++anum]=angle(0,1);
                        a[++anum]=angle(ange,-1);
                        a[++anum]=angle(angs,1);
                        a[++anum]=angle(2*pi,-1);
                    }
                    else
                    {
                        a[++anum]=angle(angs,1);
                        a[++anum]=angle(ange,-1);
                    }
                }
            sort(a+2,a+anum+1);
            a[++anum]=angle(2*pi,0);
            int k=0;
            for(int j=1;j<anum;j++)
            {
                k+=a[j].bra;
                if(!k)
                {
                    ans+=c[i].area(a[j+1].ang-a[j].ang);
                    ans+=crossP(c[i].getpo(a[j].ang),c[i].getpo(a[j+1].ang))/2;
                }
            }
        }
    printf("%.3lf
",ans);
    return 0;
}

POJ1279

题意是求多边形核的面积，直接上半平面交。

这里用的是ZZY的排序增量法，效率O(nlgn)，这个算法像极了水平可见直线的做法，只不过是用双端队列而不是栈来维护所有的半平面。

一开始我曾经尝试用ax+by+c≤0的表示法来弄半平面交，一下子就放弃了…计算几何还是必须上向量啊。

#include <iostream>
#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;
const double eps=1e-10;
struct point
{
    double x,y;
    point(){x=y=0;}
    point(double X,double Y){x=X,y=Y;}
    friend bool operator==(point a,point b){return abs(a.x-b.x)<eps&&abs(a.y-b.y)<eps;}
    friend bool operator!=(point a,point b){return !(a==b);}
};
inline point V(point s,point e){return point(e.x-s.x,e.y-s.y);}
inline double crossP(point a,point b){return a.x*b.y-a.y*b.x;}
double area(point p[],int size)//逆时针
{
    double ans=0;
    for(int i=1;i<size;i++)ans+=crossP(p[i],p[i+1]);
    ans+=crossP(p[size],p[1]);
    return ans/2;
}
struct line
{
    //halfplane:(在向量s->e的左侧)
    point s,e;
    double theta;
    line(){}
    line(point S,point E){s=S;e=E;theta=atan2(e.y-s.y,e.x-s.x);}
    friend bool operator<(line a,line b)
    {
        if(abs(a.theta-b.theta)<eps)return crossP(V(b.s,a.e),V(b.s,b.e))<-eps;
        else return a.theta<b.theta;
    }
    bool isout(point p){return crossP(V(s,e),V(s,p))<-eps;}
};
inline point intersection(line a,line b)
{
    double s1=crossP(V(a.s,b.e),V(a.s,b.s)),s2=crossP(V(a.e,b.s),V(a.e,b.e));
    return point((a.s.x*s2+a.e.x*s1)/(s2+s1),(a.s.y*s2+a.e.y*s1)/(s2+s1));
}
line h[5000];int hnum;
void addhalfplane(point s,point e){h[++hnum]=line(s,e);}
point ans[5000];int pnum;
line q[5000];int ql=1,qr=0;
void halfplane_intersection()
{
    int usenum=1;pnum=0;
    sort(h+1,h+hnum+1);
    for(int i=2;i<=hnum;i++)
        if(abs(h[i].theta-h[i-1].theta)>eps)h[++usenum]=h[i];
    //for(int i=1;i<=usenum;i++)cerr<<"("<<h[i].s.x<<","<<h[i].s.y<<") -> ("<<h[i].e.x<<","<<h[i].e.y<<") "<<(h[i].theta/3.1415926535)*180<<endl;
    ql=1;qr=0;
    q[++qr]=h[1];q[++qr]=h[2];
    for(int i=3;i<=usenum;i++)
    {
        while(ql<qr&&h[i].isout(intersection(q[qr],q[qr-1])))qr--;
        while(ql<qr&&h[i].isout(intersection(q[ql],q[ql+1])))ql++;
        q[++qr]=h[i];
    }
    while(ql<qr&&q[ql].isout(intersection(q[qr],q[qr-1])))qr--;
    while(ql<qr&&q[qr].isout(intersection(q[ql],q[ql+1])))ql++;
    if(qr-ql<=1)return ;
    int tpnum=0;
    for(int i=ql;i<qr;i++)ans[++tpnum]=intersection(q[i],q[i+1]);
    ans[++tpnum]=intersection(q[ql],q[qr]);
    for(int i=1;i<=tpnum;i++)
    //printf("(%lf,%lf)
",ans[i].x,ans[i].y);
    pnum=1;
    for(int i=2;i<=tpnum;i++)
        if(ans[i]!=ans[i-1])ans[++pnum]=ans[i];
    return ;
}
point p[5000];

int main(int argc, char *argv[])
{
    //freopen("D.DAT","r",stdin);
    //freopen("1.out","w",stdout);
    int T;scanf("%d",&T);
    while(T--)
    {
        hnum=0;
        int n;scanf("%d",&n);
        for(int i=1;i<=n;i++)scanf("%lf%lf",&p[i].x,&p[i].y);
        if(area(p,n)<-eps)
        {
            for(int i=1;i<n;i++)addhalfplane(p[i+1],p[i]);
            addhalfplane(p[1],p[n]);
        }
        else 
        {
            for(int i=1;i<n;i++)addhalfplane(p[i],p[i+1]);
            addhalfplane(p[n],p[1]);
        }
        halfplane_intersection();
        printf("%.2lf
",area(ans,pnum));
    }
    return 0;
}

BZOJ1038: [ZJOI2008]瞭望塔

半平面交的经典题。山的每一个面可以被看到当且仅当瞭望塔那个面所对应的半平面之内。

/**************************************************************
    Problem: 1038
    User: zhuohan123
    Language: C++
    Result: Accepted
    Time:40 ms
    Memory:4360 kb
****************************************************************/
 
#include <iostream>
#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;
const double eps=1e-10,limit=1e13;
inline double imin(double a,double b){return a<b?a:b;}
struct point
{
    double x,y;
    point(){x=y=0;}
    point(double X,double Y){x=X,y=Y;}
};
point V(point s,point e){return point(e.x-s.x,e.y-s.y);}
double crossP(point a,point b){return a.x*b.y-a.y*b.x;}
struct line
{
    //halfplane: s->e left
    point s,e;
    double ang;
    line(){}
    line(point S,point E){s=S;e=E;ang=atan2(e.y-s.y,e.x-s.x);}
    friend bool operator<(line a,line b)
    {
        if(abs(a.ang-b.ang)<eps)return crossP(V(b.s,b.e),V(b.s,a.e))>eps;
        else return a.ang<b.ang+eps;
    }
    bool isin(point po){return crossP(V(s,e),V(s,po))>-eps;}
    double gety(double x){return (x-s.x)*((e.y-s.y)/(e.x-s.x))+s.y;}
};
point intersection(line a,line b)
{
    double s1=crossP(V(b.s,b.e),V(b.s,a.s)),s2=crossP(V(b.s,a.e),V(b.s,b.e));
    return point((s2*a.s.x+s1*a.e.x)/(s1+s2),(s2*a.s.y+s1*a.e.y)/(s1+s2));
}
line h[5100];int hnum;
void addhalfplane(point s,point e){h[++hnum]=line(s,e);}
line q[51000];int ql,qr;
point hans[51000];int ansnum=0;
void halfplane_intersection()
{
    ql=1;qr=0;ansnum=0;
    sort(h+1,h+hnum+1);
    int unum=1;
    for(int i=2;i<=hnum;i++)
        if(abs(h[i].ang-h[i-1].ang)>eps)h[++unum]=h[i];
    q[++qr]=h[1];q[++qr]=h[2];
    for(int i=3;i<=unum;i++)
    {
        while(ql<qr&&!h[i].isin(intersection(q[qr],q[qr-1])))qr--;
        while(ql<qr&&!h[i].isin(intersection(q[ql],q[ql+1])))ql++;
        q[++qr]=h[i];
    }
    while(ql<qr&&!h[ql].isin(intersection(q[qr],q[qr-1])))qr--;
    while(ql<qr&&!h[qr].isin(intersection(q[ql],q[ql+1])))ql++;
    if(qr-ql<2)return ;
    int tansnum=0;
    hans[++tansnum]=intersection(q[ql],q[qr]);
    for(int i=ql;i<qr;i++)hans[++tansnum]=intersection(q[i],q[i+1]);
    ansnum=1;
    for(int i=2;i<=tansnum;i++)
        if(abs(hans[i].x-hans[i-1].x)>eps||abs(hans[i].y-hans[i-1].y)>eps)
        hans[++ansnum]=hans[i];
    return ;
}
point p[5100];
int main(int argc, char *argv[])
{
    //freopen("1.in","r",stdin);
    //freopen("1.out","w",stdout);
    int n;scanf("%d",&n);
    for(int i=1;i<=n;i++)scanf("%lf",&p[i].x);
    for(int i=1;i<=n;i++)scanf("%lf",&p[i].y);
    addhalfplane(point(-limit,limit),point(-limit,-limit));
    addhalfplane(point(-limit,-limit),point(limit,-limit));
    addhalfplane(point(limit,-limit),point(limit,limit));
    addhalfplane(point(limit,limit),point(-limit,limit));
    for(int i=1;i<n;i++)addhalfplane(p[i],p[i+1]);
    halfplane_intersection();
    double ans=1e10;
    for(int i=1,j=1;i<=ansnum;i++)
    {
        if(hans[i].x+eps<p[1].x||hans[i].x+eps>p[n].x+eps)continue ;
        while(j<n&&p[j+1].x<hans[i].x+eps)j++;
        ans=imin(ans,hans[i].y-line(p[j],p[j+1]).gety(hans[i].x));
    }
    for(int i=1,j=1;i<=n;i++)
    {
        while(j<ansnum&&hans[j+1].x<p[i].x+eps)j++;
        ans=imin(ans,line(hans[j],hans[j+1]).gety(p[i].x)-p[i].y);
    }
    printf("%.3lf
",ans);
    return 0;
}

poj1070

能做到这道题完全是意外…

题意很简单，做法很直接，就是很难搞！

考点只有向量的旋转。本来每次旋转的角度应该二分的，像我一样用三角函数是很不靠谱的做法……

这题有几个特殊的点要考虑的：

　　1.可能一开始是朝右转的。

　　2.不能输出-0.000。

#include <iostream>
#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;
const double pi=3.14159265358979,eps=1e-9;
int n;
struct point
{
    double x,y;
    point(){x=y=0;}
    point(double X,double Y){x=X,y=Y;}
    point rotate(point po,double theta)
    {
        double tx=po.x-x,ty=po.y-y;
        double Sin=sin(theta),Cos=cos(theta);
        double ax=tx*Cos-ty*Sin,ay=ty*Cos+tx*Sin;
        return point(x+ax,y+ay);
    }
    friend point operator+(point a,point b){return point(a.x+b.x,a.y+b.y);}
    inline point normal(){return point(-y,x);}
    inline double len(){return sqrt(x*x+y*y);}
}p[20],zx,hig;
inline point V(point s,point e){return point(e.x-s.x,e.y-s.y);}
inline double dotP(point a,point b){return a.x*b.x+a.y*b.y;}
inline double crossP(point a,point b){return a.x*b.y-a.y*b.x;}
inline double dis(point a,point b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
struct line
{
    point s,e;
    line(){}
    line(point S,point E){s=S,e=E;}
    inline bool ison(point a)
    {
        return a.x+eps>s.x&&a.x<e.x+eps&&fabs(crossP(V(s,a),V(a,e)))<eps;
    }
}l[110];
int lnum;
double len[110],k[110];

point getleftintersection(point o,double r,line l,double lpdis)
{
    double ang1=atan2(l.e.y-l.s.y,l.e.x-l.s.x);
    double ang2=asin(lpdis/r);
    double ang=ang1+ang2+pi;
    return point(o.x+r*cos(ang),o.y+r*sin(ang));
}
int goleft(int o)
{
    double ang=2*pi;
    for(int i=1;i<=n;i++)
        if(i!=o)
        {
            double r=dis(p[i],p[o]);
            for(int j=lnum;j;j--)
            {
                double lpdis=dotP(V(l[j].s,p[o]),V(l[j].s,l[j].e).normal())/V(l[j].s,l[j].e).len();
                if(lpdis>r+eps)continue;
                point inte=getleftintersection(p[o],r,l[j],lpdis);
                if(l[j].ison(inte))
                {
                    double tang=2*asin(dis(p[i],inte)/(2*r));
                    if(tang+eps<ang&&tang>eps)ang=tang;
                }
            }
        }
    for(int i=1;i<=n;i++)
        if(i!=o) p[i]=p[o].rotate(p[i],ang);
    zx=p[o].rotate(zx,ang);
    for(int i=n;i>=1;i--)
        if(i!=o)
            for(int j=1;j<=lnum;j++)
                if(l[j].ison(p[i]))
                return i;
}

point getrightintersection(point o,double r,line l,double lpdis)
{
    double ang1=atan2(l.e.y-l.s.y,l.e.x-l.s.x);
    double ang2=asin(lpdis/r);
    double ang=ang1-ang2;
    return point(o.x+r*cos(ang),o.y+r*sin(ang));
}
int goright(int o)
{
    //cerr<<o<<endl;
    double ang=2*pi;
    for(int i=1;i<=n;i++)
        if(i!=o)
        {
            double r=dis(p[i],p[o]);
            for(int j=lnum;j;j--)
            {
                double lpdis=dotP(V(l[j].s,p[o]),V(l[j].s,l[j].e).normal())/V(l[j].s,l[j].e).len();
                if(lpdis>r+eps)continue;
                point inte=getrightintersection(p[o],r,l[j],lpdis);
                if(l[j].ison(inte))
                {
                    double tang=2*asin(dis(p[i],inte)/(2*r));
                    if(tang+eps<ang&&tang>eps)ang=tang;
                    break ;
                }
            }
        }
    for(int i=1;i<=n;i++)
        if(i!=o) p[i]=p[o].rotate(p[i],2*pi-ang);
    zx=p[o].rotate(zx,2*pi-ang);
    for(int i=n;i>=1;i--)
        if(i!=o)
            for(int j=1;j<=lnum;j++)
                if(l[j].ison(p[i]))
                return i;
}

int main(int argc, char *argv[])
{
    //freopen("wheel.in","r",stdin);
    //freopen("wheel.out","w",stdout);
    int T;scanf("%d",&T);
    while(T--)
    {
        scanf("%d",&n);
        for(int i=1;i<=n;i++)scanf("%lf%lf",&p[i].x,&p[i].y);
        scanf("%lf%lf",&zx.x,&zx.y);
        lnum=0;
        do
        {
            lnum++;
            scanf("%lf%lf",&len[lnum],&k[lnum]);
        }
        while(abs(k[lnum])>eps);
        scanf("%lf%lf",&hig.x,&hig.y);
        l[1].e=hig;    
        for(int i=1;i<=lnum;i++)
        {
            double tx=len[i]/sqrt(1+k[i]*k[i]);
            l[i].s=l[i+1].e=point(l[i].e.x-tx,l[i].e.y-k[i]*tx);
        }
        //for(int i=1;i<=lnum;i++)printf("(%.4lf,%.4lf)->(%.4lf,%.4lf)
",l[i].s.x,l[i].s.y,l[i].e.x,l[i].e.y);
        int start[3]={0},ssum=0;
        for(int i=1;i<=n;i++)
            for(int j=1;j<=lnum;j++)
                if(l[j].ison(p[i]))
                    start[++ssum]=i;
        if(ssum==2&&p[start[2]].x<p[start[1]].x)swap(start[1],start[2]);
        if(ssum==1)start[2]=start[1];
        if(p[start[1]].x+eps>zx.x)
        {
            while(p[start[1]].x+eps>zx.x)
                start[1]=goleft(start[1]);            
        }
        else if(p[start[2]].x+eps<zx.x)
        {
            while(p[start[2]].x+eps<zx.x)
                start[2]=goright(start[2]);
        }
        printf("%.3lf %.3lf
",zx.x,zx.y+eps);
    }
    return 0;
}