typedef double db;
const db eps=1e-6;
//点.向量
struct point
{
db x,y,len2,angle;
point() {}
point(db _x,db _y){set(_x,_y);}
void set(db _x,db _y)
{
x=_x,y=_y;
len2=x*x+y*y;
angle=atan2(y,x);
}
void print(){printf("(%.6f,%.6f)",(double)x,(double)y);}
void println(){print();puts("");}
//逆时针旋转
point rotate(db alpha)
{
db _x=x*cos(alpha)-y*sin(alpha);
db _y=x*sin(alpha)+y*cos(alpha);
return point(_x,_y);
}
//向量夹角
db get_angle(point &p)
{
return acos((*this)*p/length()/p.length());
}
db length(){return sqrt(len2);}
point operator+(const point &p){return point(x+p.x,y+p.y);}
point operator-(const point &p){return point(x-p.x,y-p.y);}
point operator*(const db k){return point(k*x,k*y);}
point operator/(const db k){return point(x/k,y/k);}
db operator*(const point &p){return x*p.x+y*p.y;}
db operator^(const point &p){return x*p.y-y*p.x;}
bool operator==(const point &p){return fabs(x-p.x)<eps && fabs(y-p.y)<eps;}
};
//线
struct line
{
point st,ed,vec;
db len2,angle;
line() {}
line(point s,point e) {set(s,e);}
void set(point s,point e)
{
st=s,ed=e,vec=e-s;
len2=vec.len2;
angle=atan2(vec.y,vec.x);
}
//延长k倍
line operator*(db k){return line(st,st+vec*k);}
//靠近st的1/k分点
point operator/(db k){return st+vec/k;}
db length(){return sqrt(len2);}
//点在直线左侧
bool on_left(point p){return (vec^(p-st))>eps;}
//点在直线上
bool on_line(point p){return fabs(vec^(p-st))<eps;}
//点在线段上
bool on_segment(point p)
{
if(p.x+eps<min(st.x,ed.x) || p.x-eps>max(st.x,ed.x))
return false;
if(p.y+eps<min(st.y,ed.y) || p.y-eps>max(st.y,ed.y))
return false;
return on_line(p);
}
//直线平行
bool parallel(line l){return fabs(vec^l.vec)<eps;}
//直线重合
bool coincidence(line l){return on_line(l.st) && on_line(l.ed);}
//直线交点
point get_intersection(line l)
{
db x1=st.x,y1=st.y,x2=ed.x,y2=ed.y;
db x3=l.st.x,y3=l.st.y,x4=l.ed.x,y4=l.ed.y;
db k1=(x4-x3)*(y2-y1),k2=(x2-x1)*(y4-y3);
db x=(k1*x1-k2*x3+(y3-y1)*(x2-x1)*(x4-x3))/(k1-k2);
db y=(k2*y1-k1*y3+(x3-x1)*(y2-y1)*(y4-y3))/(k2-k1);
return point(x,y);
}
//线段相交
bool segment_intersection(line l)
{
if(coincidence(l))
return on_segment(l.st) || on_segment(l.ed)
|| l.on_segment(st) || l.on_segment(ed);
return ((l.st-st)^vec)*((l.ed-st)^vec)<=0
&& ((st-l.st)^l.vec)*((ed-l.st)^l.vec)<=0;
}
//线段中垂线
line vertical_bisector()
{
db x1=st.x,y1=st.y,x2=ed.x,y2=ed.y;
db xm=(x1+x2)/2.0,ym=(y1+y2)/2.0;
return line(point(xm+ym-y1,ym-xm+x1),point(xm-ym+y1,ym+xm-x1));
}
//直线旋转
line rotate(db alpha)
{
point v=vec.rotate(alpha);
return line(st,st+v);
}
//垂线
line vertical(point p)
{
point l(vec.y,-vec.x);
l=l*(distance(p)/l.length());
if(on_line(p+l))return line(p,p+l);
else return line(p,p-l);
}
//直线夹角
db get_angle(line l){return vec.get_angle(l.vec);}
//点到直线的距离
db distance(point p){return fabs((p-st)^(p-ed))/(length());}
//中点
point midpoint(){return point((st.x+ed.x)/2,(st.y+ed.y)/2);}
void print(){st.print();putchar('-');ed.print();}
void println(){print();puts("");}
};
//三角形
struct triangle
{
point v[4];line e[4];
triangle(point p1,point p2,point p3)
{
point p[4]={point(),p1,p2,p3};
set(p);
}
triangle(point p[]) {set(p);}
//构造逆时针三角形(注意特判三点共线)
void set(point p[])
{
for(int i=1;i<=3;++i)v[i]=p[i];
if(!line(v[1],v[2]).on_left(v[3]))
swap(v[2],v[3]);
for(int i=1;i<=2;++i)e[i].set(v[i],v[i+1]);
e[3].set(v[3],v[1]);
}
bool inside(point p)
{
return e[1].on_left(p)==e[2].on_left(p)
&& e[1].on_left(p)==e[3].on_left(p);
}
//重心(中线交点)
point gravity(){return (v[1]+v[2]+v[3])/3;}
//外心
point circum_center()
{
line l1=e[1].vertical_bisector();
line l2=e[2].vertical_bisector();
return l1.get_intersection(l2);
}
//内心
point inscribed_center()
{
db a=e[2].length(),b=e[3].length(),c=e[1].length();
db x=(a*v[1].x+b*v[2].x+c*v[3].x)/(a+b+c);
db y=(a*v[1].y+b*v[2].y+c*v[3].y)/(a+b+c);
return point(x,y);
}
//垂心
point orthocenter()
{
line l1=e[1].vertical(v[3]);
line l2=e[2].vertical(v[1]);
return l1.get_intersection(l2);
}
//周长
db perimeter()
{
db sum=0.0;
for(int i=1;i<=3;++i)sum+=e[i].length();
return sum;
}
//面积
db area(){return fabs((v[3]-v[1])^(v[2]-v[1]))/2;}
void print(){for(int i=1;i<=3;++i){v[i].print();if(i!=3)putchar('-');}}
void println(){print();puts("");}
};
//矩形
struct rectangle
{
point v[5];line e[5];
rectangle(point p1,point p2)
{
point p[5]={point(),point(p1.x,p1.y),point(p1.x,p2.y),
point(p2.x,p1.y),point(p2.x,p2.y)};
set(p);
}
rectangle(point p1,point p2,point p3,point p4)
{
point p[5]={point(),p1,p2,p3,p4};
set(p);
}
rectangle(point p[]) {set(p);}
//构造出顺时针或逆时针的矩形(注意特判四点共线)
void set(point p[])
{
for(int i=1;i<=4;++i)v[i].set(p[i].x,p[i].y);
if(!line(v[1],v[2]).parallel(line(v[3],v[4])))
swap(v[2],v[3]);
if(!line(v[1],v[4]).parallel(line(v[2],v[3])))
swap(v[3],v[4]);
for(int i=1;i<=3;++i)e[i].set(v[i],v[i+1]);
e[4].set(v[4],v[1]);
}
//点在矩形内或矩形上
bool inside(point p)
{
bool flag=true;
for(int i=1;i<=4;++i)
{
if(e[i].on_segment(p))return true;
if(e[i].on_left(p)!=e[1].on_left(p))
flag=false;
}
return flag;
}
bool inside(line l){return inside(l.st) && inside(l.ed);}
bool intersection(line l)
{
for(int i=1;i<=4;++i)
if(e[i].segment_intersection(l))
return true;
return false;
}
void print(){for(int i=1;i<=4;++i){v[i].print();if(i!=4)putchar('-');}}
void println(){print();puts("");}
};
//圆
struct circle
{
point c;db r;
circle() {}
circle(point p,db _r){set(p,_r);}
void set(point p,db _r){c=p,r=_r;}
//外接圆
static circle circum_circle(triangle t)
{
point p=t.circum_center();
return circle(p,line(p,t.v[1]).length());
}
//内接圆
static circle inscribed_circle(triangle t)
{
point p=t.inscribed_center();
return circle(p,t.area()*2/t.perimeter());
}
//点在圆内或圆上
bool inside(point p){return line(p,c).len2<=r*r;}
void print(){c.print(),printf(" r:%.2f",(double)r);}
void println(){print();puts("");}
};
struct polygon
{
const static int __=1e3+5;
point v[__];int n;
polygon() {}
polygon(point p[],int n)
{
for(int i=1;i<=n;++i)v[i]=p[i];
v[0]=p[n],v[n+1]=p[1];
this->n=n;
}
db area()
{
db res=0.0;
for(int i=1;i<=n;++i)
res+=v[i]^v[i+1];
return fabs(res)/2;
}
};
球体积交 && 球体积并
const ld pi=acos(-1);
ld pow2(ld x){return x*x;}
ld pow3(ld x){return x*x*x;}
ld dis(ld x1,ld y1,ld z1,ld x2,ld y2,ld z2)
{
return pow2(x1-x2)+pow2(y1-y2)+pow2(z1-z2);
}
ld cos(ld a,ld b,ld c){return (b*b+c*c-a*a)/(2*b*c);}
ld cap(ld r,ld h){return pi*(r*3-h)*h*h/3;}
//2球体积交
ld sphere_intersect(ld x1,ld y1,ld z1,ld r1,ld x2,ld y2,ld z2,ld r2)
{
ld d=dis(x1,y1,z1,x2,y2,z2);
//相离
if(d>=pow2(r1+r2))return 0;
//包含
if(d<=pow2(r1-r2))return pow3(min(r1,r2))*4*pi/3;
//相交
ld h1=r1-r1*cos(r2,r1,sqrt(d)),h2=r2-r2*cos(r1,r2,sqrt(d));
return cap(r1,h1)+cap(r2,h2);
}
//2球体积并
ld sphere_union(ld x1,ld y1,ld z1,ld r1,ld x2,ld y2,ld z2,ld r2)
{
ld d=dis(x1,y1,z1,x2,y2,z2);
//相离
if(d>=pow2(r1+r2))return (pow3(r1)+pow3(r2))*4*pi/3;
//包含
if(d<=pow2(r1-r2))return pow3(max(r1,r2))*4*pi/3;
//相交
ld h1=r1+r1*cos(r2,r1,sqrt(d)),h2=r2+r2*cos(r1,r2,sqrt(d));
return cap(r1,h1)+cap(r2,h2);
}
网友评论