算法训练 Pollution Solution

问题描述
　　作为水污染管理部门的一名雇员，你需要监控那些被有意无意倒入河流、湖泊和海洋的污染物。你的其中一项工作就是估计污染物对不同的水生态系统（珊瑚礁、产卵地等等）造成的影响。
!

你计算所使用的模型已经在图1中被说明。海岸线（图1中的水平直线）为x轴，污染源位于原点(0, 0)。污染的蔓延呈半圆形，多边形代表了被波及的生态系统。你需要计算出生态系统被污染的面积，也就是图中深蓝色部分。
输入格式
　　输入文件包含仅包含一组测试数据。
　　每组测试数据第一行为两个整数n (3 <= n <= 100), r (1 <= r <= 1000)，n表示了多边形的顶点个数，r表示了污染区域的半径；
　　接下来n行，每行包含两个整数xi (-1500 <= xi <= 1500), yi (0 <= yi <=1500)，表示每个顶点的坐标，以逆时针顺序给出；
　　数据保证多边形不自交或触及自身，没有顶点会位于圆弧上。
输出格式
　　输出多边形被圆心位于原点、半径为r的半圆覆盖的面积。
　　答案的绝对误差不得超过10^-3。
样例输入
6 10
-8 2
8 2
8 14
0 14
0 6
-8 14
样例输出
101.576437872
数据规模和约定
　　存在约30%的数据，n = 3，r <= 20；
　　存在另外约30%的数据，n <= 10，r <= 100，坐标范围不超过100；
　　存在另外约10%的数据，n <= 100，r <= 150，坐标范围不超过250；
　　存在另外约30%的数据，n <= 100，r <= 1000，数据存在梯度；
　　对于100%的数据，满足题目所示数据范围。
　
思路：
利用叉乘计算任意多边形面积的方法。
对每一对点先判断三种情况：
（1）两个点都在圆内；（2）一个在圆内，一个在圆外；（3）两个点在圆外。
（1）可以直接叉乘计算；（2）找到与圆的交点，分为圆内一点和交点的面积+交点和圆外点形成的扇形面积。 （3）先判断到原点最短距离的点C是否在圆内，若不在圆内直接计算扇形面积；若在圆内，递归计算Area(A,C)+Area(C,B);（第三种情况可能是1.两个点均在圆外而且两点的连线不与半圆相交2.两点均在圆外但是两点的连线与半圆相交）

代码：

#include<iostream>
#include<cmath>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<map>
#include<queue>
using namespace std;
const int N=1e2+5;
const int M=N/2;
int n,r;
int X[N],Y[N];
struct P{
	double x,y;
	double getlength(){
		return sqrt(x*x+y*y);
	}
	bool incircle(){
		return x*x+y*y<=r*r;
	}
	double cross(P &b){
		return x*b.y-y*b.x;
	}
}; 
double getArea(P &a,P &b){
	double degree=a.cross(b)/a.getlength()/b.getlength();
	degree=asin(degree);
	return r*r*degree/2;
}
double cal(P &a,P &b){
	bool in1=a.incircle();
	bool in2=b.incircle();
	if(in1&&in2){                               //情况1
		return a.cross(b)/2;
	} 
	else if(in1!=in2){                          //情况2
		P l=a;
		P r=b;
		P mid;
		for(int i=0;i<40;i++){
			mid=P{(l.x+r.x)/2,(l.y+r.y)/2};
			if(mid.incircle()==in1)
			    l=mid;
			else
			    r=mid; 
		}
		if(in1)
		    return a.cross(mid)/2+getArea(mid,b);
		else
		    return getArea(a,mid)+mid.cross(b)/2; 
	}
	else{                                     //情况3
		P l=a;
		P r=b;
		P mid;
		P midr;
		for(int i=0;i<40;i++){
			mid=P{(l.x+r.x)/2,(l.y+r.y)/2};
			midr=P{(l.x+r.x)/2+(r.x-l.x)*0.0001,(l.y+r.y)/2+(r.y-l.y)*0.0001};
		if(mid.getlength()<midr.getlength())
		    r=mid;
		else
		    l=mid;	
		}
		if(mid.incircle())
		    return cal(a,mid)+cal(mid,b);
		else
		    return getArea(a,b); 
	}
}
int main(){
    cin>>n>>r;
    for(int i=0;i<n;i++)
        cin>>X[i]>>Y[i];
    Y[n]=Y[0];
    X[n]=X[0];
    double ans=0;
    for(int i=0;i<n;i++){
    	P a=P{X[i],Y[i]};
    	P b=P{X[i+1],Y[i+1]};
    	ans+=cal(a,b);
	}
	printf("%lf\n",ans);0
	return 0;
}