HDU 3982 半平面交+圆和凸多边形面积并

从这里看到的这个题。。很容易想到正解。

http://blog.csdn.net/zxy_snow/article/details/6739561

思路大概一样，就是不知道为什么我被卡精度了，，，

acos精度本来就不好，然后题目还要求输出百分比+五位小数，直接把精度卡了。

反正我写出来的当半径很大的时候误差会非常大，会达到3%左右。哎，查了一个下午，用几何画板模拟。真是恶心死了！

我早知道卡精度就不做了。

不知道“小媛”姐姐是怎么做的，实在不想看这个题了。

感觉我的思维和代码都是挺清晰的。。

尽管wa了，还是贴出来吧。。

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <cmath>
#include <iostream>

#define N 222222
#define SIDE 11111
#define EPS 1e-6
#define BUG system("pause")

const double PI=acos(-1.0);

using namespace std;

struct PO
{
    double x,y,ag;
    bool fg;
    void prt() {printf("%lf     %lf    %d\n",x,y,fg);}
}p[N],o,my,dp[N];

struct LI
{
    PO a,b; double ag;
    void in(double x1,double y1,double x2,double y2)
    {a.x=x1; a.y=y1; b.x=x2; b.y=y2;}
    void prt() {printf("%lf     %lf     %lf     %lf     %lf\n",a.x,a.y,b.x,b.y,ag);}
}li[N],deq[N],qs;

double r,ah;
int n,cnt,tn,m;

inline int dc(double x)
{
    if(x>EPS) return 1;
    else if(x<-EPS) return -1;
    return 0;
}

inline PO operator -(PO a,PO b)
{
    PO c; c.x=a.x-b.x; c.y=a.y-b.y;
    return c;
}

inline PO operator +(PO a,PO b)
{
    PO c; c.x=a.x+b.x; c.y=a.y+b.y;
    return c;
}

inline PO operator *(PO a,double k)
{
    a.x*=k; a.y*=k;
    return a;
}

inline PO operator /(PO a,double k)
{
    a.x/=k; a.y/=k;
    return a;
}

inline double cross(PO a,PO b,PO c)
{
    return (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x);
}

inline double dot(PO a,PO b,PO c)
{
    return (b.x-a.x)*(c.x-a.x)+(b.y-a.y)*(c.y-a.y);
}

inline double getlen(const LI &a)
{
    return sqrt((a.b.x-a.a.x)*(a.b.x-a.a.x)+(a.b.y-a.a.y)*(a.b.y-a.a.y));
}

inline double getlen(const PO &a)
{
    return sqrt(a.x*a.x+a.y*a.y);
}

