除草，把最近的10k代码拉出来表演
#include<bits/stdc++.h>
#define deb(x) cerr<< #x <<" = " <<x<<endl
using namespace std;
const int maxn=1010;
const double eps=5e-9;
const double pi=acos(-1);
double form(double x){
while(x>=2*pi)x-=2*pi;
while(x<0)x+=2*pi;
return x;
}
int sgn(double x){return (x>eps)-(x<-eps);}
int sgn(double a,double b){return sgn(a-b);}
double sqr(double x){return x*x;}
namespace Union{
struct P{
double x,y;
P(){}
P(double x,double y):x(x),y(y){}
double len2(){
return sqr(x)+sqr(y);
}
double len(){
return sqrt(len2());
}
void print(){
fprintf(stderr,"(%.3f,%.3f)\n",x,y);
}
P turn90(){return P(-y,x);}
P norm(){return P(x/len(),y/len());}
P rot(double rad){
return P(x*cos(rad)-y*sin(rad),x*sin(rad)+y*cos(rad));
}
void read(){
scanf("%lf%lf",&x,&y);
}
};
bool operator<(P a,P b){
return sgn(a.x,b.x) ? a.x<b.x : a.y<b.y;
}
bool operator==(P a,P b){
return !sgn(a.x-b.x) and !sgn(a.y-b.y);
}
P operator+(P a,P b){
return P(a.x+b.x,a.y+b.y);
}
P operator-(P a,P b){
return P(a.x-b.x,a.y-b.y);
}
P operator*(P a,double b){
return P(a.x*b,a.y*b);
}
P operator/(P a,double b){
return P(a.x/b,a.y/b);
}
double operator^(P a,P b){
return a.x*b.x + a.y*b.y;
}
double operator*(P a,P b){
return a.x*b.y - a.y*b.x;
}
double det(P a,P b,P c){
return (b-a)*(c-a);
}
double dis(P a,P b){
return (b-a).len();
}
double Area(vector<P>poly){
double ans=0;
for(int i=1;i<poly.size();i++)
ans+=(poly[i]-poly[0])*(poly[(i+1)%poly.size()]-poly[0]);
return ans/2;
}
struct L{
P a,b;
L(){}
L(P a,P b):a(a),b(b){}
P v(){return b-a;}
};
bool onLine(P p,L l){
return sgn((l.a-p)*(l.b-p))==0;
}
bool onSeg(P p,L s){
return onLine(p,s) and sgn((s.b-s.a)^(p-s.a))>=0 and sgn((s.a-s.b)^(p-s.b))>=0;
}
bool parallel(L l1,L l2){
return sgn(l1.v()*l2.v())==0;
}
P intersect(L l1,L l2){
double s1=det(l1.a,l1.b,l2.a);
double s2=det(l1.a,l1.b,l2.b);
return (l2.a*s2-l2.b*s1)/(s2-s1);
}
P project(P p,L l){
return l.a+l.v()*((p-l.a)^l.v())/l.v().len2();
}
double dis(P p,L l){
return fabs((p-l.a)*l.v())/l.v().len();
}
struct C{
P o;
double r;
C(){}
C(P _o,double _r):o(_o),r(_r){}
void read(){
o.read();
scanf("%lf",&r);
}
P at(double rad){
return o+P(r,0).rot(rad);
}
double ang(P p){
double t=atan2(p.y-o.y,p.x-o.x);
return form(t);
}
};
bool operator==(C a,C b){
return a.o==b.o && !sgn(a.r-b.r);
}
bool inPoly(P p,vector<P>poly,bool strict){
int cnt=0;
for(int i=0;i<poly.size();i++){
P a=poly[i],b=poly[(i+1)%poly.size()];
if(onLine(p,L(a,b)))
return !strict;
int x=sgn(det(a,p,b));
int y=sgn(a.y-p.y);
int z=sgn(b.y-p.y);
cnt+=(x>0&&y<=0&&z>0);
cnt-=(x<0&&z<=0&&y>0);
}
return cnt;
}
bool inC(P p,C c,bool strict){
double d=dis(p,c.o);
if(strict){
return sgn(d,c.r)<0;
}else{
return sgn(d,c.r)<=0;
}
}
// 求圆与直线的交点
//turn90() P(-y,x)
double lfix(double x){return sgn(x)?x:0;}
bool intersect(C a, L l, P &p1, P &p2) {
double x = ((l.a - a.o)^ (l.b - l.a)),
y = (l.b - l.a).len2(),
d = x * x - y * ((l.a - a.o).len2() - a.r * a.r);
if (sgn(d) < 0) return false;
d = max(d, 0.0);
P p = l.a - ((l.b - l.a) * (x / y)), delta = (l.b - l.a) * (sqrt(d) / y);
p1 = p + delta, p2 = p - delta;
return true;
}
// 求圆与圆的交点，注意调用前要先判定重圆
bool intersect(C a, C b, P &p1, P &p2) {
double s1 = (a.o - b.o).len();
if (sgn(s1 - a.r - b.r) > 0 || sgn(s1 - fabs(a.r - b.r)) < 0) return false;
double s2 = (a.r * a.r - b.r * b.r) / s1;
double aa = (s1 + s2) * 0.5, bb = (s1 - s2) * 0.5;
P o = (b.o - a.o) * (aa / (aa + bb)) + a.o;
P delta = (b.o - a.o).norm().turn90() * sqrt(lfix(a.r * a.r - aa * aa));
p1 = o + delta, p2 = o - delta;
return true;
}
// 求点到圆的切点，按关于点的顺时针方向返回两个点
bool tang(const C &c, const P &p0, P &p1, P &p2) {
double x = (p0 - c.o).len2(), d = x - c.r * c.r;
if (d < eps) return false; // 点在圆上认为没有切点
P p = (p0 - c.o) * (c.r * c.r / x);
P delta = ((p0 - c.o) * (-c.r * sqrt(d) / x)).turn90();
p1 = c.o + p + delta;
p2 = c.o + p - delta;
return true;
}
// 求圆到圆的外共切线，按关于 c1.o 的顺时针方向返回两条线
vector<L> extan(const C &c1, const C &c2) {
vector<L> ret;
if (sgn(c1.r - c2.r) == 0) {
P dir = c2.o - c1.o;
dir = (dir * (c1.r / dir.len())).turn90();
ret.push_back(L(c1.o + dir, c2.o + dir));
ret.push_back(L(c1.o - dir, c2.o - dir));
} else {
P p = (c1.o * -c2.r + c2.o * c1.r) / (c1.r - c2.r);
P p1, p2, q1, q2;
if (tang(c1, p, p1, p2) && tang(c2, p, q1, q2)) {
if (c1.r < c2.r) swap(p1, p2), swap(q1, q2);
ret.push_back(L(p1, q1));
ret.push_back(L(p2, q2));
}
}
return ret;
}
// 求圆到圆的内共切线，按关于 c1.o 的顺时针方向返回两条线
vector<L> intan(const C &c1, const C &c2) {
vector<L> ret;
P p = (c1.o * c2.r + c2.o * c1.r) / (c1.r + c2.r);
P p1, p2, q1, q2;
if (tang(c1, p, p1, p2) && tang(c2, p, q1, q2)) { // 两圆相切认为没有切线
ret.push_back(L(p1, q1));
ret.push_back(L(p2, q2));
}
return ret;
}
C c[maxn];
vector<P> poly[maxn];
int n,m;
vector<P> intersect(vector<P>poly,C c){
vector<P>ans;
for(int i=0;i<poly.size();i++){
L l=L(poly[i],poly[(i+1)%poly.size()]);
if(!sgn(l.v().len()))
continue;
P p1,p2;
if(intersect(c,l,p1,p2)){
if(onSeg(p1,l)){
ans.push_back(p1);
//cerr<<"ins"<<endl;
//p1.print();
}
if(onSeg(p2,l))
ans.push_back(p2);
}
}
return ans;
}
double calcCir(C cir){
vector<double>ang;
ang.push_back(0);
ang.push_back(pi);
double ans=0;
for(int i=1;i<=n;i++){
if(cir==c[i])continue;
P p1,p2;
if(intersect(cir,c[i],p1,p2)){
ang.push_back(form(cir.ang(p1)));
ang.push_back(form(cir.ang(p2)));
}
}
for(int i=1;i<=m;i++){
vector<P>tmp;
tmp=intersect(poly[i],cir);
//deb(tmp.size());
for(int j=0;j<tmp.size();j++){
// tmp[j].print();
ang.push_back(form(cir.ang(tmp[j])));
}
}
sort(ang.begin(),ang.end());
for(int i=0;i<ang.size();i++){
double t1=ang[i],t2=(i+1==ang.size()?ang[0]+2*pi:ang[i+1]);
P p=cir.at((t1+t2)/2);
int ok=1;
for(int j=1;j<=n;j++){
if(cir==c[j])continue;
if(inC(p,c[j],true)){
ok=0;
break;
}
}
for(int j=1;j<=m&&ok;j++){
if(inPoly(p,poly[j],true)){
ok=0;
break;
}
}
if(ok){
double r=cir.r,x0=cir.o.x,y0=cir.o.y;
ans+=(r*r*(t2-t1)+r*x0*(sin(t2)-sin(t1))-r*y0*(cos(t2)-cos(t1)))/2;
/*cerr<<"---"<<endl;
cerr<<t1<<" "<<t2<<endl;
cerr<<x0<<" "<<y0<<endl;
cerr<<ans<<endl;*/
}
}
return ans;
}
P st;
bool bySt(P a,P b){
return dis(a,st)<dis(b,st);
}
double calcSeg(L l){
double ans=0;
vector<P>pt;
pt.push_back(l.a);
pt.push_back(l.b);
/*cerr<<"line "<<endl;
l.a.print();
l.b.print();
*/
for(int i=1;i<=n;i++){
P p1,p2;
if(intersect(c[i],l,p1,p2)){
if(onSeg(p1,l))
pt.push_back(p1);
if(onSeg(p2,l))
pt.push_back(p2);
}
}
st=l.a;
sort(pt.begin(),pt.end(),bySt);
for(int i=0;i+1<pt.size();i++){
P p1=pt[i],p2=pt[i+1];
P p=(p1+p2)/2;
int ok=1;
for(int j=1;j<=n;j++){
if(sgn(dis(p,c[j].o),c[j].r)<0){
ok=0;
break;
}
}
if(ok){
double x1=p1.x,y1=p1.y,x2=p2.x,y2=p2.y;
// fprintf(stderr,"(%.4f,%.4f)->(%.4f,%.4f)\n",x1,y1,x2,y2);
double res=(x1*y2-x2*y1)/2;
//deb(res);
ans+=res;
}
}
//cerr<<endl;
return ans;
}
C cc[maxn];
double calc(){
int nn=0;
m=1;
for(int i=1;i<=n;i++){
int ok=1;
for(int j=i+1;j<=n;j++){
if(c[i]==c[j]){
ok=0;
}
}
if(ok)
cc[++nn]=c[i];
}
n=nn;
for(int i=1;i<=n;i++)
c[i]=cc[i];
// deb(poly[1].size());
/*deb(Area(poly[1]));
for(int i=0;i<poly[1].size();i++)
poly[1][i].print();*/
if(Area(poly[1])<0){
reverse(poly[1].begin(),poly[1].end());
}
double ans=0;
for(int i=1;i<=n;i++){
ans+=calcCir(c[i]);
//cerr<<sqr(c[i].r)*pi<<endl;
}
double res=0;
for(int i=0;i<poly[1].size();i++){
P p1=poly[1][i],p2=poly[1][(i+1)%poly[1].size()];
if(p1==p2)continue;
res+=calcSeg(L(p1,p2));
}
// cerr<<res<<endl;
// cerr<<"Area"<<endl;
ans+=res;
return ans;
}
vector<P> convex(vector<P>p){
sort(p.begin(),p.end());
/*deb(p.size());
for(int i=0;i<p.size();i++)
p[i].print();*/
vector<P>ans,S;
for(int i=0;i<p.size();i++){
while(S.size()>=2
&& sgn(det(S[S.size()-2],S.back(),p[i]))<=0)
S.pop_back();
S.push_back(p[i]);
}//dw
ans=S;
S.clear();
for(int i=(int)p.size()-1;i>=0;i--){
while(S.size()>=2
&& sgn(det(S[S.size()-2],S.back(),p[i]))<=0)
S.pop_back();
S.push_back(p[i]);
}//up
for(int i=1;i+1<S.size();i++)
ans.push_back(S[i]);
/*deb(ans.size());
for(int i=0;i<ans.size();i++)
ans[i].print();
exit(0);*/
return ans;
}
}
struct P{
double x,y,z;
P(){}
P(double _x,double _y,double _z):x(_x),y(_y),z(_z){}
void read(){
scanf("%lf%lf%lf",&x,&y,&z);
}
double len(){
return sqrt(x*x+y*y+z*z);
}
}p[maxn];
double r[maxn];
double operator^(P a,P b){
return a.x*b.x+a.y*b.y+a.z*b.z;
}
P operator*(P a,P b){
return P(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x);
}
P operator/(P a,double b){
return P(a.x/b,a.y/b,a.z/b);
}
int main(){
int n,m;
P D;
while(cin>>n>>m){
Union::n=n;
D.read();
if(!n&&!m&&!D.x&&!D.y&&!D.z)break;
for(int i=1;i<=n;i++){
p[i].read();
scanf("%lf",&r[i]);
}
D=D/D.len();
P X,Y,Z=D;
if(sgn((P(1,0,0)*Z).len())){
X=P(1,0,0)*Z;
Y=Z*X;
}else{
X=P(0,1,0)*Z;
Y=Z*X;
}
X=X/X.len();
Y=Y/Y.len();
double rad=acos(Z^P(0,0,-1));
double rat=cos(rad);
for(int i=1;i<=n;i++){
Union::c[i].o.x=(X^p[i]);
Union::c[i].o.y=(Y^p[i]);
Union::c[i].r=r[i];
//Union::c[i].o.print();
//cerr<<Union::c[i].r<<endl;
}
//cerr<<endl;
vector<Union::P>pt,ch;
pt.resize(m);
for(int i=0;i<m;i++){
P p;p.read();
pt[i].x=(X^p);
pt[i].y=(Y^p);
}
ch=Union::convex(pt);
//for(int i=0;i<ch.size();i++)
// ch[i].print();
Union::poly[1]=ch;
// cerr<<Union::calc()<<endl;
printf("%.4f\n",Union::calc()/rat);
}
return 0;
}