既见君子

20160218考题 kao

si kao xiong系列第二题kao

题目详见https://121.42.201.137/problem.pdf

题目大意：

  求一个简单多边形和一个半圆的交的面积，多边形点数<=100，|坐标|<=1500

题解：

辛普森积分（奈何当时不会写QAQ）考虑到一个点如果在多边形内部那么射线法会有奇数个交点，所以我们把这些线段关于某个x的所有y坐标求出来，从第奇数个到第偶数个的这一段是在多边形内部需要加起来的，现在仍然有两三个点要挂精度我实在拯救不了它了，唔虽然lpx和yjq都有O（n）做法，但是因为从来没写过辛普森就练手来了

code：

#include<cstdio>
#include<iostream>
#include<cmath>
#include<algorithm>
using namespace std;
#define eps 1e-8
double f(double x);
int n,r;
double simpson(double l,double r){
        double mid=(l+r)/2;
        return (r-l)/6*(f(l)+4*f(mid)+f(r));
}
int sgn(double v){
        if(fabs(v)<eps) return 0;
        return v>0?1:-1;
}
double area(double a,double b,double ans,double Eps){
        double c=(a+b)/2;
        double L=simpson(a,c);
        double R=simpson(c,b);
        if(fabs(R+L-ans)<=15*Eps) return R+L;
        else{
                return area(a,c,L,Eps/2)+area(c,b,R,Eps/2);
        }
}
struct Point{
        double x,y;
        Point(double a=0,double b=0){
                x=a,y=b;
        }
        double operator *(const Point &w){
                return x*w.x+y*w.y;
        }
        Point operator -(const Point &w){
                return Point(x-w.x,y-w.y);
        }
        bool operator <(const Point &w)const{
                return sgn(y-w.y)<=0;
        }
} p[105],t[105];
bool isparallel(Point a,Point b){
        Point v=b-a;
        return sgn(v.x)==0;
}
double high(Point a,Point b,double x){
        double A=a.y-b.y;
        double B=-(a.x-b.x);
        double C=0-A*a.x-B*a.y;
        return (-C-A*x)/B;
}
bool check(Point a,Point b,double x,double y){
        return sgn(y-min(a.y,b.y))>=0&&sgn(y-max(a.y,b.y))<=0&&sgn(x-min(a.x,b.x))>=0&&sgn(x-max(a.x,b.x))<=0;
}
Point Insect(double x){
        double h=fabs(x);
        if(sgn((double)r*r-h*h)<=0) return Point(x,0);
        return Point(x,sqrt((double)r*r-h*h));
}
double f(double x){
        double ans=0;
        int tcnt=0;
        for(int i=1;i<=n;++i){
                if(isparallel(p[i],p[i%n+1])) continue;
                double y=high(p[i],p[i%n+1],x);
                if(check(p[i],p[i%n+1],x,y)){
                        t[++tcnt]=Point(x,y);
                }
        }
        Point c=Insect(x);
        sort(t+1,t+tcnt+1);
        for(int i=1;i<=tcnt;++i)
                if(sgn(t[i].y-c.y)>=0){
                        tcnt=i-1;
                        break;
                }
        t[++tcnt]=c;
        for(int i=1;i<=tcnt-1;i+=2){
                ans+=t[i+1].y-t[i].y;
        }
        return ans;
}
int main(){
        freopen("kao.in","r",stdin);
        freopen("kao.out","w",stdout);
        scanf("%d%d",&n,&r);
        for(int i=1;i<=n;++i) scanf("%lf%lf",&p[i].x,&p[i].y);
        printf("%.4f\n",area(-r,0,simpson(-r,0),eps)+area(0,r,simpson(0,r),eps));
        return 0;
}