转载来自csdn资源下载，这里做记录。
#include<cstdlib>
#include<cmath>
#include<cstdio>
#include<algorithm>
#define max(a,b) (((a)>(b))?(a):(b))
#define min(a,b) (((a)>(b))?(b):(a))
#define sign(x) ((x)>eps?1:((x)<-eps?(-1):(0)))
using namespace std;
const int MAXN=1000;
const double eps=1e-8,inf=1e50 ,Pi=acos(-1.0);
struct point {
double x,y;
point() {}
point(double _x,double _y) {
x=_x;
y=_y;
}
point operator-(const point &ne)const {
return point(x-ne.x,y-ne.y);
}
point operator+(const point ne)const {
return point(x+ne.x,y+ne.y);
}
point operator*(const double t)const{
return point(x*t,y*t);
}
point operator/(const double t)const{
if(sign(t)==0)exit(1);
return point(x/t,y/t);
}
};
struct line{
point a,b;
line(){}
line(point _a,point _b){a=_a;b=_b;}
};
struct line2 {
double a,b,c;
line2() {}
line2(double _a,double _b,double _c) {
a=_a;
b=_b;
c=_c;
}
};
struct circle {
point o;
double r;
circle(){}
circle(point _o,double _r){
o=_o;
r=_r;
}
};
struct rectangle{
point a,b,c,d;
rectangle(){}
rectangle(point _a,point _b,point _c,point _d){
a=_a;
b=_b;
c=_c;
d=_d;
}
};
struct polygon{
point p[MAXN];
int n;
};
inline double xmult(point a,point b){
return a.x*b.y-a.y*b.x;
}
inline double xmult(point o,point a,point b) {
return (a.x-o.x)*(b.y-o.y)-(b.x-o.x)*(a.y-o.y);
}
inline double xmult(double x1,double y1,double x2,double y2) {
return x1*y2-x2*y1;
}
inline double dmult(point o,point a,point b) {
return (a.x-o.x)*(b.x-o.x)+(a.y-o.y)*(b.y-o.y);
}
inline double dmult(point a,point b) {
return a.x*b.x+a.y*b.y;
}
inline double lenth(point a){
return sqrt(dmult(a,a));
}
inline double dist(point a,point b){
return lenth(b-a);
}
inline double dist2(point a,point b){
return dmult(b-a,b-a);
}
//直线一般式转两点式
line toline(double a,double b,double c) {
if(sign(b)==0)exit(1);
point A(0,-c/b),B(-c/a,0);
return line(A,B);
}
//直线两点式转一般式
line2 toline2(point a,point b) {
double A=b.y-a.y,B=a.x-b.x,C=-B*a.y-A*a.x;
return line2(A,B,C);
}
//点p绕o逆时针旋转alpha
point rotate(point o,point p,double alpha) {
point tp;
p.x-=o.x;
p.y-=o.y;
tp.x=p.x*cos(alpha)-p.y*sin(alpha)+o.x;
tp.y=p.y*cos(alpha)+p.x*sin(alpha)+o.y;
return tp;
}
//向量u的倾斜角
double angle(point u) {
return atan2(u.y,u.x);
}
//oe与os的夹角，夹角正负满足叉积
double angle(point o,point s,point e) {
point os=s-o,oe=e-o;
double bot=sqrt(dmult(os,os)*dmult(oe,oe));
double top=dmult(os,oe);
double cosfi=top/bot;
if (cosfi >= 1.0 ) return 0;
if (cosfi <= -1.0 ) return -Pi;
double fi=acos(cosfi);
if(xmult(o,s,e)>0)return fi;
else return -fi;
}
//p在l上的投影与l关系
double relation(point p,line l) {
line tl(l.a,p);
return dmult(tl.b-l.a,l.b-l.a)/dist2(l.a,l.b);
}
//p在l上的垂足
point perpendicular(point p,line l) {
double r=relation(p,l);
return l.a+((l.b-l.a)*r);
}
//求点p到线段l的最短距离,并返回线段上距该点最近的点np
double dist_p_to_seg(point p,line l,point &np) {
double r=relation(p,l);
if(r<0) {
np=l.a;
return dist(p,l.a);
}
if(r>1) {
np=l.b;
return dist(p,l.b);
}
np=perpendicular(p,l);
return dist(p,np);
}
//求点p到直线l的最短距离
double dist_p_to_line(point p,line l) {
return abs(xmult(p-l.a,l.b-l.a))/dist(l.a,l.b);
}
//线段之间最短距离
inline double dist_seg_to_seg(line p,line q) {
return min(min(dist_p_to_seg(p.a,q),dist_p_to_seg(p.b,q)),min(dist_p_to_seg(q.a,p),dist_p_to_seg(q.b,p)));
}
//求矢量线段夹角的余弦
double cosine(line u,line v) {
point pu=u.b-u.a,pv=v.b-v.a;
return dmult(pu,pv)/sqrt(dmult(pu,pu)*dmult(pv,pv));
}
//求矢量的夹角的余弦
double cosine(point a,point b){
return dmult(a,b)/sqrt(dmult(a,a)*dmult(b,b));
}
//求线段夹角
double lsangle(line u,line v) {
point o(0,0),a=u.b-u.a,b=v.b-v.a;
return angle(o,a,b);
}
//求直线斜率
double slope(line2 l) {
if(abs(l.a) < 1e-20)
return 0;
if(abs(l.b) < 1e-20)
return inf;
return -(l.a/l.b);
}
//直线倾斜角[0,Pi]
double alpha(line2 l) {
if(abs(l.a)< eps)
return 0;
if(abs(l.b)< eps)
return Pi/2;
double k=slope(l);
if(k>0)
return atan(k);
else
return Pi+atan(k);
}
//点关于直线的对称点
point symmetry(line2 l,point p){
point tp;
tp.x=((l.b*l.b-l.a*l.a)*p.x-2*l.a*l.b*p.y-2*l.a*l.c)/(l.a*l.a+l.b*l.b);
tp.y=((l.a*l.a-l.b*l.b)*p.y-2*l.a*l.b*p.x-2*l.b*l.c)/(l.a*l.a+l.b*l.b);
return tp;
}
//判多边形是否逆时针
bool is_unclock(polygon pg) {
int n=pg.n;
pg.p[n]=pg.p[0];
double area=0;
for(int i=0; i<n; i++)
area+=xmult(pg.p[i].x,pg.p[i].y,pg.p[i+1].x,pg.p[i+1].y);
return area>-eps;
}
//改变多边形时针顺序
void to_unclock(polygon &pg) {
for(int i=0,j=pg.n-1; i<j; i++,j--)
swap(pg.p[i],pg.p[j]);
}
//判定凸多边形,顶点按顺时针或逆时针给出,允许相邻边共线
int is_convex(point p[],int n) {
int i,s[3]= {1,1,1};
for (i=0; i<n&&s[1]|s[2]; i++)
s[(sign(xmult(p[(i+1)%n],p[(i+2)%n],p[i]))+3)%3]=0;
return s[1]|s[2];
}
//判定凸多边形,顶点按顺时针或逆时针给出,不允许相邻边共线
int is_convex_v2(point p[],int n) {
int i,s[3]= {1,1,1};
for (i=0; i<n&&s[0]&&s[1]|s[2]; i++)
s[(sign(xmult(p[(i+1)%n],p[(i+2)%n],p[i]))+3)%3]=0;
return s[0]&&s[1]|s[2];
}
//判点在凸多边形内或多边形边上,顶点按顺时针或逆时针给出
int inside_convex(point q,point p[],int n) {
int i,s[3]= {1,1,1};
for (i=0; i<n&&s[1]|s[2]; i++)
s[(sign(xmult(p[(i+1)%n],q,p[i]))+3)%3]=0;
return s[1]|s[2];
}
//判点在凸多边形内,顶点按顺时针或逆时针给出,在多边形边上返回0
int inside_convex2(point q,point p[],int n) {
int i,s[3]= {1,1,1};
for (i=0; i<n&&s[0]&&s[1]|s[2]; i++)
s[(sign(xmult(p[(i+1)%n],q,p[i]))+3)%3]=0;
return s[0]&&s[1]|s[2];
}
//判点在线段上
inline int p_on_seg(point a,point p1,point p2) {
if(fabs(xmult(a,p1,p2))<=eps&&(a.x-p1.x)*(a.x-p2.x)<eps&&(a.y-p1.y)*(a.y-p2.y)<eps)
return 1;
return 0;
}
//判点在线段端点左方
inline int p_on_segvex(point s,point p) {
return fabs(p.y-s.y)<eps&&(p.x<=s.x+eps);
}
//判线段相交 <=:不规范相交
inline int seg_inter(line s,line p) {
double minx1=min(s.a.x,s.b.x),maxx1=max(s.a.x,s.b.x);
double minx2=min(p.a.x,p.b.x),maxx2=max(p.a.x,p.b.x);
double miny1=min(s.a.y,s.b.y),maxy1=max(s.a.y,s.b.y);
double miny2=min(p.a.y,p.b.y),maxy2=max(p.a.y,p.b.y);
if((minx1>maxx2+eps)||(minx2>maxx1+eps)||
(miny1>maxy2+eps)||(miny2>maxy1+eps))
return 0;
else
return sign(xmult(s.a,s.b,p.a)*xmult(s.a,s.b,p.b))<=0&&
sign(xmult(p.a,p.b,s.a)*xmult(p.a,p.b,s.b))<=0;
}
//判点在多边形内部
inline int p_in_polygon(point a,point p[],int n) {
int count = 0;
line s,ps;
ps.a = a,ps.b = a;
ps.b.x = inf;
for(int i = 0; i < n; i++) {
s.a = p[i];
if(i + 1 < n)s.b = p[i+1];
else s.b = p[0];
if (s.a.y > s.b.y)swap(s.a,s.b);
if (p_on_seg(a,s.a,s.b))return 2;
if ((fabs(s.a.y-s.b.y)>eps)) {
if (p_on_segvex(s.b,a)) count++;
else if (seg_inter(ps,s))count++;
}
}
if (count%2)return 1;
return 0;
}
//多边形内部最长线段
point stk[MAXN];
double seg_max_len(line u,polygon &pg){
double ans=0.0,tmp=inf;
pg.p[pg.n]=pg.p[0];
int n=pg.n,top=0;
for(int i=0;i<n;i++){
line v(pg.p[i],pg.p[i+1]);
int s1=sign(xmult(u.a,u.b,v.a)),s2=sign(xmult(u.a,u.b,v.b));
if(s1*s2<=0&&(s1!=0||s2!=0)){
stk[top++]=line_intersection(u,v);
}
}
stk[top++]=u.a;
stk[top++]=u.b;
sort(stk,stk+top);
top=unique(stk,stk+top)-stk;
point mp,lp=stk[0];
for(int i=1;i<top;i++){
mp=(lp+stk[i])*0.5;
if(!p_in_polygon(mp,pg))lp=stk[i];
ans=max(ans,dist2(lp,stk[i]));
}
return sqrt(ans);
}
double maxlenth(polygon &pg){
double ans=0.0;
for(int i=0;i<pg.n-1;i++){
for(int j=i+1;j<pg.n;j++){
ans=max(ans,seg_max_len(line(pg.p[i],pg.p[j]),pg));
ans=max(ans,seg_max_len(line(pg.p[j],pg.p[i]),pg));
}
}
return ans;
}
//凸包对踵点长度
double opposite_lenth(polygon &pg){
double ans=inf;
int a,b,c;
pg.p[pg.n]=pg.p[0];
for(a=0,b=1,c=2;a<pg.n;a++,b++){
while(getarea(pg.p[a],pg.p[b],pg.p[c])<getarea(pg.p[a],pg.p[b],pg.p[(c+1)%pg.n]))c=(c+1)%pg.n;
ans=min(ans,dist_p_to_line(pg.p[c],line(pg.p[a],pg.p[b])));
}
return ans;
}
//判p1,p2是否在l1,l2两侧
inline int opposite_side(point p1,point p2,point l1,point l2) {
return xmult(l1,l2,p1)*xmult(l1,l2,p2)<-eps;
}
//判线段在任意多边形内,顶点按顺时针或逆时针给出,与边界相交返回1
int seg_in_polygon(point l1,point l2,point p[],int n) {
point t[MAXN],tt;
int i,j,k=0;
if (!p_in_polygon(l1,p,n)||!p_in_polygon(l2,p,n))//不在内部
return 0;
for (i=0; i<n; i++)
if (opposite_side(l1,l2,p[i],p[(i+1)%n])&&opposite_side(p[i],p[(i+1)%n],l1,l2))
return 0;
else if (p_on_seg(l1,p[i],p[(i+1)%n]))
t[k++]=l1;
else if (p_on_seg(l2,p[i],p[(i+1)%n]))
t[k++]=l2;
else if (p_on_seg(p[i],l1,l2))
t[k++]=p[i];
for (i=0; i<k; i++)
for (j=i+1; j<k; j++) {
tt.x=(t[i].x+t[j].x)/2;
tt.y=(t[i].y+t[j].y)/2;
if (!p_in_polygon(tt,p,n))
return 0;
}
return 1;
}
//求直线交点，必须存在交点，或者预判断【解析几何方法】
point line_intersection(line u,line v) {
double a1=u.b.y-u.a.y,b1=u.a.x-u.b.x;
double c1=u.b.y*(-b1)-u.b.x*a1;
double a2=v.b.y-v.a.y,b2=v.a.x-v.b.x;
double c2=v.b.y*(-b2)-v.b.x*a2;
double D=xmult(a1,b1,a2,b2);
return point(xmult(b1,c1,b2,c2)/D,xmult(c1,a1,c2,a2)/D);
}
//求线段交点，必须存在交点，或者预判断【平面几何方法】
point line_intersection2(line u,line v) {
point ret=u.a;
double t=xmult(u.a-v.a,v.b-v.a)/xmult(u.b-u.a,v.b-v.a);
t=fabs(t);
ret.x+=(u.b.x-u.a.x)*t;
ret.y+=(u.b.y-u.a.y)*t;
return ret;
}
//三角形重心
point barycenter(point a,point b,point c){
return (a+b+c)/3.0;
}
//多边形重心
point barycenter(point p[],int n) {
point ret,t;
double t1=0,t2;
int i;
ret.x=ret.y=0;
for (i=1; i<n-1; i++)
if (fabs(t2=xmult(p[i+1],p[0],p[i]))>eps) {
t=barycenter(p[0],p[i],p[i+1]);
ret.x+=t.x*t2;
ret.y+=t.y*t2;
t1+=t2;
}
if (fabs(t1)>eps)
ret.x/=t1,ret.y/=t1;
return ret;
}
//求多边形面积
inline double getarea(point pg[],int n) {
double area=0;
pg[n]=pg[0];
for(int i=0; i<n; i++)
area+=xmult(pg[i].x,pg[i].y,pg[i+1].x,pg[i+1].y);
return fabs(area)/2.0;
}
//解方程ax^2+bx+c=0
int equaltion(double a,double b,double c,double &x1,double &x2) {
double der=b*b-4*a*c;
switch(sign(der)) {
case -1:
return 0;
case 0:
x1=-b/(2*a);
return 1;
case 1:
der=sqrt(der);
x1=(-b-der)/(2*a);
x2=(-b+der)/(2*a);
return 2;
}
}
//线段与圆交点
int line_circle_intersection(line u,circle c,point &p1,point &p2) {
double dis=lenth(u.b-u.a);
point d=(u.b-u.a)/dis;
point E=c.o-u.a;
double a=dmult(E,d);
double a2=a*a;
double e2=dmult(E,E);
double r2=c.r*c.r;
if((r2-e2+a2)<0){
return 0;
} else {
double f=sqrt(r2 - e2 + a2);
double t=a-f;
int cnt=0;
if(t>-eps&&t-dis<eps){//去掉后面变成射线
p1=u.a+(d*t);
cnt++;
}
t=a+f;
if(t>-eps&&t-dis<eps) {
p2=u.a+(d*t);
cnt++;
}
return cnt;
}
}
//给出在任意多边形内部的一个点
point a_point_in_polygon(polygon pg) {
point v,a,b,r;
int i,index;
v=pg.p[0];
index=0;
for(i=1; i<pg.n; i++) {
if(pg.p[i].y<v.y) {
v=pg.p[i];
index=i;
}
}
a=pg.p[(index-1+pg.n)%pg.n];
b=pg.p[(index+1)%pg.n];
point q;
polygon tri;
tri.n=3;
tri.p[0]=a;
tri.p[1]=v;
tri.p[2]=b;
double md=inf;
int in1=index;
bool bin=false;
for(i=0; i<pg.n; i++) {
if(i == index)continue;
if(i == (index-1+pg.n)%pg.n)continue;
if(i == (index+1)%pg.n)continue;
if(!inside_convex2(pg.p[i],tri.p,3))continue;
bin=true;
if(dist(v,pg.p[i])<md) {
q=pg.p[i];
md=dist(v,q);
}
}
if(!bin) {
r.x=(a.x+b.x)/2;
r.y=(a.y+b.y)/2;
return r;
}
r.x=(v.x+q.x)/2;
r.y=(v.y+q.y)/2;
return r;
}
//求在多边形外面的点到凸包的切点
void p_cut_polygon(point p,polygon pg,point &rp,point &lp) {
line ep,en;
bool blp,bln;
rp=pg.p[0];
lp=pg.p[0];
for(int i=1; i<pg.n; i++) {
ep.a=pg.p[(i+pg.n-1)%pg.n];
ep.b=pg.p[i];
en.a=pg.p[i];
en.b=pg.p[(i+1)%pg.n];
blp=xmult(ep.b-ep.a,p-ep.a)>=0;
bln=xmult(en.b-en.a,p-en.a)>=0;
if(!blp&&bln) {
if(xmult(pg.p[i]-p,rp-p)>0)
rp=pg.p[i];
}
if(blp&&!bln) {
if(xmult(lp-p,pg.p[i]-p)>0)
lp=pg.p[i];
}
}
return ;
}
//判断点p在圆c内
bool p_in_circle(point p,circle c) {
return c.r*c.r>dist2(p,c.o);
}
//求矩形第4个点
point rect4th(point a,point b,point c) {
point d;
if(abs(dmult(a-c,b-c))<eps) {
d=a+b-c;
}
if(abs(dmult(a-b,c-b))<eps) {
d=a+c-b;
}
if(abs(dmult(c-a,b-a))<eps) {
d=c+b-a;
}
return d;
}
//判两圆关系
int CircleRelation(circle c1,circle c2){
double d=lenth(c1.o-c2.o);
if( fabs(d-c1.r-c2.r) < eps )
return 2;//外切
if( fabs(d-fabs(c1.r-c2.r)) < eps )
return 4;//内切
if( d > c1.r+c2.r )
return 1;//相离
if( d < fabs(c1.r-c2.r) )
return 5;//内含
if( fabs(c1.r-c2.r) < d && d < c1.r+c2.r)
return 3;//相交
return 0; //error!
}
//判圆与矩形关系,矩形水平
bool Circle_In_Rec(circle c,rectangle r) {
if( r.a.x < c.o.x && c.o.x < r.b.x && r.c.y < c.o.y && c.o.y < r.b.y ) {
line line1(r.a, r.b);
line line2(r.b, r.c);
line line3(r.c, r.d);
line line4(r.d, r.a);
if(c.r<dist_p_to_line(c.o,line1)&&c.r<dist_p_to_line(c.o,line2)&&c.r<dist_p_to_line(c.o,line3)&&c.r<dist_p_to_line(c.o,line4))
return true;
}
return false;
}
//射线关于平面的反射
void reflect(line2 u,line2 v,line2 &l){
double n,m;
double tpb,tpa;
tpb=u.b*v.b+u.a*v.a;
tpa=v.a*u.b-u.a*v.b;
m=(tpb*u.b+tpa*u.a)/(u.b*u.b+u.a*u.a);
n=(tpa*u.b-tpb*u.a)/(u.b*u.b+u.a*u.a);
if(fabs(u.a*v.b-v.a*u.b)<1e-20) {
l.a=v.a;
l.b=v.b;
l.c=v.c;
return;
}
double xx,yy; //(xx,yy)是入射线与镜面的交点。
xx=(u.b*v.c-v.b*u.c)/(u.a*v.b-v.a*u.b);
yy=(v.a*u.c-u.a*v.c)/(u.a*v.b-v.a*u.b);
l.a=n;
l.b=-m;
l.c=m*yy-xx*n;
}
//两圆交点（预判断不相交情况）
void c2point(circle c1,circle c2,point &rp1,point &rp2) {
double a,b,r;
a=c2.o.x-c1.o.x;
b=c2.o.y-c1.o.y;
r=(a*a+b*b+c1.r*c1.r-c2.r*c2.r)/2;
if(a==0&&b!=0) {
rp1.y=rp2.y=r/b;
rp1.x=sqrt(c1.r*c1.r-rp1.y*rp1.y);
rp2.x=-rp1.x;
} else if(a!=0&&b==0) {
rp1.x=rp2.x=r/a;
rp1.y=sqrt(c1.r*c1.r-rp1.x*rp2.x);
rp2.y=-rp1.y;
} else if(a!=0&&b!=0) {
double delta;
delta=b*b*r*r-(a*a+b*b)*(r*r-c1.r*c1.r*a*a);
rp1.y=(b*r+sqrt(delta))/(a*a+b*b);
rp2.y=(b*r-sqrt(delta))/(a*a+b*b);
rp1.x=(r-b*rp1.y)/a;
rp2.x=(r-b*rp2.y)/a;
}
rp1=rp1+c1.o;
rp2=rp2+c1.o;
}
//圆外一点引圆的切线
void cutpoint(circle c,point sp,point &rp1,point &rp2) {
circle c2;
c2.o=(c.o+sp)/2.0;
c2.r=lenth(c2.o-sp);
c2point(c,c2,rp1,rp2);
}
//圆c1上，与c2的外切点
void c2cuto(circle c1, circle c2, point &p1, point &p2) {
double d = dist(c1.o, c2.o), dr = c1.r - c2.r;
double b = acos(dr / d);
double a = angle(c2.o-c1.o);
double a1 = a - b, a2 = a + b;
p1=point(cos(a1) * c1.r, sin(a1) * c1.r) + c1.o;
p2=point(cos(a2) * c1.r, sin(a2) * c1.r) + c1.o;
}
//圆c1上，与c2的内切点
void c2cuti(circle c1,circle c2,point &p1,point &p2){
point dr=c2.o-c1.o;
dr=dr/lenth(dr);
point a=c1.o-(dr*c1.r),b=c1.o+(dr*c1.r);
point c=c2.o-(dr*c2.r),d=c2.o+(dr*c2.r);
circle E((a+c)/2.0,lenth(c-a)/2.0),F((b+d)/2.0,lenth(d-b)/2.0);
point q1,q2;
c2point(E,F,q1,q2);
point L=line_intersection2(line(c1.o,c2.o),line(q1,q2));
circle c3((c1.o+L)/2.0,lenth(L-c1.o)/2.0);
c2point(c1,c3,p1,p2);
}
时间： 2024-08-30 05:53:10