1、20Zhejiang UniversityICPC TeamRoutine Libraryby WishingBone (Dec. 2002)Last Update (Nov. 2004) by Riveria211、 几何 251.1 注意 251.2 几何公式 251.3 多边形 271.4 多边形切割 301.5 浮点函数 311.6 面积 361.7 球面 371.8 三角形 381.9 三维几何 401.10 凸包 471.11 网格 491.12 圆 491.13 整数函数 512 组合 542.1 组合公式 .542.2 排列组合生成 .542.3 生成 gray 码 .562.
2、4 置换(polya) .562.5 字典序全排列 .572.6 字典序组合 .573、结构 583.1 并查集 .583.2 堆 .593.3 线段树 .603.4 子段和 .653.5 子阵和 .654、 数论664.1 阶乘最后非 0 位 .664.2 模线性方程组 .674.3 素数 .684.4 欧拉函数 .695、 数值计算705.1 定积分计算(Romberg) .705.2 多项式求根(牛顿法) 725.3 周期性方程(追赶法) 73226、 .图论 NP 搜索746.1 最大团 .746.2 最大团(npi+pi-beta;相等应该是abs(a1-a2)pi+pi-eps;7
3、. 需要的话尽量使用 atan2,注意:atan2(0,0)=0,atan2(1,0)=pi/2,atan2(-1,0)=-pi/2,atan2(0,1)=0,atan2(0,-1)=pi.8. cross product = |u|*|v|*sin(a)dot product = |u|*|v|*cos(a)9. (P1-P0)x(P2-P0)结果的意义:正: 在顺时针(0,pi)内负: 在逆时针(0,pi)内0 : ,共线,夹角为 0 或 pi10. 误差限缺省使用 1e-8!1.2 几何公式三角形:1. 半周长 P=(a+b+c)/22. 面积 S=aHa/2=absin(C)/2=sq
4、rt(P(P-a)(P-b)(P-c)3. 中线 Ma=sqrt(2(b2+c2)-a2)/2=sqrt(b2+c2+2bccos(A)/24. 角平分线 Ta=sqrt(bc(b+c)2-a2)/(b+c)=2bccos(A/2)/(b+c)5. 高线 Ha=bsin(C)=csin(B)=sqrt(b2-(a2+b2-c2)/(2a)2)266. 内切圆半径 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. 外接圆
5、半径 R=abc/(4S)=a/(2sin(A)=b/(2sin(B)=c/(2sin(C)四边形:D1,D2 为对角线,M 对角线中点连线,A 为对角线夹角1. a2+b2+c2+d2=D12+D22+4M22. S=D1D2sin(A)/2(以下对圆的内接四边形)3. ac+bd=D1D24. S=sqrt(P-a)(P-b)(P-c)(P-d),P 为半周长正 n 边形:R 为外接圆半径,r 为内切圆半径1. 中心角 A=2PI/n2. 内角 C=(n-2)PI/n3. 边长 a=2sqrt(R2-r2)=2Rsin(A/2)=2rtan(A/2)4. 面积 S=nar/2=nr2tan
6、(A/2)=nR2sin(A)/2=na2/(4tan(A/2)圆:1. 弧长 l=rA2. 弦长 a=2sqrt(2hr-h2)=2rsin(A/2)3. 弓形高 h=r-sqrt(r2-a2/4)=r(1-cos(A/2)=atan(A/4)/24. 扇形面积 S1=rl/2=r2A/25. 弓形面积 S2=(rl-a(r-h)/2=r2(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,l 为斜
7、高,p 为底面周长3. 全面积 T=S+A棱台:1. 体积 V=(A1+A2+sqrt(A1A2)h/3,A1.A2 为上下底面积,h 为高(以下为正棱台)2. 侧面积 S=(p1+p2)l/2,p1.p2 为上下底面周长,l 为斜高3. 全面积 T=S+A1+A227圆柱:1. 侧面积 S=2PIrh2. 全面积 T=2PIr(h+r)3. 体积 V=PIr2h圆锥:1. 母线 l=sqrt(h2+r2)2. 侧面积 S=PIrl3. 全面积 T=PIr(l+r)4. 体积 V=PIr2h/3圆台:1. 母线 l=sqrt(h2+(r1-r2)2)2. 侧面积 S=PI(r1+r2)l3.
8、全面积 T=PIr1(l+r1)+PIr2(l+r2)4. 体积 V=PI(r12+r22+r1r2)h/3球:1. 全面积 T=4PIr22. 体积 V=4PIr3/3球台:1. 侧面积 S=2PIrh2. 全面积 T=PI(2rh+r12+r22)3. 体积 V=PIh(3(r12+r22)+h2)/6球扇形:1. 全面积 T=PIr(2h+r0),h 为球冠高,r0 为球冠底面半径2. 体积 V=2PIr2h/31.3 多边形#include #include #define MAXN 1000#define offset 10000#define eps 1e-8#define zer
9、o(x) (x)0?(x):-(x)eps?1:(x)eps)t=barycenter(p0,pi,pi+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;1.4 多边形切割/多边形切割/可用于半平面交#define MAXN 100#define eps 1e-8#define zero(x) (x)0?(x):-(x)eps;point intersection(point u1,point u2,point v1,point v2)point ret=u1;double
10、 t=(u1.x-v1.x)*(v1.y-v2.y)-(u1.y-v1.y)*(v1.x-v2.x)/(u1.x-u2.x)*(v1.y-v2.y)-(u1.y-u2.y)*(v1.x-v2.x);ret.x+=(u2.x-u1.x)*t;ret.y+=(u2.y-u1.y)*t;return ret;/将多边形沿 l1,l2 确定的直线切割在 side 侧切割,保证 l1,l2,side 不共线void polygon_cut(intint m=0,i;for (i=0;i#define eps 1e-8#define zero(x) (x)0?(x):-(x)eps;int same_si
11、de(point p1,point p2,point l1,point l2)return xmult(l1,p1,l2)*xmult(l1,p2,l2)eps;/判两点在线段异侧,点在线段上返回 0int opposite_side(point p1,point p2,line l)return xmult(l.a,p1,l.b)*xmult(l.a,p2,l.b)eps)return distance(p,l.a)eps)return distance(p,l1)eps)36return distance(p,l.a)eps)return distance(p,l1)struct poin
12、tdouble x,y;/计算 cross product (P1-P0)x(P2-P0)double xmult(point p1,point p2,point p0)return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);double xmult(double x1,double y1,double x2,double y2,double x0,double y0)return (x1-x0)*(y2-y0)-(x2-x0)*(y1-y0);/计算三角形面积,输入三顶点double area_triangle(point p1,poin
13、t p2,point p3)return fabs(xmult(p1,p2,p3)/2;double area_triangle(double x1,double y1,double x2,double y2,double x3,double y3)return fabs(xmult(x1,y1,x2,y2,x3,y3)/2;/计算三角形面积,输入三边长37double area_triangle(double a,double b,double c)double s=(a+b+c)/2;return sqrt(s*(s-a)*(s-b)*(s-c);/计算多边形面积,顶点按顺时针或逆时针给出
14、double area_polygon(int n,point* p)double s1=0,s2=0;int i;for (i=0;iconst double pi=acos(-1);/计算圆心角 lat 表示纬度,-90=pi+pi)dlng-=pi+pi;if (dlngpi)dlng=pi+pi-dlng;lat1*=pi/180,lat2*=pi/180;return acos(cos(lat1)*cos(lat2)*cos(dlng)+sin(lat1)*sin(lat2);/计算距离,r 为球半径double line_dist(double r,double lng1,doub
15、le lat1,double lng2,double lat2)double dlng=fabs(lng1-lng2)*pi/180;while (dlng=pi+pi)dlng-=pi+pi;if (dlngpi)dlng=pi+pi-dlng;lat1*=pi/180,lat2*=pi/180;return r*sqrt(2-2*(cos(lat1)*cos(lat2)*cos(dlng)+sin(lat1)*sin(lat2);/计算球面距离,r 为球半径inline double sphere_dist(double r,double lng1,double lat1,double l
16、ng2,double lat2)38return r*angle(lng1,lat1,lng2,lat2);1.8 三角形#include struct pointdouble x,y;struct linepoint a,b;double distance(point p1,point p2)return sqrt(p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y);point intersection(line u,line v)point ret=u.a;double t=(u.a.x-v.a.x)*(v.a.y-v.b.y)-(u.a.y-v.
17、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;return ret;/外心point circumcenter(point a,point b,point c)line u,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
18、.b.x=v.a.x-a.y+c.y;v.b.y=v.a.y+a.x-c.x;return intersection(u,v);/内心point incenter(point a,point b,point c)line u,v;double m,n;u.a=a;m=atan2(b.y-a.y,b.x-a.x);n=atan2(c.y-a.y,c.x-a.x);u.b.x=u.a.x+cos(m+n)/2);39u.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+
19、cos(m+n)/2);v.b.y=v.a.y+sin(m+n)/2);return intersection(u,v);/垂心point perpencenter(point a,point b,point c)line u,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;return intersection(u,v);/重心/到三角形三顶点距离的平方和最小的点/三角形内到三边距离之积最大的点point barycenter(point a,point
20、 b,point c)line u,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;return intersection(u,v);/费马点/到三角形三顶点距离之和最小的点point fermentpoint(point a,point b,point c)point u,v;double step=fabs(a.x)+fabs(a.y)+fabs(b.x)+fabs(b.y)+fabs(c.x)+fabs(c.y);int i,j,k;u.x=(a.x+b.x+c.x
21、)/3;u.y=(a.y+b.y+c.y)/3;while (step1e-10)40for (k=0;kdistance(v,a)+distance(v,b)+distance(v,c)u=v;return u;1.9 三维几何/三维几何函数库#include #define eps 1e-8#define zero(x) (x)0?(x):-(x)epsint dot_inplane_ex(point3 p,point3 s1,point3 s2,point3 s3)return dot_inplane_in(p,s1,s2,s3)/判两点在线段同侧,点在线段上返回 0,不共面无意义int
22、 same_side(point3 p1,point3 p2,line3 l)return dmult(xmult(subt(l.a,l.b),subt(p1,l.b),xmult(subt(l.a,l.b),subt(p2,l.b)eps;int same_side(point3 p1,point3 p2,point3 l1,point3 l2)return dmult(xmult(subt(l1,l2),subt(p1,l2),xmult(subt(l1,l2),subt(p2,l2)eps;/判两点在线段异侧,点在线段上返回 0,不共面无意义int opposite_side(point
23、3 p1,point3 p2,line3 l)return dmult(xmult(subt(l.a,l.b),subt(p1,l.b),xmult(subt(l.a,l.b),subt(p2,l.b)eps;43int same_side(point3 p1,point3 p2,point3 s1,point3 s2,point3 s3)return dmult(pvec(s1,s2,s3),subt(p1,s1)*dmult(pvec(s1,s2,s3),subt(p2,s1)eps;/判两点在平面异侧,点在平面上返回 0int opposite_side(point3 p1,point3
24、 p2,plane3 s)return dmult(pvec(s),subt(p1,s.a)*dmult(pvec(s),subt(p2,s.a)#define eps 1e-848#define zero(x) (x)0?(x):-(x)0?1:-1):(ret0?1:-1);void _graham(int n,point* p,intfor (p1=p2=p0,i=1;ieps|(zero(p1.y-pi.y)p2.x/=n,p2.y/=n;pk=p0,p0=p1;qsort(p+1,n-1,sizeof(point),graham_cp);for (ch0=p0,ch1=p1,ch2=
25、p2,s=i=3;i2;int gcd(int a,int b)return b?gcd(b,a%b):a;/多边形上的网格点个数int grid_onedge(int n,point* p)int i,ret=0;for (i=0;i#define eps 1e-8struct pointdouble x,y;double xmult(point p1,point p2,point p0)return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);double distance(point p1,point p2)return sqrt(p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y);double disptoline(point p,point l1,point l2)return fabs(xmult(p,l1,l2)/distance(l1,l2);point intersection(point u1,point u2,point v1,point v2)point ret=u1;double t=(u1.x-v1.x)*(v1.y-v2.y)-(u1.y-v1.y)*(v1.x-v2.x)
