浙江大学ACM模板

ZhejiangUniversity

ICPCTeam

RoutineLibrary

byWishingBone(Dec.2002)

LastUpdate(Nov.2004)byRivena

1

1、几何......................................................................5

11注意5

1.3多边形..............................................................7

1.4多边形切割.........................................................10

15浮点函数J]

।&二角形[8

19三维几何...........................................................20

1.10凸包..............................................................27

1.11网格...............................................................29

1.12圆................................................................29

1.13整数函数..........................................................31

2、组合34

21组合公式34

2.2排列组合生成.......................................................34

2.3生成gray码........................................................36

2.4置换(polya).........................................................................................................................36

2.5字典序全排列.......................................................37

2.6字典序组合.........................................................37

3、结构.....................................................................38

3.1并查集.............................................................38

3.2堆.................................................................39

3.3线段树.............................................................40

3.4子段和.............................................................45

3.5子阵和.............................................................45

4、数论46

4.1阶乘最后非。位.....................................................46

4.2模线性方程组.......................................................47

4.3素数...............................................................48

5、数值计算50

5.2多项式求根(牛顿法).................................................52

5.3周期性方程(追赶法).................................................53

6、图论一NP搜索...........................................................54

6.1最大团.............................................................54

6.2最大团(n<64)(faster).........................................................................................................55

7、图论一连通性57

7•1(dfsqB)•••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••57

7.3无向图的块(bfs邻接阵)..............................................59

7.4无向图连通分支(dfs/bfs邻接阵).......................................60

2

7.5有向图强连通分支（dfs/bfs邻接阵）.....................................61

7.6有向图最小点基（邻接阵）.............................................62

8、图论一匹配..............................................................63

8.1二分图最大匹配（hungary邻接表）......................................63

8.2二分图最大匹配（hungary邻接阵）......................................64

8.3二分图最大匹配（hungary正向表）......................................64

8.4二分图最佳匹配（kuhnjnunkras邻接阵）................................65

8.5一般图匹配（邻接表）.................................................66

8.6一般图匹配（邻接阵）.................................................67

8.7一般图匹配（正向表）.................................................67

9、图论一网络流............................................................68

9.1最大流（邻接阵）.....................................................68

9.2上下界最大流（邻接阵）...............................................69

93上下界最小流（邻接阵）................................................70

9.4最大流无流量（邻接阵）...............................................71

9.5最小费用最大流（邻接阵）.............................................71

10、图论一应用.............................................................72

10.1欧拉回路（邻接阵）..................................................72

10.2树的前序表转化....................................................73

10.3树的优化算法......................................................74

10.4拓扑排序（邻接阵）..................................................75

10.5最佳边割集........................................................76

10.6最佳点割集........................................................77

10.7最小边割集........................................................78

10.8最小点割集........................................................79

10.9最小路径覆盖......................................................81

Ik图论一支撑树...........................................................81

11.1最小生成树（kruskal邻接表）..........................................81

11.2最小生成树（kruskal正向表）..........................................83

11.3最小生成树（prim+binary_heap邻接表）................................84

11.4最小生成树（prim+binary_heap正向表）................................85

11.5最小生成树（prim+mapped_heap邻接表）...............................86

11.6最小生成树（prim+mapped_heap正向表）..............................88

11.7最小生成树（prim邻接阵）............................................89

11.8最小树形图（邻接阵）................................................89

12、图论一最短路径.........................................................91

12.1最短路径（单源bellman_ford邻接阵）..................................91

12.2最短路径（单源dijkstra+bfs邻接表）...................................91

12.3最短路径（单源dijkstra+bfs正向表）...................................92

12.4最短路径（单源dijkstra+binary_heap邻接表）...........................93

12.5最短路径（单源dijkstra+binary_heap正向表）...........................94

12.6最短路径（单源dijkstra+m叩ped_heap邻接表）.........................95

12.7最短路径（单源dijkstra+mapped_heap正向表）.........................96

12.8最短路径（单源dijkstra邻接阵）.......................................97

3

12.9最短路径（多源floyd_warshall邻接阵）.................................98

13、应用...................................................................98

13.1Joseph问题........................................................98

13.2N皇后构造解......................................................99

13.3布尔母函数.......................................................100

13.4第k元素.........................................................100

13.5幻方构造.........................................................101

13.6模式匹配（kmp）........................................................................................................102

13.7逆序对数.........................................................103

13.8字符串最小表示...................................................103

13.9最长公共单调子序列..............................................104

13.10最长子序列......................................................105

13.11最大子串匹配....................................................106

13.12最大子段和......................................................107

13.13最大子阵和......................................................107

14、其它...................................................................108

14.1大数（只能处理正数）...............................................108

14.2分数.............................................................114

14.3矩阵.............................................................116

14.4线性方程组.......................................................118

14.5线性相关.........................................................120

14.6日期.............................................................120

4

1、几何

1.1注意

1.注意舍入方式(0.5的舍入方向);防止输出-0.

2.几何题注意多测试不对称数据.

3.整数几何注意xmult和dmult是否会出界；

符点几何注意eps的使用.

4.避免使用斜率;注意除数是否会为0.

5.公式一定要化简后再代入.

6.判断同一个2*PI域内两角度差应该是

abs(al-a2)<betallabs(a1-a2)>pi+pi-beta;

相等应该是

abs(al-a2)<epsllabs(al-a2)>pi+pi-eps;

7.需要的话尽量使用atan2,注意:atan2(0,0)=0,

atan2(1,0)=pi/2,atan2(-1,0)=-pi/2,atan2(0,l)=0,atan2(0,-1)=pi.

8.crossproduct=lul*lvl*sin(a)

dotproduct=lul*lvl*cos(a)

9.(Pl-P0)x(P2-P0)结果的意义：

正:＜PO,P1＞在＜P0,P2＞顺时针(0,pi)内

负：＜P(),P1＞在＜P0,P2＞逆时针(0,pi)内

0:＜P0,Pl＞,＜P0,P2＞共线,夹角为0或pi

10.误差限缺省使用le-8!

1.2几何公式

三角形：

1.半周长P=(a+b+c)/2

2.面积S=aHa/2=absin(C)/2=sqrt(P(P-a)(P-b)(P-c))

3.中线Ma=sqrt(2(bA2+cA2)-aA2)/2=sqrt(bA2+cA2+2bccos(A))/2

4.角平分线Ta=sqrt(bc((b+c)A2-aA2))/(b+c)=2bccos(A/2)/(b+c)

5.高线Ha=bsin(C)=csin(B)=sqrt(bA2-((aA2+bA2-cA2)/(2a))A2)

5

6.内切圆半径r=S/P=asin(B/2)sin(C/2)/sin((B+C)/2)

=4Rsin(A/2)sin(B/2)sin(C/2)=sqrt((P-a)(P-b)(P-c)/P)

=Ptan(A/2)tan(B/2)tan(C/2)

7,外接圆半径R=abc/(4S)=a/(2sin(A))=b/(2sin(B))=c/(2sin(C))

四边形：

D1.D2为对角线,M对角线中点连线,A为对角线夹角

1.aA2+bA2+cA2+dA2=D1A2+D2A2+4MA2

2.S=DlD2sin(A)/2

(以下对圆的内接四边形)

3.ac+bd=D1D2

4.S=sqrt((P-a)(P-b)(P-c)(P-d)),P为半周长

正n边形：

R为外接圆半径,r为内切圆半径

1.中心角A=2PI/n

2.内角C=(n-2)PI/n

3.边长a=2sqrt(RA2-rA2)=2Rsin(A/2)=2rtan(A/2)

4.面积S=nar/2=nrA2tan(A/2)=nRA2sin(A)/2=naA2/(4tan(A/2))

圆：

1.弧长l=rA

2.弦长a=2sqrt(2hr-hA2)=2rsin(A/2)

3.弓形高h=r-sqrt(rA2-aA2/4)=r(1-cos(A/2))=atan(A/4)/2

4.扇形面积Sl=rl/2=rA2A/2

5.弓形面积S2=(rl-a(r-h))/2=rA2(A-sin(A))/2

棱柱：

1.体积V=Ah,A为底面积,h为高

2.侧面积S=lp,l为棱长,p为直截面周长

3.全面积T=S+2A

棱锥：

1.体积V=Ah/3,A为底面积,h为高

(以下对正棱锥)

2.侧面积S=lp/2,1为斜高,p为底面周长

3.全面积T=S+A

棱台：

1.体积V=(Al+A2+sqrt(AlA2))h/3,Al.A2为上下底面积,h为高

(以下为正棱台)

2.侧面积S=(pl+p2)l/2,pl.p2为上下底面周长,1为斜高

3.全面积T=S+A1+A2

6

圆柱：

1.侧面积S=2PIrh

2.全面积T=2PIr(h+r)

3.体积V=PIrA2h

圆锥：

1.母线l=sqrt(hA2+rA2)

2.侧面积S=PIrl

3.全面积T=PIr(l+r)

4.体积V=PIrA2h/3

圆台：

1.母线l=sqrt(hA2+(rl-r2)A2)

2.侧面积S=PI(rl+⑵1

3.全面积T=PIrl(l+rl)+PIr2(l+r2)

4.体积V=PI(rlA2+r2A2+rlr2)h/3

球：

1.全面积T=4P1rA2

2.体积V=4PIM3/3

球台：

1.侧面积S=2PIrh

2.全面积T=PI(2rh+rlA2+r2A2)

3.体积V=PIh(3(r1A2+r2A2)+hA2)/6

球扇形：

1.全面积T=PIr(2h+rO),h为球冠高,rO为球冠底面半径

2.体积V=2P1rA2h/3

1.3多边形

#include<stdlib.h>

#include<math.h>

#defineMAXN1000

#defineoffset10000

#defineepsle-8

#definezero(x)(((x)>0?(x):-(x))<eps)

#define_sign(x)((x)>eps?1:((x)<-eps?2:0))

structpoint{doublex,y;};

structline{pointa,b;};

doublexmult(pointpl,pointp2,pointp0){

return(pl.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(pl.y-p0.y);

7

〃判定凸多边形,顶点按顺时针或逆时针给出，允许相邻边共线

intis_convex(intn,point*p){

inti,s[3]={1,1,1);

for(i=0;i<n&&s[l]ls[2];i++)

sLsign(xmult(p[(i+l)%n],p[(i+2)%n],p[i]))]=0;

returns[l]ls[2];

〃判定凸多边形,顶点按顺时针或逆时针给出，不允许相邻边共线

intis_convex_v2(intn,point*p){

inti,s[3]={1,1,1);

for(i=0;i<n&&s[0]&&s[l]ls[2];i++)

s[_sign(xmult(p[(i+l)%nJ,p[(i+2)%n],p[i]))]=0;

returns[0]&&s[l]ls[2];

)

〃判点在凸多边形内或多边形边上，顶点按顺时针或逆时针给出

intinside_convex(pointq,intn,point*p){

inti,s[3]={1,1,1);

for(i=0;i<n&&s[l]ls[2];i++)

sLsign(xmult(p[(i+l)%n],q,p[i]))]=O;

returns[l]ls[2];

)

〃判点在凸多边形内，顶点按顺时针或逆时针给出，在多边形边上返回0

intinside_convex_v2(pointq,intn,point*p){

inti,s[3]={1,1,1);

for(i=0;i<n&&s[0]&&s[l]ls[2];i++)

sLsign(xmult(p[(i+l)%n],q,p[i]))]=O;

returns[0]&&s[l]ls[2];

〃判点在任意多边形内,顶点按顺时针或逆时针给出

//on_edge表示点在多边形边上时的返回值,offset为多边形坐标上限

intinside_polygon(pointq,intn,point*p,inton_edge=l){

pointq2;

inti=0,count;

while(i<n)

for(count=i=0,q2.x=rand()+offset,q2.y=rand()+offset;i<n;i4-4-)

if

(zero(xmult(q,p[i],p[(i4-l)%n]))&&(p[i].x-q.x)*(p[(i+l)%n].x-q.x)<eps&&(p[i].y-q.y)*(p[(i+l)%

n].y-q.y)<eps)

8

returnon_edge;

elseif(zero(xmult(q,q2,p[i])))

break;

elseif

(xmult(q,p[i],q2)*xmult(q,p[(i+l)%n],q2)<-eps&&xmult(p[i],q,p[(i+l)%n])*xmult(p[i],q2,p[(i+l)

%n])<-eps)

count++;

returncount&l;

)

inlineintopposite_side(pointpl,pointp2,point11,point12){

returnxmult(ll,p1,12)*xmult(ll,p2,12)<-eps;

)

inlineintdot_online_in(pointp,point11,point12){

returnzero(xmult(p,ll,12))&&(ll.x-p.x)*(12.x-p.x)<eps&&(ll.y-p.y)*(12.y-p.y)<eps;

)

〃判线段在任意多边形内,顶点按顺时针或逆时针给出，与边界相交返回1

intinside_polygon(point11,point12,intn,point*p){

pointt[MAXN],tt;

inti,j,k=O;

if(!inside_polygon(l1,n,p)II!inside_polygon(12,n,p))

return0;

fbr(i=0;i<n;i++)

if(opposite_side(l1,12,p[i],p[(i+l)%n])&&opposite_side(p[i],p[(i+l)%n],11,12))

return0;

elseif(dot_online_in(l1,p[i],p[(i+1)%n]))

t[k++]=ll;

elseif(dot_online_in(12,p[i],p[(i+l)%n]))

t[k++]=12;

elseif(dot_online_in(p|i],ll,12))

t[k++]=p[i];

for(i=0;i<k;i++)

for(j=i+l;j<k;j++){

tt.x=(t[i].x+t[j].x)/2;

tt.y=(t[i].y+t[j].y)/2;

if(!inside_polygon(tt,n,p))

return0;

)

return1;

)

pointintersection(lineu,linev){

9

pointret=u.a;

doublet=((u.a.x-v.a.x)*(v.a.y-v.b.y)-(u.a.y-v.a.y)*(v.a.x-v.b.x))

/((u.a.x-u.b.x)*(v.a.y-v.b.y)-(u.a.y-u.b.y)*(v.a.x-v.b.x));

ret.x+=(u.b.x-u.a.x)*t;

ret.y+=(u.b.y-u.a.y)*t;

returnret;

)

pointbarycenter(pointa,pointb,pointc){

lineu,v;

u.a.x=(a.x+b.x)/2;

u.a.y=(a.y+b.y)/2;

irb=c;

v.a.x=(a.x+c.x)/2;

v.a.y=(a.y+c.y)/2;

v.b=b;

returnintersection(u,v);

)

〃多边形重心

pointbarycenter(intn,point*p){

pointret,t;

doubletl=O,t2;

inti;

ret.x=ret.y=O;

fbr(i=l;i<n-l;i++)

if(fabs(t2=xmult(p[0],p[i],p[i+lj))>eps){

t=barycenter(p[0],p[iJ,p[i4-1]);

ret.x+=t.x*t2;

ret.y+=t.y*t2;

tl+=t2;

)

if(fabs(tl)>eps)

ret.x/=tl,ret.y/=tl;

returnret;

)

1.4多边形切割

〃多边形切割

〃可用于半平面交

#defineMAXN100

#defineepsle-8

#definezero(x)(((x)>0?(x):-(x))<eps)

10

structpoint{doublex,y;};

doublexmult(pointpl,pointp2,pointp0){

return(p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(pl.y-pO.y);

)

intsame_side(pointpl,pointp2,point11,point12){

returnxmult(ll,pl,12)*xmult(ll,p2,12)>eps;

)

pointintersection(pointul,pointu2,pointvl,pointv2){

pointret=u1;

doublet=((ul.x-vl.x)*(vl.y-v2.y)-(u1.y-vl.y)*(vl.x-v2.x))

/((u1.x-u2.x)*(v1.y-v2.y)-(ul.y-u2.y)*(vl.x-v2.x));

ret.x+=(u2.x-ul.x)*t;

ret.y+=(u2.y-ul.y)*t;

returnret;

)

〃将多边形沿11,12确定的直线切割在side侧切割，保证11,12,side不共线

voidpolygon_cut(int&n,point*p,point11,point12,pointside){

pointpp[100];

intm=0,i;

fbr(i=0;i<n;i++){

if(same_side(p[i],side,l1,12))

pp[m++]=p[i];

if

(!same_side(p[i],p[(i+l)%n],ll,12)&&!(zero(xmult(p[i],l1,12))&&zero(xmult(p[(i+l)%n],l1,12))))

pp[m++]=intersection(p[i],p[(i4-l)%n],ll,12);

)

for(n=i=0;i<m;i++)

if(!ill!zero(pp[i].x-pp[i-l].x)ll!zero(pp[i].y-pp[i-l].y))

p[n++]=pp[i];

if(zero(p[n-1].x-p[0].x)&&zero(p[n-1].y-p[0].y))

n-;

if(n<3)

n=0;

)

1.5浮点函数

〃浮点几何函数库

#include<math.h>

#defineepsle-8

11

#definezero(x)(((x)>0?(x):-(x))<eps)

structpoint{doublex,y;};

structline{pointa,b;};

〃计算crossproduct(Pl-P0)x(P2-P0)

doublexmult(pointpl,pointp2,pointp0){

return(pl.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(pl.y-p0.y);

)

doublexmult(doublex1,doubleyl,doublex2,doubley2,doublexO,doubley0){

return(x1-x0)*(y2-y0)-(x2-x0)*(y1-yO);

)

〃计算dotproduct(Pl-P0).(P2-P0)

doubledmult(pointpinpointp2,pointp0){

return(pl.x-p0.x)*(p2.x-p0.x)+(pl.y-p0.y)*(p2.y-p0.y);

)

doubledmult(doublexl,doubleyl,doublex2,doubley2,doublexO,doubley0){

return(x1・x0)*(x2・x0)+(y1-y0)*(y2-y0);

)

〃两点距离

doubledistance(pointpinpointp2){

returnsqrt((pl.x-p2.x)*(pl.x-p2.x)+(pl.y-p2.y)*(pl.y-p2.y));

)

doubledistance(doublexl,doubleyl,doublex2,doubley2){

returnsqrt((xl-x2)*(xl-x2)+(yl-y2)*(yl-y2));

)

〃判三点共线

intdots_inline(pointpl,pointp2,pointp3){

returnzero(xmult(p1,p2,p3));

)

intdots_inline(doublex1,doubleyl,doublex2,doubley2,doublex3,doubley3){

returnzero(xmult(xl,yl,x2,y2,x3,y3));

}

〃判点是否在线段上，包括端点

intdot_online_in(pointpjinel){

returnzero(xmult(p,l.a,l.b))&&(l.a.x-p.x)*(l.b.x-p.x)<eps&&(l.a.y-p.y)*(l.b.y-p.y)<eps;

)

intdot_online_in(pointp,point1hpoint12){

returnzero(xmult(p,l1,I2))&&(1l.x-p.x)*(12.x-p.x)<eps&&(l1.y-p.y)*(12.y-p.y)<eps;

)

intdot_online_in(doublex,doubley,doublexI,doubleyl,doublex2,doubley2){

12

returnzero(xmult(x,y,xl,yl,x2,y2))&&(xl-x)*(x2-x)<eps&&(yl-y)*(y2-y)<eps;

〃判点是否在线段上，不包括端点

intdot_online_ex(pointp,line1){

return

dot_online_in(p,l)&&(!zero(p.x・l.a.x)ll!zero(p.y・l.a.y))&&(!zero(p.x-Lb.x)ll!zero(p.y・l.b.y));

)

intdot_online_ex(pointp,point11,point12){

return

dot_online_in(p,ll,12)&&(!zero(p.x-l1.x)II!zero(p.y-l1.y))&&(!zero(p.x-12.x)ll!zero(p.y-12.y));

)

intdot_online_ex(doublex,doubley,doublexl,doubleyl,doublex2,doubley2){

return

dot_online_in(x,y,xl,yl,x2,y2)&&(!zero(x-xl)ll!zero(y-yl))&&(!zero(x-x2)ll!zero(y-y2));

)

〃判两点在线段同侧，点在线段上返回0

intsame_side(pointpl,pointp2,line1){

returnxmult(l.a,p1,Lb)：icxmult(l.a,p2,l.b)>eps;

)

intsame_side(pointpl,pointp2,point11,point12){

returnxmult(ll,p1,12)*xmult(l1,p2,12)>eps;

)

〃判两点在线段异侧，点在线段上返回0

intopposite_side(pointphpointp2,line1){

returnxmult(l.a,p1J.b)：icxmult(l.a,p2,l.b)<-eps;

)

intopposite_side(pointpinpointp2,point11,point12){

returnxmult(ll,pl,12)*xmult(ll,p2,12)<-eps;

)

〃判两直线平行

intparallel(lineu,linev){

returnzero((u.a.x-u.b.x)*(v.a.y-v.b.y)-(v.a.x-v.b.x)*(u.a.y-u.b.y));

)

intparallel(pointul,pointu2,pointvl,pointv2){

returnzero((ul.x-u2.x)*(vl.y-v2.y)-(vl.x-v2.x)*(ul.y-u2.y));

)

〃判两直线垂直

intperpendicular(lineu,linev){

returnzero((u.a.x-u.b.x)*(v.a.x-v.b.x)+(u.a.y-u.b.y)*(v.a.y-v.b.y));

13

intperpendicular(pointul,pointu2,pointvl,pointv2){

returnzero((ul.x-u2.x)*(v1.x-v2.x)+(u1.y-u2.y)*(vl.y-v2.y));

)

〃判两线段相交,包括端点和部分重合

intintersect_in(lineu,linev){

if(!dots_inline(u.a,u.b,v.a)ll!dots_inline(u.a,u.b,v.b))

return!same_side(u.a,u.b,v)&&!same_side(v.a,v.b,u);

returndot_online_in(u.a,v)lldot_online_in(u.b,v)lldot_online_in(v.a,u)lldot_online_in(v.b,u);

)

intintersect_in(pointu1,pointu2,pointvl,pointv2){

if(!dots_inline(ul,u2,vl)ll!dots_inline(uI,u2,v2))

return!same_side(ul,u2,vl,v2)&&!same_side(v1,v2,u1,u2);

return

dot_online_in(u1,v1,v2)lldot_online_in(u2,v1,v2)Ildot_online_in(v1,u1,u2)Ildot_online_in(v2,u1,u

2)；

〃判两线段相交，不包括端点和部分重合

intintersect_ex(lineu,linev){

returnopposite_side(u.a,u.b,v)&&opposite_side(v.a,v.b,u);

)

intintersect_ex(pointul,pointu2,pointvl,pointv2){

returnopposite_side(ul,u2,vl,v2)&&opposite_side(vl,v2,ul,u2);

〃计算两直线交点，注意事先判断直线是否平行！

〃线段交点请另外判线段相交(同时还是要判断是否平行!)

pointintersection(lineu,linev){

pointret=u.a;

doublet=((u.a.x-v.a,x)*(v.a.y-v.b.y)-(u.a.y-v.a.y)*(v.a.x-v.b.x))

/((u.a.x-u.b.x)*(v.a.y-v.b.y)-(u.a.y-u.b.y)*(v.a.x-v.b.x));

ret.x+=(u.b.x-u.a.x)*t;

ret.y+=(u.b.y-u.a.y)*t;

returnret;

)

pointintersection(pointul,pointu2,pointvl,pointv2){

pointret=u1;

doublet=((ul.x-vl.x)*(vI.y-v2.y)-(ul.y-vl.y)*(vl.x-v2.x))

/((ul.x-u2.x)*(vl.y-v2.y)-(ul.y-u2.y)*(vl.x-v2.x));

ret.x+=(u2.x-ul.x)*t;

ret.y+=(u2.y-ul.y)*t;

returnret;

14

〃点到直线上的最近点

pointptoline(pointp,line1){

pointt=p;

t.x+=l.a.y-l.b.y,t.y+=l.b.x-l.a.x;

returnintersection(p,t,l.a,Lb);

)

pointptoline(pointp,pointIhpoint12){

pointt=p;

t.x+=l1.y-12.y,t.y+=12.x-l1.x;

returnintersection(p,t,l1,12);

)

〃点到直线距离

doubledisptoline(pointpjine1){

returnfabs(xmult(p,l.a,l.b))/distance(l.a,l.b);

)

doubledisptoline(pointp,point11,point12){

returnfabs(xmult(p,l1,12))/distance(l1,12);

)

doubledisptoline(doublex,doubley,doublexl,doubleyl,doublex2,doubley2){

returnfabs(xmult(x,y,xl,y1,x2,y2))/distance(x1,yI,x2,y2);

)

〃点到线段上的最近点

pointptoseg(pointp,line1){

pointt=p;

t.x+=l.a.y-l.b.y,t,y+=Lb.x-l.a.x;

if(xmult(l.a,t,p)*xmult(l.b,t,p)>eps)

returndistance(p,l.a)<distance(p,l.b)?l.a:l.b;

returnintersection(p,t,l.a,l.b);

)

pointptoseg(pointp,point11,point12){

pointt=p;

t.x+=l1.y-12.y,t.y+=12.x-l1.x;

if(xmult(ll,t,p)*xmult(12,t,p)>eps)

returndistance(p,l1)<distance(p,12)711:12;

returnintersection(p,t,l1,12);

)

〃点到线段距离

doubledisptoseg(pointp,line1){

pointt=p;

15

t.x+=l.a.y-l.b.y,t.y+=l.b.x-l.a.x;

if(xmult(l.a,t,p)*xmult(l.b,t,p)>eps)

returndistance(p,l.a)<distance(p,l.b)?distance(p,l.a):distance(p,l.b);

returnfabs(xmult(p,l.a,1.b))/distance(La,l.b);

)

doubledisptoseg(pointp,point11,point12){

pointt=p;

t.x+=ll.y-12.y,t.y+=12.x-l1.x;

if(xmult(ll,t,p)*xmult(12,t,p)>cps)

returndistance(p,ll)<distance(p,12)?distance(p,ll):distance(p,12);

returnfabs(xmult(p,l1,12))/distance(l1,12);

)

//矢量V以P为顶点逆时针旋转angle并放大scale倍

pointrotate(pointv,pointp,doubleangle,doublescale){

pointret=p;

v.x-=p.x,v.y-=p.y;

p.x=scale*cos(angle);

p.y=scale*sin(angle);

ret.x+=v.x*p.x-v.y*p.y;

ret.y+=v.x*p.y+v.y*p.x;

returnret;

)

1.6面积

#include<math.h>

structpoint{doublex,y;};

〃计算crossproduct(Pl-P0)x(P2-P0)

doublexmult(pointpl,pointp2,pointp0){

return(pLx-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(pl.y-pO.y);

)

doublexmult(doublexl,doubleyl,doublex2,doubley2,doublexO,doubley0){

return(x1-x0)*(y2-y0)-(x2-x0)*(y1-yO);

)

〃计算三角形面积，输入三顶点

doublearea_triangle(pointpl,pointp2,pointp3){

returnfabs(xmult(p1,p2,p3))/2;

)

doublearea_triangle(doublex1,doubleyhdoublex2,doubley2,doublex3,doubley3){

returnfabs(xmult(xl,yl,x2,y2,x3,y3))/2;

16

〃计算三角形面积，输入三边长

doublearea_triangle(doublea,doubleb,doublec){

doubles=(a+b+c)/2;

returnsqrt(s*(s-a)*(s-b)*(s-c));

)

〃计算多边形面积，顶点按顺时针或逆时针给出

doublearea_polygon(intn,point*p){

doublesl=0,s2=0;

inti;

fbr(i=0;i<n;i++)

sl+=p[(i+l)%n].y*p[i].x,s2+=p[(i+l)%n].y*p[(i+2)%n].x;

returnfabs(sl-s2)/2;

)

1.7球面

#include<math.h>

constdoublepi=acos(-l);

〃计算圆心角lat表示纬度,-90<=w<=90,lng表示经度

〃返回两点所在大圆劣弧对应圆心角,0<=angle<=pi

doubleangle(doubleIng1,doublelat1,doublelng2,doublelat2){

doubledlng=fabs(lngl-lng2)*pi/180;

while(dlng>=pi+pi)

dlng-=pi+pi;

if(dlng>pi)

dlng=pi+pi-dlng;

latl*=pi/180,lat2*=pi/180;

returnacos(cos(latl)*cos(lat2)*cos(dlng)4-sin(latl)*sin(lat2));

)

//计算距离,r为球半径

doubleline_dist(doubler,doubleIngl,doublelatl,doubleIng2,doublelat2){

doubledlng=fabs(lngl-lng2)*pi/l80;

while(dlng>=pi+pi)

dlng-=pi+pi;

if(dlng>pi)

dlng=pi+pi-dlng;

latl*=pi/180,lat2*=pi/180;

returnr*sqrt(2-2*(cos(latl)*cos(lat2)*cos(dlng)+sin(latl)*sin(lat2)));

)

17

〃计算球面距离J为球半径

inlinedoublesphere_dist(doubler,double1ng1,double1at1,doublelng2,doublelat2){

returnr*angle(lngl,latI,lng2,lat2);

)

1.8三角形

#include<math.h>

structpoint{doublex,y;};

structline{pointa,b;};

doubledistance(pointpl,pointp2){

returnsqrt((pl.x-p2.x)*(pl.x-p2.x)+(pl.y・p2.y)*(pl.y-p2.y));

)

pointintersection(lineu,linev){

pointret=u.a;

doublet=((u.a.x-v.a.x)*(v.a.y-v.b.y)-(u.a.y-v.a.y)*(v.a.x-v.b.x))

/((u.a.x-u.b.x)*(v.a.y-v.b.y)-(u.a.y-u.b.y)*(v.a.x-v.b.x));

ret.x+=(u.b.x-u.a.x)*t;

ret.y+=(u.b.y-u.a.y)*t;

returnret;

)

//外心

pointcircumcenter(pointa,pointb,pointc){

lineu,v;

u.a.x=(a.x+b.x)/2;

u.a.y=(a.y+b.y)/2;

u.b.x=u.a.x-a.y+b.y;

u.b.y=u.a.y+a.x-b.x;

v.a.x=(a.x+c.x)/2;

v.a.y=(a.y+c.y)/2;

v上.x=v.a.x・a.y+c.y;

v.b.y=v.a.y+a.x-c.x;

returnintersection(u,v);

)

〃内心

pointincenter(pointa,pointb,pointc){

lineu,v;

doublem,n;

u.a=a;

m=atan2(b.y-a.y,b.x-a.x);

18

n=atan2(c.y-a.y,c.x-a.x);

u.b.x=u.a.x+cos((m+n)/2);

u.b.y=u.a.y+sin((m+n)/2);

v.a=b;

m=atan2(a.y-b.y,a.x-b.x);

n=atan2(c.y-b.y,c.x-b.x);

v.b.x=v.a.x+cos((m+n)/2);

v.b.y=v.a.y+sin((m+n)/2);

returnintersection(u,v);

)

〃垂心

pointperpencenter(pointa,pointb,pointc){

lineu,v;

u.a=c;

u.b.x=u.a.x-a.y+b.y;

u.b.y=u.a.y+a.x-b.x;

v.a=b;

v.b.x=v.a.x-a.y+c.y;

v.b.y=v.a.y+a.x-c.x;

returnintersection(u,v);

)

〃重心

〃到三角形三顶点距离的平方和最小的点

〃三角形内到三边距离之积最大的点

pointbarycenter(pointa,pointb,pointc){

lineu,v;

u.a.x=(a.x+b.x)/2;

u.a.y=(a.y+b.y)/2;

u.b=c;

v.a.x=(a.x+c.x)/2;

v.a.y=(a.y+c.y)/2;

v.b=b;

returnintersection(u,v);

)