}
/*（sp-op）*（ep-op）的叉积
r=multiply（sp,ep,op）,得到（sp-op）*（ep-op）的叉积
r>0:
sp在矢量opep的顺时针方向；
r=0：
opspep三点共线；
r<0:
sp在矢量opep的逆时针方向*/
doublemultiply（d_POINTsp,d_POINTep,d_POINTop）
{
return（（sp.x-op.x）*（ep.y-op.y）-（ep.x-op.x）*（sp.y-op.y））;
}
doubleamultiply（d_POINTsp,d_POINTep,d_POINTop）
{
returnfabs（（sp.x-op.x）*（ep.y-op.y）-（ep.x-op.x）*（sp.y-op.y））;
}
/*矢量（p1-op）和（p2-op）的点积
r=dotmultiply（p1,p2,op）,得到矢量（p1-op）和（p2-op）的点积如果两个矢量都非零矢量
r<0:
两矢量夹角为锐角；
r=0：
两矢量夹角为直角；
r>0:
两矢量夹角为钝角*/
doubledotmultiply（d_POINTp1,d_POINTp2,d_POINTp0）
{
return（（p1.x-p0.x）*（p2.x-p0.x）+（p1.y-p0.y）*（p2.y-p0.y））;
}
/*判断点p是否在线段l上
条件：
（p在线段l所在的直线上）&&（点p在以线段l为对角线的矩形内）*/
boolonline（LINESEGl,d_POINTp）
{
return（（multiply（l.e,p,l.s）==0）
&&（（（p.x-l.s.x）*（p.x-l.e.x）<=0）&&（（p.y-l.s.y）*（p.y-l.e.y）<=0）））;
}
//返回点p以点o为圆心逆时针旋转alpha（单位：
弧度）后所在的位置
d_POINTrotate（d_POINTo,doublealpha,d_POINTp）
{
d_POINTtp;
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;
returntp;
}
/*返回顶角在o点，起始边为os，终止边为oe的夹角（单位：
弧度）
角度小于pi，返回正值
角度大于pi，返回负值
可以用于求线段之间的夹角*/
doubleangle（d_POINTo,d_POINTs,d_POINTe）
{
doublecosfi,fi,norm;
doubledsx=s.x-o.x;
doubledsy=s.y-o.y;
doubledex=e.x-o.x;
doubledey=e.y-o.y;
cosfi=dsx*dex+dsy*dey;
norm=（dsx*dsx+dey*dey）*（dex*dex+dey*dey）;
cosfi/=sqrt（norm）;
if（cosfi>=1.0）return0;
if（cosfi<=-1.0）return-3.1415926;
fi=acos（cosfi）;
if（dsx*dey-dsy*dex>0）returnfi;//说明矢量os在矢量oe的顺时针方向
return-fi;
}
/*判断点C在线段AB所在的直线l上垂足P的与线段AB的关系
本函数是根据下面的公式写的，P是点C到线段AB所在直线的垂足
ACdotAB
r=----------------------
||AB||^2
（Cx-Ax）（Bx-Ax）+（Cy-Ay）（By-Ay）
=----------------------------------------------------
L^2
rhasthefollowingmeaning:
r=0P=A
r=1P=B
r<0PisonthebackwardextensionofAB
r>1PisontheforwardextensionofAB
0*/
doublerelation（d_POINTc,LINESEGl）
{
LINESEGtl;
tl.s=l.s;
tl.e=c;
returndotmultiply（tl.e,l.e,l.s）/（dist（l.s,l.e）*dist（l.s,l.e））;
}
//求点C到线段AB所在直线的垂足P
d_POINTperpendicular（d_POINTp,LINESEGl）
{
doubler=relation（p,l）;
d_POINTtp;
tp.x=l.s.x+r*（l.e.x-l.s.x）;
tp.y=l.s.y+r*（l.e.y-l.s.y）;
returntp;
}
/*求点p到线段l的最短距离
返回线段上距该点最近的点np注意：
np是线段l上到点p最近的点，不一定是垂足*/
doubleptolinesegdist（d_POINTp,LINESEGl,d_POINT&np）
{
doubler=relation（p,l）;
if（r<0）
{
np=l.s;
returndist（p,l.s）;
}
if（r>1）
{
np=l.e;
returndist（p,l.e）;
}
np=perpendicular（p,l）;
returndist（p,np）;
}
//求点p到线段l所在直线的距离
//请注意本函数与上个函数的区别
doubleptoldist（d_POINTp,LINESEGl）
{
returnabs（multiply（p,l.e,l.s））/dist（l.s,l.e）;
}
/*计算点到折线集的最近距离,并返回最近点.
注意：
调用的是ptolineseg（）函数*/
doubleptopointset（intvcount,d_POINTpointset[],d_POINTp,d_POINT&q）
{
inti;
doublecd=double（INF）,td;
LINESEGl;
d_POINTtq,cq;
for（i=0;i{
l.s=pointset[i];
l.e=pointset[i+1];
td=ptolinesegdist（p,l,tq）;
if（td{
cd=td;
cq=tq;
}
}
q=cq;
returncd;
}
/*判断圆是否在多边形内*/
boolCircleInsidePolygon（intvcount,d_POINTcenter,doubleradius,d_POINTpolygon[]）
{
d_POINTq;
doubled;
q.x=0;
q.y=0;
d=ptopointset（vcount,polygon,center,q）;
if（delsereturnfalse;
}
/*返回两个矢量l1和l2的夹角的余弦（-1~1）
注意：
如果想从余弦求夹角的话，注意反余弦函数的值域是从0到pi*/
doublecosine（LINESEGl1,LINESEGl2）
{
doubled=（（l1.e.x-l1.s.x）*（l2.e.x-l2.s.x）+（l1.e.y-l1.s.y）*（l2.e.y-l2.s.y））;
return（d/（dist（l1.e,l1.s）*dist（l2.e,l2.s）））;
}
//返回线段l1与l2之间的夹角
//单位：
弧度范围（-pi，pi）
doublelsangle（LINESEGl1,LINESEGl2）
{
d_POINTo,s,e;
o.x=o.y=0;
s.x=l1.e.x-l1.s.x;
s.y=l1.e.y-l1.s.y;
e.x=l2.e.x-l2.s.x;
e.y=l2.e.y-l2.s.y;
returnangle（o,s,e）;
}
//判断线段u和v相交（包括相交在端点处）
boolintersect（LINESEGu,LINESEGv）
{
return（（max（u.s.x,u.e.x）>=min（v.s.x,v.e.x））&&//排斥实验
（max（v.s.x,v.e.x）>=min（u.s.x,u.e.x））&&
（max（u.s.y,u.e.y）>=min（v.s.y,v.e.y））&&
（max（v.s.y,v.e.y）>=min（u.s.y,u.e.y））&&
（multiply（v.s,u.e,u.s）*multiply（u.e,v.e,u.s）>=0）&&//跨立实验
（multiply（u.s,v.e,v.s）*multiply（v.e,u.e,v.s）>=0））;
}
//判断线段u和v相交（不包括双方的端点）
boolintersect_A（LINESEGu,LINESEGv）
{
return（（intersect（u,v））&&
（!
online（u,v.s））&&
（!
online（u,v.e））&&
（!
online（v,u.e））&&
（!
online（v,u.s）））;
}
//判断线段v所在直线与线段u相交
//方法：
判断线段u是否跨立线段v
boolintersect_l（LINESEGu,LINESEGv）
{
returnmultiply（u.s,v.e,v.s）*multiply（v.e,u.e,v.s）>=0;
}
//已知两点坐标p1和p2，求过这两点的直线解析方程：
a*x+b*y+c=0（a>=0）
LINEmakeLine（d_POINTp1,d_POINTp2）
{
LINEtl;
intsign=1;
tl.a=p2.y-p1.y;
if（tl.a<0）
{
sign=-1;
tl.a=sign*tl.a;
}
tl.b=sign*（p1.x-p2.x）;
tl.c=sign*（p1.y*p2.x-p1.x*p2.y）;
returntl;
}
//根据直线解析方程返回直线的斜率k,水平线返回0,竖直线返回1e200
doubleslope（LINEl）
{
if（abs（l.a）<1e-20）return0;
if（abs（l.b）<1e-20）returnINF;
return-（l.a/l.b）;
}
//返回直线的倾斜角alpha（0-pi）
//注意：
atan（）返回的是-PI/2~PI/2
doublealpha（LINEl）
{
if（abs（l.a）if（abs（l.b）doublek=slope（l）;
if（k>0）
returnatan（k）;
else
returnPI+atan（k）;
}
//求点p关于直线l的对称点
d_POINTsymmetry（LINEl,d_POINTp）
{
d_POINTtp;
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）;
returntp;
}
//如果两条直线l1（a1*x+b1*y+c1=0）,l2（a2*x+b2*y+c2=0）相交，返回true，且返回交点p
boollineInterSect（LINEl1,LINEl2,d_POINT&p）//是L1，L2
{
doubled=l1.a*l2.b-l2.a*l1.b;
if（abs（d）returnfalse;
p.x=（l2.c*l1.b-l1.c*l2.b）/d;
p.y=（l2.a*l1.c-l1.a*l2.c）/d;
returntrue;
}
//如果线段l1和l2相交，返回true且交点由（inter）返回，否则返回false
boolintersection（LINESEGl1,LINESEGl2,d_POINT&inter）
{
LINEll1,ll2;
ll1=makeLine（l1.s,l1.e）;
ll2=makeLine（l2.s,l2.e）;
if（lineInterSect（ll1,ll2,inter））
returnonline（l1,inter）;
else
returnfalse;
}
//如果无特别说明，输入多边形顶点要求按逆时针排列
//返回多边形面积（signed）；
//输入顶点按逆时针排列时，返回正值；否则返回负值
doublearea_of_polygon（intvcount,d_POINTpolygon[]）
{
inti;
doubles;
if（vcount<3）
return0;
s=polygon[0].y*（polygon[vcount-1].x-polygon[1].x）;
for（i=1;is+=polygon[i].y*（polygon[（i-1）].x-polygon[（i+1）%vcount].x）;
returns/2;
}
//判断顶点是否按逆时针排列
//如果输入顶点按逆时针排列，返回true
boolisconterclock（intvcount,d_POINTpolygon[]）
{
returnarea_of_polygon（vcount,polygon）>0;
}
/*射线法判断点q与多边形polygon的位置关系
要求polygon为简单多边形，顶点时针排列
如果点在多边形内：
返回0
如果点在多边形边上：
返回1
如果点在多边形外：
返回2*/
intinsidepolygon（intvcount,d_POINTPolygon[],d_POINTq）
{
intc=0,i;
LINESEGl1,l2;
l1.s=q;l1.e=q;l1.e.x=double（INF）;
//n=vcount;
for（i=0;i{
l2.s=Polygon[i];
l2.e=Polygon[（i+1）%vcount];
//doubleee=Polygon[（i+2）%vcount].x;
//doubless=Polygon[（i+3）%vcount].y;
if（online（l2,q））
return1;
if（intersect_A（l1,l2））
c++;//相交且不在端点
if（online（l1,l2.e）&&!
online（l1,l2.s）&&l2.e.y>l2.e.y）
c++;//l2的一个端点在l1上且该端点是两端点中纵坐标较大的那个
if（!
online（l1,l2.e）&&online（l1,l2.s）&&l2.e.yc++;//忽略平行边
}
if（c%2==1）
return0;
else
return2;
}
//判断点q在凸多边形polygon内
//点q是凸多边形polygon内[包括边上]时，返回true
//注意：
多边形polygon一定要是凸多边形
boolInsideConvexPolygon（intvcount,d_POINTpolygon[],d_POINTq）
{
d_POINTp;
LINESEGl;
inti;
p.x=0;p.y=0;
for（i=0;i多边形顶点平均值
{
p.x+=polygon[i].x;
p.y+=polygon[i].y;
}
p.x/=vcount;
p.y/=vcount;
for（i=0;i{
l.s=polygon[i];
l.e=polygon[（i+1）%vcount];
if（multiply（p,l.e,l.s）*multiply（q,l.e,l.s）<0）
/*点p和点q在边l的两侧，说明点q肯定在多边形外*/
returnfalse;
}
returntrue;
}
/*寻找凸包的graham扫描法
PointSet为输入的点集；
ch为输出的凸包上的点集，按照逆时针方向排列;
n为PointSet中的点的数目
len为输出的凸包上的点的个数*/
voidGraham_scan（d_POINTPointSet[],d_POINTch[],intn,int&len）
{
inti,j,k=0,top=2;
d_POINTtmp;
//选取PointSet中y坐标最小的点PointSet[k]，如果这样的点有多个，则取最左边的一个
for（i=1;iif（PointSet[i].y&&（PointSet[i].xk=i;
tmp=PointSet[0];
PointSet[0]=PointSet[k];
PointSet[k]=tmp;//现在PointSet中y坐标最小的点在PointSet[0]
for（i=1;i的按照距离PointSet[0]从近到远进行排序*/
{
k=i;
for（j=i+1;jif（multiply（PointSet[j],PointSet[k],PointSet[0]）>0||//极角更小
（multiply（PointSet[j],PointSet[k],PointSet[0]）==0）&&/*极角相等，距离更短*/dist（PointSet[0],PointSet[j]）k=j;
tmp=PointSet[i];
PointSet[i]=PointSet[k];
PointSet[k]=tmp;
}
ch[0]=PointSet[0];
ch[1]=PointSet[1];
ch[2]=PointSet[2];
for（i=3;i{
while（multiply（PointSet[i],ch[top],ch[top-1]）>=0）top--;
ch[++top]=PointSet[i];
}
len=top+1;
}
//求凸多边形的重心,要求输入多边形按逆时针排序
d_POINTgravitycenter（intvcount,d_POINTpolygon[]）
{
d_POINTtp;
doublex,y,s,x0,y0,cs,k;
x=0;y=0;s=0;
for（inti=1;i{
x0=（polygon[0].x+polygon[i].x+polygon[i+1].x）/3;
y0=（polygon[0].y+polygon[i].y+polygon[i+1].y）/3;//求当前三角形的重心
cs=multiply（polygon[i],polygon[i+1],polygon[0]）/2;
//三角形面积可以直接利用该公式求解
if（abs（s）<1e-20）
{
x=x0;y=y0;s+=cs;continue;
}
k=cs/s;//求面积比例
x=（x+k*x0）/（1+k）;
y=（y+k*y0）/（1+k）;
s+=cs;
}
tp.x=x;
tp.y=y;
returntp;
}
/*所谓凸多边形的直径，即凸多边形任两个顶点的最大距离。
下面的算法
仅耗时O（n），是一个优秀的算法。
输入必须是一个凸多边形，且顶点
必须按顺序（顺时针、逆时针均可）依次输入。
若输入不是凸多边形
而是一般点集，则要先求其凸包。
就是先求出所有跖对，然后求出每
个跖对的距离，取最大者。
点数要多于5个*/
voidDiameter（d_POINTch[],intn,double&dia）
{
intznum=0,i,j,k=1;
intzd[MAXV][2];
doubletmp
升级会员 