inline double getdis(const PO &a,const PO &b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

inline double getdisl(const PO &a,const LI &b)
{
    PO t1=a-b.a,t2=a-b.b; 
    return fabs(cross(o,t1,t2))/getlen(b);
}

inline void read()
{
    scanf("%lf%d",&r,&tn);
    double a,b,c,d;
    n=0;
    for(int i=1;i<=tn;i++)
    {
        ++n;
        scanf("%lf%lf%lf%lf",&li[i].a.x,&li[i].a.y,&li[i].b.x,&li[i].b.y);
        if(dc(getdisl(o,li[i])-r)>0) n--; 
    }
    scanf("%lf%lf",&my.x,&my.y);
}

inline void getdir()//规范每条线的方向 
{
    for(int i=1;i<=n;i++)
        if(dc(cross(li[i].a,li[i].b,my))<0) swap(li[i].a,li[i].b);
}

inline void getangle()//计算所有直线的极角 
{
    for(int i=1;i<=n;i++)
        li[i].ag=atan2(li[i].b.y-li[i].a.y,li[i].b.x-li[i].a.x);
    li[n+1].in(-SIDE,-SIDE,SIDE,-SIDE);
    li[n+2].in(SIDE,-SIDE,SIDE,SIDE);
    li[n+3].in(SIDE,SIDE,-SIDE,SIDE);
    li[n+4].in(-SIDE,SIDE,-SIDE,-SIDE);
    for(int i=n+1;i<=n+4;i++) li[i].ag=atan2(li[i].b.y-li[i].a.y,li[i].b.x-li[i].a.x);
    n+=4;
}

inline bool cmp(const LI &a,const LI &b)//极角排序，内侧在前 
{
    if(dc(a.ag-b.ag)==0) return dc(cross(a.a,a.b,b.a))<0;
    return a.ag>b.ag;
}

inline PO getpoint(const LI &a,const LI &b)//直线交点 
{
    double k1=cross(a.a,b.b,b.a),k2=cross(a.b,b.a,b.b);
    PO c;
    c=a.a+(a.b-a.a)*k1/(k1+k2);
    return c;
}

inline void getcut()//半平面交 
{
    sort(li+1,li+1+n,cmp);
    m=1;
    for(int i=2;i<=n;i++)
        if(dc(li[i].ag-li[m].ag)!=0)
            li[++m]=li[i];
    deq[1]=li[1]; deq[2]=li[2];
    int bot=1,top=2;
    for(int i=3;i<=m;i++)
    {
        while(bot<top&&dc(cross(li[i].a,li[i].b,getpoint(deq[top],deq[top-1])))<0) top--;
        while(bot<top&&dc(cross(li[i].a,li[i].b,getpoint(deq[bot],deq[bot+1])))<0) bot++;
        deq[++top]=li[i];
    }
    while(bot<top&&dc(cross(deq[bot].a,deq[bot].b,getpoint(deq[top],deq[top-1])))<0) top--;
    while(bot<top&&dc(cross(deq[top].a,deq[top].b,getpoint(deq[bot],deq[bot+1])))<0) bot++;
    if(bot==top) return;
    cnt=0;
    for(int i=bot;i<top;i++) p[++cnt]=getpoint(deq[i],deq[i+1]);
    if(top-1>bot) p[++cnt]=getpoint(deq[bot],deq[top]);
}

inline double getg(const PO &s,const PO &t)//弓形面积 
{
    double res;
    double cc=dot(o,s,t)/getlen(s)/getlen(t);//夹角的余弦值，不能用正弦 asin的范围是[-π/2, π/2] acos的范围是[0, π] 
    cc=acos(cc);
    if(dc(cross(o,s,t))<0) {cc=2*PI-cc;res=fabs(cross(o,s,t))/2;} //三角形 
    else res=-fabs(cross(o,s,t)/2);
    double la=r*cc; //弧长 
    double ss=r*la/2;//扇形面积 
    return res+ss;
}

inline PO getf(const PO &a,const PO &b)//ab的右旋法向量 
{
    PO zt=b-a,rt;
    rt.x=zt.y; rt.y=-zt.x;
    return rt;
}

inline PO rotate(const PO &a,double sss,double ccc)
{
    PO rt;
    rt.x=a.x*ccc-a.y*sss;
    rt.y=a.x*sss+a.y*ccc;
    return rt;
}

inline void getcpoint(const LI &a,PO &s,PO &t)//直线与圆的交点 
{
    double h=getdisl(o,a);
    double tt=sqrt(r*r-h*h);
    double sss=tt/r,ccc=h/r;
    PO f;
    if(dc(cross(o,a.a,a.b))<0) f=getf(a.b,a.a);//判断顺逆时针旋转 
    else f=getf(a.a,a.b);
    f=f/getlen(a)*r;
    s=rotate(f,-sss,ccc);
    t=rotate(f,sss,ccc);
}

inline bool onseg(const PO &a,const LI &b)//a在直线b上，求a是否在线段b上 
{
    PO s1=b.a-a,s2=b.b-a;
    if(dc(dot(o,s1,s2))<=0) return true;
    return false;
}

inline bool cmp1(const PO &a,const PO &b)
{
    return a.ag<b.ag;
}

inline bool cmp2(const PO &a,const PO &b)//unique的去重函数 
{
    if(dc(a.ag-b.ag)==0) return true;
    return false;
}

inline void cir()
{
    for(int i=1;i<=(cnt>>1);i++) swap(p[i],p[cnt-i+1]);//半平面交是顺时针的，转换为逆时针 
    n=0;
    for(int i=1;i<=cnt;i++)
    {
        int dis=dc(getdis(o,p[i])-r);
        if(dis==-1)//凸多边形在圆内的顶点 
        {
            dp[++n]=p[i];
            dp[n].fg=0;
        }
        else if(dis==0)//在圆上的顶点 
        {
            dp[++n]=p[i];
            dp[n].fg=1;
        }
    }
    p[cnt+1]=p[1];
    PO s,t;LI pl;
    for(int i=1;i<=cnt;i++)//凸多边形和圆的交点 
    {
        pl.a=p[i]; pl.b=p[i+1];
        if(dc(getdisl(o,pl)-r)<0)//在圆内 
        {
            getcpoint(pl,s,t);
            if(onseg(s,pl))
            {
                dp[++n]=s;
                dp[n].fg=1;
            }
            if(onseg(t,pl))
            {
                dp[++n]=t;
                dp[n].fg=1;
            }
        }
    }
    for(int i=1;i<=n;i++) dp[i].ag=atan2(dp[i].y-my.y,dp[i].x-my.x);
    sort(dp+1,dp+1+n,cmp1); 
    n=unique(dp+1,dp+1+n,cmp2)-dp-1;
    dp[n+1]=dp[1];
    ah=0;
    for(int i=1;i<=n;i++)
        if(dp[i].fg&&dp[i+1].fg) ah+=getg(dp[i],dp[i+1]);
}

inline double getarea()
{
    double res=0;
    dp[n+1]=dp[1];
    for(int i=1;i<=n;i++) res+=cross(o,dp[i],dp[i+1]);
    return res/2;
}

inline void go()
{
    if(n==0)
    {
        printf("100.00000%%\n");
        return;
    }
    getdir();
    getangle();
    getcut();
    cir();
    double area=fabs(getarea());
    printf("%.5lf%%\n",(area+ah)/(PI*r*r)*100);
    
}

int main()
{
    int cas; scanf("%d",&cas);
    for(int i=1;i<=cas;i++)
    {
        printf("Case %d: ",i);
        read(),go();
    }
    return 0;
}

（至今写的最长的计算几何题了。。。）