思路:考试的时候我非常地**，写了圆并，然后还TM写了半平面交和三角剖分，虽然只有30分。。但是看在我写了500行的份上还是挂着吧。。
#include<cstdio>
#include<cmath>
#include<cstring>
#include<iostream>
#include<algorithm>
const double Pi=acos(-);
const double eps=1e-;
int n,m,cnt,tot,num;
bool pd[],check;
int read(){
int t=,f=;char ch=getchar();
while (ch<''||ch>'') {if (ch=='-') f=-;ch=getchar();}
while (''<=ch&&ch<='') {t=t*+ch-'';ch=getchar();}
return t*f;
}
struct node{
double l,r;
}L[];
struct Point{
double x,y;
Point(){}
Point(double x0,double y0):x(x0),y(y0){}
}p[],c[],d[];
struct Line{
Point s,e;
double slop;
Line(){}
Line(Point s0,Point e0):s(s0),e(e0){}
}l[],q[];
struct Circle{
Point p;
double r;
Circle(){}
Circle(double x0,double y0,double r0):p(Point(x0,y0)),r(r0){}
}a[];
bool cmp2(node a,node b){
if (fabs(a.l-b.l)<eps) return a.r<b.r;
else return (a.l<b.l);
}
int sgn(double x){if (x>eps) return ;if (x<-eps) return -;return ;}
double operator *(Point p1,Point p2){return p1.x*p2.y-p1.y*p2.x;}
double operator /(Point p1,Point p2){return p1.x*p2.x+p1.y*p2.y;}
Point operator -(Point p1,Point p2){return Point(p1.x-p2.x,p1.y-p2.y);}
Point operator +(Point p1,Point p2){return Point(p1.x+p2.x,p1.y+p2.y);}
Point operator *(Point p1,double x){return Point(p1.x*x,p1.y*x);}
Point operator /(Point p1,double x){return Point(p1.x/x,p1.y/x);}
double sqr(double x){return x*x;}
double llen(Point p1){return sqrt(sqr(p1.x)+sqr(p1.y)); }
double dis(Point p1,Point p2){return llen(p1-p2);}
double dis(Point p1){return llen(p1);}
Point evector(Point p1){double t=llen(p1);if (t<eps) return Point(,);else return p1/t;}
bool cmp3(Line p1,Line p2){
if (p1.slop!=p2.slop) return p1.slop<p2.slop;
else return (p1.e-p1.s)*(p2.e-p1.s)>;
}
bool cmp(Point p1,Point p2){
double t=(p1-p[])*(p2-p[]);
if (fabs(t)<eps) return dis(p[],p1)<dis(p[],p2);
return t>;
}
bool conclude(Circle p1,Circle p2){
double Len=dis(p1.p,p2.p);
if (fabs(p1.r-p2.r)>=Len) return ;
return ;
}
bool intersect(Circle p1,Circle p2){
double Len=dis(p1.p,p2.p);
if (fabs(p1.r-p2.r)<=Len&&Len<=p1.r+p2.r) return ;
return ;
}
double dist_line(Line p){
double A,B,C,dist;
A=p.s.y-p.e.y;
B=p.s.x-p.e.x;
C=p.s.x*p.e.y-p.s.y*p.e.x;
dist=fabs(C)/sqrt(sqr(A)+sqr(B));
return dist;
}
double dist_line(Point p1,Line p){
p.s.x-=p1.x;
p.e.x-=p1.x;
p.s.y-=p1.y;
p.e.y-=p1.y;
return dist_line(p);
}
double get_cos(double a,double b,double c){
return (b*b+c*c-a*a)/(*b*c);
}
double get_angle(Point p1,Point p2){
if (!sgn(dis(p1))||!sgn(dis(p2))) return 0.0;
double A,B,C;
A=dis(p1);
B=dis(p2);
C=dis(p1,p2);
if (C<=eps) return 0.0;
return fabs(acos(get_cos(C,A,B)));
}
Point get_point(Point p){
double T=sqr(p.x)+sqr(p.y);
return Point(sgn(p.x)*sqrt(sqr(p.x)/T),sgn(p.y)*sqrt(sqr(p.y)/T));
}
bool bigconclude(Circle p1,Circle p2,Point p){
Point e1=p2.p-p1.p;
Point e2=p-p1.p;
if (sgn(e1/e2)<) return ;
else return ;
}
double S(Point p1,Point p2,Point p3){
return fabs((p2-p1)*(p3-p1))/;
}
double work(Point p1,Point p2,double R){
Point O=Point(,);
double f=sgn(p1*p2),res=;
if (!sgn(f)||!sgn(dis(p1))||!sgn(dis(p2))) return 0.0;
double l=dist_line(Line(p1,p2));
double a=dis(p1);
double b=dis(p2);
double c=dis(p1,p2);
if (a<=R&&b<=R){
return fabs(p1*p2)/2.0*f;
}
if (a>=R&&b>=R&&l>=R){
double ang=get_angle(p1,p2);
return fabs((ang/(2.0))*(R*R))*f;
}
if (a>=R&&b>=R&&l<=R&&(get_cos(a,b,c)<=||get_cos(b,a,c)<=)){
double ang=get_angle(p1,p2);
return fabs((ang/(2.0))*(R*R))*f;
}
if (a>=R&&b>=R&&l<=R&&(get_cos(a,b,c)>&&get_cos(b,a,c)>)){
double dist=dist_line(Line(p1,p2));
double len=sqrt(sqr(R)-sqr(dist))*2.0;
double ang1=get_angle(p1,p2);
double cos2=get_cos(len,R,R);
res+=fabs(len*dist/2.0);
double ang2=ang1-acos(cos2);
res+=fabs((ang2/())*(R*R));
return res*f;
}
if ((a>=R&&b<R)||(a<R&&b>=R)){
if (b>a) std::swap(a,b),std::swap(p1,p2);
double T=sqr(p1.x-p2.x)+sqr(p1.y-p2.y);
Point u=Point(sgn(p1.x-p2.x)*sqrt(sqr(p1.x-p2.x)/T),sgn(p1.y-p2.y)*sqrt(sqr(p1.y-p2.y)/T));
double dist=dist_line(Line(p1,p2));
double len=sqrt(R*R-dist*dist);
double len2=sqrt(sqr(dis(p2))-sqr(dist));
if (fabs(dis(p2+u*len2)-dist)<=eps) len+=len2;
else len-=len2;
Point p=p2+u*len;
res+=S(O,p2,p);
double ang=get_angle(p1,p);
res+=fabs((ang/2.0)*R*R);
return res*f;
}
return ;
}
Point inter(Line p1,Line p2){
double t1=(p2.e-p1.s)*(p1.e-p1.s);
double t2=(p1.e-p1.s)*(p2.s-p1.s);
if (fabs(t1+t2)<eps) check=;
double k=t1/(t1+t2);
Point p;
p.x=p2.e.x+(p2.s.x-p2.e.x)*k;
p.y=p2.e.y+(p2.s.y-p2.e.y)*k;
return p;
}
bool jud(Line p1,Line p2,Line p3){
Point p=inter(p1,p2);
return (p3.e-p3.s)*(p-p3.s)<;
}
double inter(Circle P){
double res=;p[n+]=p[];
for (int i=;i<=n+;i++)
c[i].x=p[i].x-P.p.x,c[i].y=p[i].y-P.p.y;
for (int i=;i<=n;i++)
res+=work(c[i],c[i+],P.r);
return res;
}
void hpi(){
int L=,R=;tot=;
for (int i=;i<=cnt;i++){
if (l[i].slop!=l[i-].slop) tot++;
l[tot]=l[i];
}
q[++R]=l[];q[++R]=l[];
int all=tot;
for (int i=;i<=all;i++){
while (L<R&&jud(q[R],q[R-],l[i])) R--;
while (L<R&&jud(q[L],q[L+],l[i])) L++;
q[++R]=l[i];
}
while (L<R&&jud(q[R],q[R-],q[L])) R--;
while (L<R&&jud(q[L],q[L+],q[R])) L++;
tot=;
for (int i=L;i<R;i++)
d[++tot]=inter(q[i],q[i+]);
}
bool judge(Point p1,Point p2){
return (fabs(p1.x-p2.x)<eps&&fabs(p1.y-p2.y)<eps);
}
double gworking1(Point c1,Point c2,Point p1,Point p2){
double rev=sgn(c1*c2);
if (fabs(rev)<eps) return ;
num=;
Point O=Point(,);
if ((c1*c2)<) std::swap(c1,c2);
l[++num].s=O,l[num].e=c1;if (judge(l[num].s,l[num].e)) num--;
l[++num].s=c1,l[num].e=c2;if (judge(l[num].s,l[num].e)) num--;
l[++num].s=c2,l[num].e=O;if (judge(l[num].s,l[num].e)) num--;
if ((p1*p2)<) std::swap(p1,p2);
l[++num].s=O,l[num].e=p1;if (judge(l[num].s,l[num].e)) num--;
l[++num].s=p1,l[num].e=p2;if (judge(l[num].s,l[num].e)) num--;
l[++num].s=p2,l[num].e=O;if (judge(l[num].s,l[num].e)) num--;
for (int i=;i<=num;i++)
l[i].slop=atan2(l[i].e.y-l[i].s.y,l[i].e.x-l[i].s.x);
std::sort(l+,l++num,cmp3);
check=;
hpi();
if (check) return ;
d[tot+]=d[];
double res=;
for (int i=;i<=tot;i++)
res+=(d[i]*d[i+])/2.0;
return fabs(res)*rev;
}
double Gworking1(double a1,double a2,Circle w){
double r=w.r;
Point p1=Point(r*cos(a1),r*sin(a1));
Point p2=Point(r*cos(a2),r*sin(a2));
p[n+]=p[];
for (int i=;i<=n+;i++)
c[i]=Point(p[i].x-w.p.x,p[i].y-w.p.y);
double res=;
for (int i=;i<=n;i++)
res+=gworking1(c[i],c[i+],p1,p2);
return res;
}
Point inter(Circle p,Line l,Point p1){
double Dist=dist_line(l);
double dis1=llen(p1);
double len=sqrt(sqr(dis1)-sqr(Dist));
Point e=evector(l.e-l.s);
Point sx=p1;
double R=p.r;
double len2=sqrt(sqr(R)-sqr(Dist));
sx=sx-e*len;
sx=sx+e*len2;
return sx;
}
double gworking2(Point c1,Point c2,Point p1,Point p2,double R){
double rev=sgn(c1*c2),res=;Point O=Point(,);
if ((p1*p2)<) std::swap(p1,p2);
if ((c1*c2)<) std::swap(c1,c2);
if (fabs(rev)<eps) return ;
if (!sgn(dis(p1))||!sgn(dis(p2))) return 0.0;
double l=dist_line(Line(p1,p2));
double a=dis(p1);
double b=dis(p2);
double c=dis(p1,p2);
if (a<=R&&b<=R){
res=fabs(p1*p2)/2.0;
}
else
if (a>=R&&b>=R&&l>=R){
double ang=get_angle(p1,p2);
res=abs((ang/(2.0))*(R*R));
}
else
if (a>=R&&b>=R&&l<=R&&(get_cos(a,b,c)<=||get_cos(b,a,c)<=)){
double ang=get_angle(p1,p2);
res=fabs((ang/(2.0))*(R*R));
}
else
if (a>=R&&b>=R&&l<=R&&(get_cos(a,b,c)>&&get_cos(b,a,c)>)){
double dist=dist_line(Line(p1,p2));
double len=sqrt(sqr(R)-sqr(dist))*2.0;
double ang1=get_angle(p1,p2);
double cos2=get_cos(len,R,R);
res+=fabs(len*dist/2.0);
double ang2=ang1-acos(cos2);
res+=fabs((ang2/())*(R*R));
}
else
if ((a>=R&&b<R)||(a<R&&b>=R)){
if (b>a) std::swap(a,b),std::swap(p1,p2);
double T=sqr(p1.x-p2.x)+sqr(p1.y-p2.y);
Point u=Point(sgn(p1.x-p2.x)*sqrt(sqr(p1.x-p2.x)/T),sgn(p1.y-p2.y)*sqrt(sqr(p1.y-p2.y)/T));
double dist=dist_line(Line(p1,p2));
double len=sqrt(R*R-dist*dist);
double len2=sqrt(sqr(dis(p2))-sqr(dist));
if (fabs(dis(p2+u*len2)-dist)<=eps) len+=len2;
else len-=len2;
Point p=p2+u*len;
res+=S(O,p2,p);
double ang=get_angle(p1,p);
res+=fabs((ang/2.0)*R*R);
}
int in1=((sgn(c1*p1)*sgn(p1*c2))>=);
int in2=((sgn(c1*p2)*sgn(p2*c2))>=);
if (in1&&in2) return res*rev;
if (!in1){
if (dis(p1)<=R){
Point sb=inter(Line(O,c1),Line(p1,p2));
double Di=dis(sb);
if (Di<=R) res-=S(O,sb,p1);
else{
Point sx=inter(Circle(,,R),Line(p1,p2),p1);
double ang=get_angle(c1,sx);
res-=fabs((ang/2.0)*R*R);
res-=S(O,sx,p1);
}
}else{
Point sb=inter(Line(O,c1),Line(p1,p2));
double Di=dis(sb);
if (Di>=R){
double ang=get_angle(c1,p1);
res-=fabs((ang/2.0)*R*R);
}else{
Point sx=inter(Circle(,,R),Line(p1,p2),p1);
double ang2=get_angle(p1,sx);
res-=fabs((ang2/2.0)*R*R);
res-=S(sb,sx,O);
}
}
}
if (!in2){
if (dis(p2)<=R){
Point sb=inter(Line(O,c2),Line(p1,p2));
double Di=dis(sb);
if (Di<=R) res-=S(O,sb,p2);
else{
Point sx=inter(Circle(,,R),Line(p2,p1),p2);
double ang=get_angle(c2,sx);
res-=fabs((ang/2.0)*R*R);
res-=S(O,sx,p2);
}
}else{
Point sb=inter(Line(O,c2),Line(p1,p2));
double Di=dis(sb);
if (Di>=R){
double ang=get_angle(c2,p2);
res-=fabs((ang/2.0)*R*R);
}else{
Point sx=inter(Circle(,,R),Line(p2,p1),p2);
double ang2=get_angle(p2,sx);
res-=fabs((ang2/2.0)*R*R);
res-=S(sb,sx,O);
}
}
}
return res*rev;
}
double Gworking2(double a1,double a2,Circle w){
double res=,r=w.r;
Point p1=Point(r*cos(a1),r*sin(a1));
Point p2=Point(r*cos(a2),r*sin(a2));
p[n+]=p[];
for (int i=;i<=n+;i++)
c[i]=Point(p[i].x-w.p.x,p[i].y-w.p.y);
for (int i=;i<=n;i++)
res+=gworking2(c[i],c[i+],p1,p2,r);
return res;
}
void graham(){
int k=;
for (int i=;i<=n;i++)
if (p[i].y<p[k].y||(p[i].x<p[k].x&&p[i].y==p[k].y)) k=i;
std::sort(p+,p++n,cmp);
int top=;c[top]=p[];
for (int i=;i<=n;i++){
while (top>&&(c[top]-c[top-])*(p[i]-c[top])<=) top--;
c[++top]=p[i];
}
n=top;
for (int i=;i<=top;i++)
p[i]=c[i];
}
void Work(Circle p1,Point a1,Point a2){
Point O=p1.p;double r=p1.r;
a1.x-=O.x;a2.x-=O.x;
a1.y-=O.y;a2.y-=O.y;
double x1=atan2(a1.y,a1.x),x2=atan2(a2.y,a2.x);
if (x1<) x1+=*Pi;
if (x2<) x2+=*Pi;
if (x1<=x2){
cnt++;
L[cnt].l=x1;
L[cnt].r=x2;
}else{
cnt++;
L[cnt].l=x1;
L[cnt].r=*Pi;
cnt++;
L[cnt].l=0.0;
L[cnt].r=x2;
}
}
void add(Circle p1,Circle p2){
double x1=p1.p.x,x2=p2.p.x,y1=p1.p.y,y2=p2.p.y,r1=p1.r,r2=p2.r;
double A=-*(x1+x2),B=-*(y1+y2),C=sqr(x1)+sqr(x2)+sqr(y1)+sqr(y2)-sqr(r1)-sqr(r2);
Point S,E;
if (fabs(A)<eps){
S=Point(,(C)/B);
E=Point(,(C)/B);
}else
if (fabs(B)<eps){
S=Point((C)/A,);
E=Point((C)/A,);
}else{
S=Point(,(C-A)/B);
E=Point(,(C-*A)/B);
}
Line l=Line(S,E);
double Dist1=dist_line(p1.p,l);
double Dist2=llen(p1.p-S);
Point e=evector(E-S);
double Dist3=sqrt(sqr(Dist2)-sqr(Dist1));
e=e*Dist3;
Point a1,a2;
if (fabs(dis(S+e,p1.p)-Dist1)<eps){
S=S+e;
}
else S=S-e;
e=e/Dist3;
double Dist4=sqrt(sqr(p1.r)-sqr(llen(S-p1.p)));
e=e*Dist4;
a1=S+e;
a2=S-e;
if ((a1-p1.p)*(a2-p1.p)<) std::swap(a1,a2);
//a1放在a2的顺时针方向、
if (bigconclude(p1,p2,a1)) std::swap(a1,a2);
//a1到a2必须是逆时针。
Work(p1,a1,a2);
}
int main(){
freopen("aerolite.in","r",stdin);
freopen("aerolite.out","w",stdout);
scanf("%d",&m);
for (int i=;i<=m;i++){
scanf("%lf%lf%lf",&a[i].p.x,&a[i].p.y,&a[i].r);
}
scanf("%d",&n);
for (int i=;i<=n;i++){
scanf("%lf%lf",&p[i].x,&p[i].y);
}
graham();
for (int i=;i<=m;i++) pd[i]=;
for (int i=;i<=m;i++){
for (int j=;j<=m;j++)
if (i!=j){
if (conclude(a[i],a[j])) {pd[i]=;break;}
}
}
double ans=;
for (int i=;i<=m;i++)
if (pd[i]){
bool mark=;
for (int j=;j<=m;j++)
if (i!=j&&pd[j]){
if (intersect(a[i],a[j])) {mark=;break;}
}
if (!mark) ans+=inter(a[i]),pd[i]=;
}
int j=;
for (int i=;i<=m;i++)
if (pd[i]) a[++j]=a[i];
for (int i=;i<=m;i++) pd[i]=;
for (int i=;i<=m;i++){
cnt=;
for (int j=;j<=m;j++)
if (i!=j){
add(a[i],a[j]);
}
std::sort(L+,L++cnt,cmp2);
L[].r=0.0;
L[cnt+].l=*Pi;
double Reach=0.0,Last=0.0;
int ww=;
for (int j=;j<=cnt;j++)
if (L[j].l>Reach){
if (j==) ww++;
ans+=Gworking1(Last,Reach,a[i]);
Last=L[j].l;
ans+=Gworking2(Reach,L[j].l,a[i]);
Reach=L[j].r;
}else{
Reach=L[j].r;
}
if (fabs(Reach-*Pi)<eps){
ans+=Gworking1(Last,Reach,a[i]);
}else
if (fabs(Reach-*Pi)>=eps){
ans+=Gworking1(Last,Reach,a[i]);
ans+=Gworking2(Reach,*Pi,a[i]);
}
}
printf("%.8lf\n",ans);
return ;
}