C语言算法整理

Mr_扛扛 贡献于2012-08-05

作者 Riveria  创建于2004-11-10 04:34:00   修改者雨林木风  修改于2010-04-14 19:36:00字数64914

文档摘要:1. 注意舍入方式(0.5的舍入方向);防止输出-0. 2. 几何题注意多测试不对称数据. 3. 整数几何注意xmult和dmult是否会出界; 符点几何注意eps的使用. 4. 避免使用斜率;注意除数是否会为0. 5. 公式一定要化简后再代入.
关键词:

ACM Routine Library by Wood_K (Apr. 2010) Last Update (Apr. 2010) by Wood_K 1、 几何 25 1.1 注意 25 1.2 几何公式 25 1.3 多边形 27 1.5 浮点函数 31 1.6 面积 36 1.7 球面 37 1.8 三角形 38 1.9 三维几何 40 1.10 凸包 47 1.11 网格 49 1.12 圆 49 1.13 整数函数 51 2、 组合 54 2.1 组合公式 54 2.2 排列组合生成 54 2.3 生成gray码 56 2.4 置换(polya) 56 2.5 字典序全排列 57 2.6 字典序组合 57 3、 结构 58 3.1 并查集 58 3.2 堆 59 3.3 线段树 60 3.4 子段和 65 3.5 子阵和 65 4、 数论 66 4.1 阶乘最后非0位 66 4.2 模线性方程组 67 4.3 素数 68 4.4 欧拉函数 69 5、 数值计算 70 5.1 定积分计算(Romberg) 70 5.2 多项式求根(牛顿法) 72 5.3 周期性方程(追赶法) 73 6、 图论—NP搜索 74 6.1 最大团 74 6.2 最大团(n<64)(faster) 75 7、 图论—连通性 77 7.1 无向图关键点(dfs邻接阵) 77 7.2 无向图关键边(dfs邻接阵) 78 7.3 无向图的块(bfs邻接阵) 79 7.4 无向图连通分支(dfs/bfs邻接阵) 80 7.5 有向图强连通分支(dfs/bfs邻接阵) 81 7.6 有向图最小点基(邻接阵) 82 8、 图论—匹配 83 8.1 二分图最大匹配(hungary邻接表) 83 8.2 二分图最大匹配(hungary邻接阵) 84 8.3 二分图最大匹配(hungary正向表) 84 8.4二分图最佳匹配(kuhn_munkras邻接阵) 85 8.5 一般图匹配(邻接表) 86 8.6 一般图匹配(邻接阵) 87 8.7 一般图匹配(正向表) 87 9、 图论—网络流 88 9.1 最大流(邻接阵) 88 9.2 上下界最大流(邻接阵) 89 9.3 上下界最小流(邻接阵) 90 9.4 最大流无流量(邻接阵) 91 9.5 最小费用最大流(邻接阵) 91 10、 图论—应用 92 10.1 欧拉回路(邻接阵) 92 10.2 树的前序表转化 93 10.3 树的优化算法 94 10.4 拓扑排序(邻接阵) 95 10.5 最佳边割集 96 10.6 最佳点割集 97 10.7 最小边割集 98 10.8 最小点割集 99 10.9 最小路径覆盖 101 11、 图论—支撑树 101 11.1 最小生成树(kruskal邻接表) 101 11.2 最小生成树(kruskal正向表) 103 11.3 最小生成树(prim+binary_heap邻接表) 104 11.4 最小生成树(prim+binary_heap正向表) 105 11.5 最小生成树(prim+mapped_heap邻接表) 106 11.6 最小生成树(prim+mapped_heap正向表) 108 11.7 最小生成树(prim邻接阵) 109 11.8 最小树形图(邻接阵) 109 12、 图论—最短路径 111 12.1 最短路径(单源bellman_ford邻接阵) 111 12.2 最短路径(单源dijkstra+bfs邻接表) 111 12.3 最短路径(单源dijkstra+bfs正向表) 112 12.4 最短路径(单源dijkstra+binary_heap邻接表) 113 12.5 最短路径(单源dijkstra+binary_heap正向表) 114 12.6 最短路径(单源dijkstra+mapped_heap邻接表) 115 12.7 最短路径(单源dijkstra+mapped_heap正向表) 116 12.8 最短路径(单源dijkstra邻接阵) 117 12.9 最短路径(多源floyd_warshall邻接阵) 118 13、 应用 118 13.1 Joseph问题 118 13.2 N皇后构造解 119 13.3 布尔母函数 120 13.4 第k元素 120 13.5 幻方构造 121 13.6 模式匹配(kmp) 122 13.7 逆序对数 123 13.8 字符串最小表示 123 13.9 最长公共单调子序列 124 13.10 最长子序列 125 13.11 最大子串匹配 126 13.12 最大子段和 127 13.13 最大子阵和 127 14、 其它 128 14.1 大数(只能处理正数) 128 14.2 分数 134 14.3 矩阵 136 14.4 线性方程组 138 14.5 线性相关 140 14.6 日期 140 1、 几何 1.1 注意 1. 注意舍入方式(0.5的舍入方向);防止输出-0. 2. 几何题注意多测试不对称数据. 3. 整数几何注意xmult和dmult是否会出界; 符点几何注意eps的使用. 4. 避免使用斜率;注意除数是否会为0. 5. 公式一定要化简后再代入. 6. 判断同一个2*PI域内两角度差应该是 abs(a1-a2)pi+pi-beta; 相等应该是 abs(a1-a2)pi+pi-eps; 7. 需要的话尽量使用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或pi 10. 误差限缺省使用1e-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(b^2+c^2)-a^2)/2=sqrt(b^2+c^2+2bccos(A))/2 4. 角平分线 Ta=sqrt(bc((b+c)^2-a^2))/(b+c)=2bccos(A/2)/(b+c) 5. 高线 Ha=bsin(C)=csin(B)=sqrt(b^2-((a^2+b^2-c^2)/(2a))^2) 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. a^2+b^2+c^2+d^2=D1^2+D2^2+4M^2 2. S=D1D2sin(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(R^2-r^2)=2Rsin(A/2)=2rtan(A/2) 4. 面积 S=nar/2=nr^2tan(A/2)=nR^2sin(A)/2=na^2/(4tan(A/2)) 圆: 1. 弧长 l=rA 2. 弦长 a=2sqrt(2hr-h^2)=2rsin(A/2) 3. 弓形高 h=r-sqrt(r^2-a^2/4)=r(1-cos(A/2))=atan(A/4)/2 4. 扇形面积 S1=rl/2=r^2A/2 5. 弓形面积 S2=(rl-a(r-h))/2=r^2(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为斜高,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+A2 圆柱: 1. 侧面积 S=2PIrh 2. 全面积 T=2PIr(h+r) 3. 体积 V=PIr^2h 圆锥: 1. 母线 l=sqrt(h^2+r^2) 2. 侧面积 S=PIrl 3. 全面积 T=PIr(l+r) 4. 体积 V=PIr^2h/3 圆台: 1. 母线 l=sqrt(h^2+(r1-r2)^2) 2. 侧面积 S=PI(r1+r2)l 3. 全面积 T=PIr1(l+r1)+PIr2(l+r2) 4. 体积 V=PI(r1^2+r2^2+r1r2)h/3 球: 1. 全面积 T=4PIr^2 2. 体积 V=4PIr^3/3 球台: 1. 侧面积 S=2PIrh 2. 全面积 T=PI(2rh+r1^2+r2^2) 3. 体积 V=PIh(3(r1^2+r2^2)+h^2)/6 球扇形: 1. 全面积 T=PIr(2h+r0),h为球冠高,r0为球冠底面半径 2. 体积 V=2PIr^2h/3 1.3 多边形 #include #include #define MAXN 1000 #define offset 10000 #define eps 1e-8 #define zero(x) (((x)>0?(x):-(x))eps?1:((x)<-eps?2:0)) struct point{double x,y;}; struct line{point a,b;}; 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); } //判定凸多边形,顶点按顺时针或逆时针给出,允许相邻边共线 int is_convex(int n,point* p){ int i,s[3]={1,1,1}; for (i=0;ieps){ 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; } 1.4 浮点函数 //浮点几何函数库 #include #define eps 1e-8 #define zero(x) (((x)>0?(x):-(x))eps; } int same_side(point p1,point p2,point l1,point l2){ return xmult(l1,p1,l2)*xmult(l1,p2,l2)>eps; } //判两点在线段异侧,点在线段上返回0 int opposite_side(point p1,point p2,line l){ return xmult(l.a,p1,l.b)*xmult(l.a,p2,l.b)<-eps; } int opposite_side(point p1,point p2,point l1,point l2){ return xmult(l1,p1,l2)*xmult(l1,p2,l2)<-eps; } //判两直线平行 int parallel(line u,line v){ return zero((u.a.x-u.b.x)*(v.a.y-v.b.y)-(v.a.x-v.b.x)*(u.a.y-u.b.y)); } int parallel(point u1,point u2,point v1,point v2){ return zero((u1.x-u2.x)*(v1.y-v2.y)-(v1.x-v2.x)*(u1.y-u2.y)); } //判两直线垂直 int perpendicular(line u,line v){ return zero((u.a.x-u.b.x)*(v.a.x-v.b.x)+(u.a.y-u.b.y)*(v.a.y-v.b.y)); } int perpendicular(point u1,point u2,point v1,point v2){ return zero((u1.x-u2.x)*(v1.x-v2.x)+(u1.y-u2.y)*(v1.y-v2.y)); } //判两线段相交,包括端点和部分重合 int intersect_in(line u,line v){ if (!dots_inline(u.a,u.b,v.a)||!dots_inline(u.a,u.b,v.b)) return !same_side(u.a,u.b,v)&&!same_side(v.a,v.b,u); return dot_online_in(u.a,v)||dot_online_in(u.b,v)||dot_online_in(v.a,u)||dot_online_in(v.b,u); } int intersect_in(point u1,point u2,point v1,point v2){ if (!dots_inline(u1,u2,v1)||!dots_inline(u1,u2,v2)) return !same_side(u1,u2,v1,v2)&&!same_side(v1,v2,u1,u2); return dot_online_in(u1,v1,v2)||dot_online_in(u2,v1,v2)||dot_online_in(v1,u1,u2)||dot_online_in(v2,u1,u2); } //判两线段相交,不包括端点和部分重合 int intersect_ex(line u,line v){ return opposite_side(u.a,u.b,v)&&opposite_side(v.a,v.b,u); } int intersect_ex(point u1,point u2,point v1,point v2){ return opposite_side(u1,u2,v1,v2)&&opposite_side(v1,v2,u1,u2); } //计算两直线交点,注意事先判断直线是否平行! //线段交点请另外判线段相交(同时还是要判断是否平行!) 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.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 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)) /((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; } //点到直线上的最近点 point ptoline(point p,line l){ point t=p; t.x+=l.a.y-l.b.y,t.y+=l.b.x-l.a.x; return intersection(p,t,l.a,l.b); } point ptoline(point p,point l1,point l2){ point t=p; t.x+=l1.y-l2.y,t.y+=l2.x-l1.x; return intersection(p,t,l1,l2); } //点到直线距离 double disptoline(point p,line l){ return fabs(xmult(p,l.a,l.b))/distance(l.a,l.b); } double disptoline(point p,point l1,point l2){ return fabs(xmult(p,l1,l2))/distance(l1,l2); } double disptoline(double x,double y,double x1,double y1,double x2,double y2){ return fabs(xmult(x,y,x1,y1,x2,y2))/distance(x1,y1,x2,y2); } //点到线段上的最近点 point ptoseg(point p,line l){ point t=p; 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) return distance(p,l.a)eps) return distance(p,l1)eps) return distance(p,l.a)eps) return distance(p,l1) struct point{double 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,point 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; } //计算三角形面积,输入三边长 double area_triangle(double a,double b,double c){ double s=(a+b+c)/2; return sqrt(s*(s-a)*(s-b)*(s-c)); } //计算多边形面积,顶点按顺时针或逆时针给出 double area_polygon(int n,point* p){ double s1=0,s2=0; int i; for (i=0;i const double pi=acos(-1); //计算圆心角lat表示纬度,-90<=w<=90,lng表示经度 //返回两点所在大圆劣弧对应圆心角,0<=angle<=pi double angle(double lng1,double lat1,double lng2,double lat2){ double dlng=fabs(lng1-lng2)*pi/180; while (dlng>=pi+pi) dlng-=pi+pi; if (dlng>pi) 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,double lat1,double lng2,double lat2){ double dlng=fabs(lng1-lng2)*pi/180; while (dlng>=pi+pi) dlng-=pi+pi; if (dlng>pi) 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 lng2,double lat2){ return r*angle(lng1,lat1,lng2,lat2); } 1.7 三角形 #include struct point{double x,y;}; struct line{point 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.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.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); 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); 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 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)/3; u.y=(a.y+b.y+c.y)/3; while (step>1e-10) for (k=0;k<10;step/=2,k++) for (i=-1;i<=1;i++) for (j=-1;j<=1;j++){ v.x=u.x+step*i; v.y=u.y+step*j; if (distance(u,a)+distance(u,b)+distance(u,c)>distance(v,a)+distance(v,b)+distance(v,c)) u=v; } return u; } 1.8 三维几何 //三维几何函数库 #include #define eps 1e-8 #define zero(x) (((x)>0?(x):-(x))eps&& vlen(xmult(subt(p,s.b),subt(p,s.c)))>eps&&vlen(xmult(subt(p,s.c),subt(p,s.a)))>eps; } int dot_inplane_ex(point3 p,point3 s1,point3 s2,point3 s3){ return dot_inplane_in(p,s1,s2,s3)&&vlen(xmult(subt(p,s1),subt(p,s2)))>eps&& vlen(xmult(subt(p,s2),subt(p,s3)))>eps&&vlen(xmult(subt(p,s3),subt(p,s1)))>eps; } //判两点在线段同侧,点在线段上返回0,不共面无意义 int 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(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 opposite_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 same_side(point3 p1,point3 p2,plane3 s){ return dmult(pvec(s),subt(p1,s.a))*dmult(pvec(s),subt(p2,s.a))>eps; } int 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; } //判两点在平面异侧,点在平面上返回0 int opposite_side(point3 p1,point3 p2,plane3 s){ return dmult(pvec(s),subt(p1,s.a))*dmult(pvec(s),subt(p2,s.a))<-eps; } int opposite_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; } //判两直线平行 int parallel(line3 u,line3 v){ return vlen(xmult(subt(u.a,u.b),subt(v.a,v.b))) #define eps 1e-8 #define zero(x) (((x)>0?(x):-(x))0?1:-1):(ret>0?1:-1); } void _graham(int n,point* p,int& s,point* ch){ int i,k=0; for (p1=p2=p[0],i=1;ieps||(zero(p1.y-p[i].y)&&p1.x>p[i].x)) p1=p[k=i]; p2.x/=n,p2.y/=n; p[k]=p[0],p[0]=p1; qsort(p+1,n-1,sizeof(point),graham_cp); for (ch[0]=p[0],ch[1]=p[1],ch[2]=p[2],s=i=3;i2&&xmult(ch[s-2],p[i],ch[s-1])<-eps;s--); } //构造凸包接口函数,传入原始点集大小n,点集p(p原有顺序被打乱!) //返回凸包大小,凸包的点在convex中 //参数maxsize为1包含共线点,为0不包含共线点,缺省为1 //参数clockwise为1顺时针构造,为0逆时针构造,缺省为1 //在输入仅有若干共线点时算法不稳定,可能有此类情况请另行处理! //不能去掉点集中重合的点 int graham(int n,point* p,point* convex,int maxsize=1,int dir=1){ point* temp=new point[n]; int s,i; _graham(n,p,s,temp); for (convex[0]=temp[0],n=1,i=(dir?1:(s-1));dir?(i0?(x):-(x)) struct point{int x,y;}; 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-8 struct point{double 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)) /((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; } //判直线和圆相交,包括相切 int intersect_line_circle(point c,double r,point l1,point l2){ return disptoline(c,l1,l2)-eps||t2>-eps; t.x+=l1.y-l2.y; t.y+=l2.x-l1.x; return xmult(l1,c,t)*xmult(l2,c,t)fabs(r1-r2)-eps; } //计算圆上到点p最近点,如p与圆心重合,返回p本身 point dot_to_circle(point c,double r,point p){ point u,v; if (distance(p,c)0?1:(((a)<0?-1:0))) struct point{int x,y;}; struct line{point a,b;}; //计算cross product (P1-P0)x(P2-P0) int xmult(point p1,point p2,point p0){ return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y); } int xmult(int x1,int y1,int x2,int y2,int x0,int y0){ return (x1-x0)*(y2-y0)-(x2-x0)*(y1-y0); } //计算dot product (P1-P0).(P2-P0) int dmult(point p1,point p2,point p0){ return (p1.x-p0.x)*(p2.x-p0.x)+(p1.y-p0.y)*(p2.y-p0.y); } int dmult(int x1,int y1,int x2,int y2,int x0,int y0){ return (x1-x0)*(x2-x0)+(y1-y0)*(y2-y0); } //判三点共线 int dots_inline(point p1,point p2,point p3){ return !xmult(p1,p2,p3); } int dots_inline(int x1,int y1,int x2,int y2,int x3,int y3){ return !xmult(x1,y1,x2,y2,x3,y3); } //判点是否在线段上,包括端点和部分重合 int dot_online_in(point p,line l){ return !xmult(p,l.a,l.b)&&(l.a.x-p.x)*(l.b.x-p.x)<=0&&(l.a.y-p.y)*(l.b.y-p.y)<=0; } int dot_online_in(point p,point l1,point l2){ return !xmult(p,l1,l2)&&(l1.x-p.x)*(l2.x-p.x)<=0&&(l1.y-p.y)*(l2.y-p.y)<=0; } int dot_online_in(int x,int y,int x1,int y1,int x2,int y2){ return !xmult(x,y,x1,y1,x2,y2)&&(x1-x)*(x2-x)<=0&&(y1-y)*(y2-y)<=0; } //判点是否在线段上,不包括端点 int dot_online_ex(point p,line l){ return dot_online_in(p,l)&&(p.x!=l.a.x||p.y!=l.a.y)&&(p.x!=l.b.x||p.y!=l.b.y); } int dot_online_ex(point p,point l1,point l2){ return dot_online_in(p,l1,l2)&&(p.x!=l1.x||p.y!=l1.y)&&(p.x!=l2.x||p.y!=l2.y); } int dot_online_ex(int x,int y,int x1,int y1,int x2,int y2){ return dot_online_in(x,y,x1,y1,x2,y2)&&(x!=x1||y!=y1)&&(x!=x2||y!=y2); } //判两点在直线同侧,点在直线上返回0 int same_side(point p1,point p2,line l){ return sign(xmult(l.a,p1,l.b))*xmult(l.a,p2,l.b)>0; } int same_side(point p1,point p2,point l1,point l2){ return sign(xmult(l1,p1,l2))*xmult(l1,p2,l2)>0; } //判两点在直线异侧,点在直线上返回0 int opposite_side(point p1,point p2,line l){ return sign(xmult(l.a,p1,l.b))*xmult(l.a,p2,l.b)<0; } int opposite_side(point p1,point p2,point l1,point l2){ return sign(xmult(l1,p1,l2))*xmult(l1,p2,l2)<0; } //判两直线平行 int parallel(line u,line v){ return (u.a.x-u.b.x)*(v.a.y-v.b.y)==(v.a.x-v.b.x)*(u.a.y-u.b.y); } int parallel(point u1,point u2,point v1,point v2){ return (u1.x-u2.x)*(v1.y-v2.y)==(v1.x-v2.x)*(u1.y-u2.y); } //判两直线垂直 int perpendicular(line u,line v){ return (u.a.x-u.b.x)*(v.a.x-v.b.x)==-(u.a.y-u.b.y)*(v.a.y-v.b.y); } int perpendicular(point u1,point u2,point v1,point v2){ return (u1.x-u2.x)*(v1.x-v2.x)==-(u1.y-u2.y)*(v1.y-v2.y); } //判两线段相交,包括端点和部分重合 int intersect_in(line u,line v){ if (!dots_inline(u.a,u.b,v.a)||!dots_inline(u.a,u.b,v.b)) return !same_side(u.a,u.b,v)&&!same_side(v.a,v.b,u); return dot_online_in(u.a,v)||dot_online_in(u.b,v)||dot_online_in(v.a,u)||dot_online_in(v.b,u); } int intersect_in(point u1,point u2,point v1,point v2){ if (!dots_inline(u1,u2,v1)||!dots_inline(u1,u2,v2)) return !same_side(u1,u2,v1,v2)&&!same_side(v1,v2,u1,u2); return dot_online_in(u1,v1,v2)||dot_online_in(u2,v1,v2)||dot_online_in(v1,u1,u2)||dot_online_in(v2,u1,u2); } //判两线段相交,不包括端点和部分重合 int intersect_ex(line u,line v){ return opposite_side(u.a,u.b,v)&&opposite_side(v.a,v.b,u); } int intersect_ex(point u1,point u2,point v1,point v2){ return opposite_side(u1,u2,v1,v2)&&opposite_side(v1,v2,u1,u2); } 2、 组合 2.1 组合公式 1. C(m,n)=C(m,m-n) 2. C(m,n)=C(m-1,n)+C(m-1,n-1) derangement D(n) = n!(1 - 1/1! + 1/2! - 1/3! + ... + (-1)^n/n!) = (n-1)(D(n-2) - D(n-1)) Q(n) = D(n) + D(n-1) 求和公式,k = 1..n 1. sum( k ) = n(n+1)/2 2. sum( 2k-1 ) = n^2 3. sum( k^2 ) = n(n+1)(2n+1)/6 4. sum( (2k-1)^2 ) = n(4n^2-1)/3 5. sum( k^3 ) = (n(n+1)/2)^2 6. sum( (2k-1)^3 ) = n^2(2n^2-1) 7. sum( k^4 ) = n(n+1)(2n+1)(3n^2+3n-1)/30 8. sum( k^5 ) = n^2(n+1)^2(2n^2+2n-1)/12 9. sum( k(k+1) ) = n(n+1)(n+2)/3 10. sum( k(k+1)(k+2) ) = n(n+1)(n+2)(n+3)/4 12. sum( k(k+1)(k+2)(k+3) ) = n(n+1)(n+2)(n+3)(n+4)/5 2.2 排列组合生成 //gen_perm产生字典序排列P(n,m) //gen_comb产生字典序组合C(n,m) //gen_perm_swap产生相邻位对换全排列P(n,n) //产生元素用1..n表示 //dummy为产生后调用的函数,传入a[]和n,a[0]..a[n-1]为一次产生的结果 #define MAXN 100 int count; #include void dummy(int* a,int n){ int i; cout<=0;k*=n-(i--)) for (j=i+1;j=0;i--) p[i]=t%(n-i),t/=n-i; for (i=n-1;i;i--) for (j=i-1;j>=0;j--) if (p[j]<=p[i]) p[i]++; } 2.6 字典序组合 //字典序组合与序号的转换 //comb为组合数C(n,m),必要时换成大数,注意处理C(n,m)=0|n(k=comb(n-j,m-i-1));t-=k,j++); } 数论 4.1 阶乘最后非0位 //求阶乘最后非零位,复杂度O(nlogn) //返回该位,n以字符串方式传入 #include #define MAXN 10000 int lastdigit(char* buf){ const int mod[20]={1,1,2,6,4,2,2,4,2,8,4,4,8,4,6,8,8,6,8,2}; int len=strlen(buf),a[MAXN],i,c,ret=1; if (len==1) return mod[buf[0]-'0']; for (i=0;i=0;i--) c=c*10+a[i],a[i]=c/5,c%=5; } return ret+ret%2*5; } 4.4 欧拉函数 int gcd(int a,int b){ return b?gcd(b,a%b):a; } inline int lcm(int a,int b){ return a/gcd(a,b)*b; } //求1..n-1中与n互质的数的个数 int eular(int n){ int ret=1,i; for (i=2;i*i<=n;i++) if (n%i==0){ n/=i,ret*=i-1; while (n%i==0) n/=i,ret*=i; } if (n>1) ret*=n-1; return ret; } 3、 数值计算 5.2 多项式求根(牛顿法) /* 牛顿法解多项式的根 输入:多项式系数c[],多项式度数n,求在[a,b]间的根 输出:根 要求保证[a,b]间有根 */ double fabs( double x ) { return (x<0)? -x : x; } double f(int m, double c[], double x) { int i; double p = c[m]; for (i=m; i>0; i--) p = p*x + c[i-1]; return p; } int newton(double x0, double *r, double c[], double cp[], int n, double a, double b, double eps) { int MAX_ITERATION = 1000; int i = 1; double x1, x2, fp, eps2 = eps/10.0; x1 = x0; while (i < MAX_ITERATION) { x2 = f(n, c, x1); fp = f(n-1, cp, x1); if ((fabs(fp)<0.000000001) && (fabs(x2)>1.0)) return 0; x2 = x1 - x2/fp; if (fabs(x1-x2)b) return 0; *r = x2; return 1; } x1 = x2; i++; } return 0; } double Polynomial_Root(double c[], int n, double a, double b, double eps) { double *cp; int i; double root; cp = (double *)calloc(n, sizeof(double)); for (i=n-1; i>=0; i--) { cp[i] = (i+1)*c[i+1]; } if (a>b) { root = a; a = b; b = root; } if ((!newton(a, &root, c, cp, n, a, b, eps)) && (!newton(b, &root, c, cp, n, a, b, eps))) newton((a+b)*0.5, &root, c, cp, n, a, b, eps); free(cp); if (fabs(root) max) { max = size; for (i = 0; i < size; ++ i) res[i] = rr[i]; bb = 1; } } int maxclique(int n, int mat[][MAXN], int *ret) { int max = 0, bb, c[MAXN], i, j; int vn, v[MAXN], rr[MAXN]; for (c[i = n - 1] = 0; i >= 0; -- i) { for (vn = 0, j = i + 1; j < n; ++ j) if (mat[i][j]) v[vn ++] = j; bb = 0; rr[0] = i; clique(vn, v, mat, 1, max, bb, ret, rr, c); c[i] = max; } return max; } 6.2 最大团(n<64)(faster) /** * WishingBone's ACM/ICPC Routine Library * * maximum clique solver */ #include using std::vector; // clique solver calculates both size and consitution of maximum clique // uses bit operation to accelerate searching // graph size limit is 63, the graph should be undirected // can optimize to calculate on each component, and sort on vertex degrees // can be used to solve maximum independent set class clique { public: static const long long ONE = 1; static const long long MASK = (1 << 21) - 1; char* bits; int n, size, cmax[63]; long long mask[63], cons; // initiate lookup table clique() { bits = new char[1 << 21]; bits[0] = 0; for (int i = 1; i < 1 << 21; ++i) bits[i] = bits[i >> 1] + (i & 1); } ~clique() { delete bits; } // search routine bool search(int step, int size, long long more, long long con); // solve maximum clique and return size int sizeClique(vector >& mat); // solve maximum clique and return constitution vector consClique(vector >& mat); }; // search routine // step is node id, size is current solution, more is available mask, cons is constitution mask bool clique::search(int step, int size, long long more, long long cons) { if (step >= n) { // a new solution reached this->size = size; this->cons = cons; return true; } long long now = ONE << step; if ((now & more) > 0) { long long next = more & mask[step]; if (size + bits[next & MASK] + bits[(next >> 21) & MASK] + bits[next >> 42] >= this->size && size + cmax[step] > this->size) { // the current node is in the clique if (search(step + 1, size + 1, next, cons | now)) return true; } } long long next = more & ~now; if (size + bits[next & MASK] + bits[(next >> 21) & MASK] + bits[next >> 42] > this->size) { // the current node is not in the clique if (search(step + 1, size, next, cons)) return true; } return false; } // solve maximum clique and return size int clique::sizeClique(vector >& mat) { n = mat.size(); // generate mask vectors for (int i = 0; i < n; ++i) { mask[i] = 0; for (int j = 0; j < n; ++j) if (mat[i][j] > 0) mask[i] |= ONE << j; } size = 0; for (int i = n - 1; i >= 0; --i) { search(i + 1, 1, mask[i], ONE << i); cmax[i] = size; } return size; } // solve maximum clique and return constitution // calls sizeClique and restore cons vector clique::consClique(vector >& mat) { sizeClique(mat); vector ret; for (int i = 0; i < n; ++i) if ((cons & (ONE << i)) > 0) ret.push_back(i); return ret; } 5、 图论—连通性 7.1 无向图关键点(dfs邻接阵) //无向图的关键点,dfs邻接阵形式,O(n^2) //返回关键点个数,key[]返回点集 //传入图的大小n和邻接阵mat,不相邻点边权0 #define MAXN 110 void search(int n,int mat[][MAXN],int* dfn,int* low,int now,int& ret,int* key,int& cnt,int root,int& rd,int* bb){ int i; dfn[now]=low[now]=++cnt; for (i=0;i=dfn[now]){ if (now!=root&&!bb[now]) key[ret++]=now,bb[now]=1; else if(now==root) rd++; } } else if (dfn[i]1&&!bb[i]) key[ret++]=i,bb[i]=1; } return ret; } 7.2 无向图关键边(dfs邻接阵) //无向图的关键边,dfs邻接阵形式,O(n^2) //返回关键边条数,key[][2]返回边集 //传入图的大小n和邻接阵mat,不相邻点边权0 #define MAXN 100 void search(int n,int mat[][MAXN],int* dfn,int* low,int now,int& cnt,int key[][2]){ int i; for (low[now]=dfn[now],i=0;idfn[now]) key[cnt][0]=i,key[cnt++][1]=now; if (low[i] void dummy(int n,int* a){ for (int i=0;i=dfn[now]){ for (st[sp]=-1,a[0]=now,m=1;st[sp]!=i;a[m++]=st[--sp]); dummy(m,a); } } else if (dfn[i]next) if (!v[t=e->to]){ if (match[t]<0) match[now]=t,match[t]=now,ret=1; else{ v[t]=1; if (aug(n,list,match,v,match[t])) match[now]=t,match[t]=now,ret=1; v[t]=0; } if (ret) break; } v[now]=0; return ret; } int graph_match(int n,edge_t* list[],int* match){ int v[MAXN],i,j; for (i=0;i=2;) if (match[i]<0&&aug(n,list,match,v,i)) i=0,j-=2; else i++; for (i=j=0;i=0); return j/2; } 7、 图论—网络流 9.1 最大流(邻接阵) //求网络最大流,邻接阵形式 //返回最大流量,flow返回每条边的流量 //传入网络节点数n,容量mat,源点source,汇点sink #define MAXN 100 #define inf 1000000000 int max_flow(int n,int mat[][MAXN],int source,int sink,int flow[][MAXN]){ int pre[MAXN],que[MAXN],d[MAXN],p,q,t,i,j; if (source==sink) return inf; for (i=0;i0) flow[pre[i]-1][i]+=d[sink],i=pre[i]-1; else flow[i][-pre[i]-1]-=d[sink],i=-pre[i]-1; } for (j=i=0;ic[j])) j=i; if (j<0) return ret; if (j==sink) break; for (v[j]=1,i=0;ic[i]&&c[j]>c[i]) c[i]=mat[j][i]0) flow[pre[i]-1][i]+=d[sink],i=pre[i]-1; else flow[i][-pre[i]-1]-=d[sink],i=-pre[i]-1; } for (j=i=0;i=0;i--) while (mat[now][i]){ mat[now][i]--,mat[i][now]--; find_path_u(n,mat,i,step,path); } path[step++]=now; } void find_path_d(int n,int mat[][MAXN],int now,int& step,int* path){ int i; for (i=n-1;i>=0;i--) while (mat[now][i]){ mat[now][i]--; find_path_d(n,mat,i,step,path); } path[step++]=now; } int euclid_path(int n,int mat[][MAXN],int start,int* path){ int ret=0; find_path_u(n,mat,start,ret,path); // find_path_d(n,mat,start,ret,path); return ret; } 10.2 树的前序表转化 //将用边表示的树转化为前序表示的树 //传入节点数n和邻接表list[],邻接表必须是双向的,会在函数中释放 //pre[]返回前序表,map[]返回前序表中的节点到原来节点的映射 #define MAXN 10000 struct node{ int to; node* next; }; void prenode(int n,node* list[],int* pre,int* map,int* v,int now,int last,int& id){ node* t; int p=id++; for (v[map[p]=now]=1,pre[p]=last;list[now];){ t=list[now],list[now]=t->next; if (!v[t->to]) prenode(n,list,pre,map,v,t->to,p,id); } } void makepre(int n,node* list[],int* pre,int* map){ int v[MAXN],id=0,i; for (i=0;i=0;i--) if (!c[i]){ set[i]=1; if (pre[i]!=-1) c[pre[i]]=1; ret++; } return ret; } //最大边独立集 int max_edge_independent(int n,int* pre,int* set){ int c[MAXN],i,ret=0; for (i=0;i=0;i--) if (!c[i]&&pre[i]!=-1&&!c[pre[i]]){ set[i]=1; c[pre[i]]=1; ret++; } return ret; } //最小顶点覆盖集 int min_node_cover(int n,int* pre,int* set){ int c[MAXN],i,ret=0; for (i=0;i=0;i--) if (!c[i]&&pre[i]!=-1&&!c[pre[i]]){ set[i]=1; c[pre[i]]=1; ret++; } return ret; } //最小顶点支配集 int min_node_dominant(int n,int* pre,int* set){ int c[MAXN],i,ret=0; for (i=0;i=0;i--) if (!c[i]&&(pre[i]==-1||!set[pre[i]])){ if (pre[i]!=-1){ set[pre[i]]=1; c[pre[i]]=1; if (pre[pre[i]]!=-1) c[pre[pre[i]]]=1; } else set[i]=1; ret++; } return ret; } 10.4 拓扑排序(邻接阵) //拓扑排序,邻接阵形式,复杂度O(n^2) //如果无法完成排序,返回0,否则返回1,ret返回有序点列 //传入图的大小n和邻接阵mat,不相邻点边权0 #define MAXN 100 int toposort(int n,int mat[][MAXN],int* ret){ int d[MAXN],i,j,k; for (i=0;ic[j])) j=i; if (j<0) return ret; if (j==sink) break; for (v[j]=1,i=0;ic[i]&&c[j]>c[i]) c[i]=mat[j][i]c[j])) j=i; if (j<0) return ret; if (j==sink) break; for (v[j]=1,i=0;ic[i]&&c[j]>c[i]) c[i]=mat[j][i]c[j])) j=i; if (j<0) return ret; if (j==sink) break; for (v[j]=1,i=0;ic[i]&&c[j]>c[i]) c[i]=mat[j][i]c[j])) j=i; if (j<0) return ret; if (j==sink) break; for (v[j]=1,i=0;ic[i]&&c[j]>c[i]) c[i]=mat[j][i] #define MAXN 310 #define _clr(x) memset(x,0xff,sizeof(int)*n) int hungary(int n,int mat[][MAXN],int* match1,int* match2){ int s[MAXN],t[MAXN],p,q,ret=0,i,j,k; for (_clr(match1),_clr(match2),i=0;i=0)) for (_clr(t),s[p=q=0]=i;p<=q&&match1[i]<0;p++) for (k=s[p],j=0;j=0;j=p) match2[j]=k=t[j],p=match1[k],match1[k]=j; } return ret; } inline int path_cover(int n,int mat[][MAXN],int* pre,int* next){ return n-hungary(n,mat,next,pre); } 9、 应用 13.1 Joseph问题 // Joseph's Problem // input: n,m -- the number of persons, the inteval between persons // output: -- return the reference of last person int josephus0(int n, int m) { if (n == 2) return (m%2) ? 2 : 1; int v = (m+josephus0(n-1,m)) % n; if (v == 0) v = n; return v; } int josephus(int n, int m) { if (m == 1) return n; if (n == 1) return 1; if (m >=n) return josephus0(n,m); int l = (n/m)*m; int j = josephus(n - (n/m), m); if (j <= n-l) return l+j; j -= n-l; int t = (j/(m-1))*m; if ((j % (m-1)) == 0) return t-1; return t + (j % (m-1)); } 13.2 N皇后构造解 //N皇后构造解,n>=4 void even1(int n,int *p){ int i; for (i=1;i<=n/2;i++) p[i-1]=2*i; for (i=n/2+1;i<=n;i++) p[i-1]=2*i-n-1; } void even2(int n,int *p){ int i; for (i=1;i<=n/2;i++) p[i-1]=(2*i+n/2-3)%n+1; for (i=n/2+1;i<=n;i++) p[i-1]=n-(2*(n-i+1)+n/2-3)%n; } void generate(int,int*); void odd(int n,int *p){ generate(n-1,p),p[n-1]=n; } void generate(int n,int *p){ if (n&1) odd(n,p); else if (n%6!=2) even1(n,p); else even2(n,p); } 13.3 布尔母函数 //布尔母函数 //判m[]个价值为w[]的货币能否构成value //适合m[]较大w[]较小的情况 //返回布尔量 //传入货币种数n,个数m[],价值w[]和目标值value #define MAXV 100000 int genfunc(int n,int* m,int* w,int value){ int i,j,k,c; char r[MAXV]; for (r[0]=i=1;i<=value;r[i++]=0); for (i=0;i>1];ij) l=j+1; else r=j; } return a[k]; } 13.5 幻方构造 //幻方构造(l!=2) #define MAXN 100 void dllb(int l,int si,int sj,int sn,int d[][MAXN]){ int n,i=0,j=l/2; for (n=1;n<=l*l;n++){ d[i+si][j+sj]=n+sn; if (n%l){ i=(i)?(i-1):(l-1); j=(j==l-1)?0:(j+1); } else i=(i==l-1)?0:(i+1); } } void magic_odd(int l,int d[][MAXN]){ dllb(l,0,0,0,d); } void magic_4k(int l,int d[][MAXN]){ int i,j; for (i=0;i=0&&!_match(pat[i+1],pat[j]);i=fail[i]); fail[j]=(_match(pat[i+1],pat[j])?i+1:-1); } for (i=j=0;i #define MAXN 1000000 #define _cp(a,b) ((a)<=(b)) typedef int elem_t; elem_t _tmp[MAXN]; int inv(int n,elem_t* a){ int l=n>>1,r=n-l,i,j; int ret=(r>1?(inv(l,a)+inv(r,a+l)):0); for (i=j=0;i<=l;_tmp[i+j]=a[i],i++) for (ret+=j;j int MinString(vector &str) { int i, j, k; vector ss(str.size() << 1); for (i = 0; i < str.size(); i ++) ss[i] = ss[i + str.size()] = str[i]; for (i = k = 0, j = 1; k < str.size() && i < str.size() && j < str.size(); ) { for (k = 0; k < str.size() && ss[i + k] == ss[j + k]; k ++); if (k < str.size()) { if (ss[i + k] > ss[j + k]) i += k + 1; else j += k + 1; if (i == j) j ++; } } return i < j ? i : j; } 13.9 最长公共单调子序列 // 最长公共递增子序列, 时间复杂度O(n^2 * logn),空间 O(n^2) /** * n为a的大小, m为b的大小 * 结果在ans中 * "define _cp(a,b) ((a)<(b))"求解最长严格递增序列 */ #define MAXN 1000 #define _cp(a,b) ((a)<(b)) typedef int elem_t; elem_t DP[MAXN][MAXN]; int num[MAXN], p[1<<20]; int LIS(int n, elem_t *a, int m, elem_t *b, elem_t *ans){ int i, j, l, r, k; DP[0][0] = 0; num[0] = (b[0] == a[0]); for(i = 1; i < m; i++) { num[i] = (b[i] == a[0]) || num[i-1]; DP[i][0] = 0; } for(i = 1; i < n; i++){ if(b[0] == a[i] && !num[0]) { num[0] = 1; DP[0][0] = i<<10; } for(j = 1; j < m; j++){ for(k=((l=0)+(r=num[j-1]-1))>>1; l<=r; k=(l+r)>>1) if(_cp(a[DP[j-1][k]>>10], a[i])) l=k+1; else r=k-1; if(l < num[j-1] && i == (DP[j-1][l]>>10) ){ if(l >= num[j]) DP[j][num[j]++] = DP[j-1][l]; else DP[j][l] = _cp(a[DP[j][l]>>10],a[i]) ? DP[j][l] : DP[j-1][l]; } if(b[j] == a[i]){ for(k=((l=0)+(r=num[j]-1))>>1; l<=r; k=(l+r)>>1) if(_cp(a[DP[j][k]>>10], a[i])) l=k+1; else r=k-1; DP[j][l] = (i<<10) + j; num[j] += (l>=num[j]); p[DP[j][l]] = l ? DP[j][l-1] : -1; } } } for (k=DP[m-1][i=num[m-1]-1];i>=0;ans[i--]=a[k>>10],k=p[k]); return num[m-1]; } 13.10 最长子序列 //最长单调子序列,复杂度O(nlogn) //注意最小序列覆盖和最长序列的对应关系,例如 //"define _cp(a,b) ((a)>(b))"求解最长严格递减序列,则 //"define _cp(a,b) (!((a)>(b)))"求解最小严格递减序列覆盖 //可更改元素类型和比较函数 #define MAXN 10000 #define _cp(a,b) ((a)>(b)) typedef int elem_t; int subseq(int n,elem_t* a){ int b[MAXN],i,l,r,m,ret=0; for (i=0;iret)) for (m=((l=1)+(r=ret))>>1;l<=r;m=(l+r)>>1) if (_cp(a[b[m]],a[i])) l=m+1; else r=m-1; return ret; } int subseq(int n,elem_t* a,elem_t* ans){ int b[MAXN],p[MAXN],i,l,r,m,ret=0; for (i=0;iret)) for (m=((l=1)+(r=ret))>>1;l<=r;m=(l+r)>>1) if (_cp(a[b[m]],a[i])) l=m+1; else r=m-1; for (m=b[i=ret];i;ans[--i]=a[m],m=p[m]); return ret; } 13.11 最大子串匹配 //最大子串匹配,复杂度O(mn) //返回最大匹配值,传入两个串和串的长度,重载返回一个最大匹配 //注意做字符串匹配是串末的'\0'没有置! //可更改元素类型,更换匹配函数和匹配价值函数 #include #define MAXN 100 #define max(a,b) ((a)>(b)?(a):(b)) #define _match(a,b) ((a)==(b)) #define _value(a,b) 1 typedef char elem_t; int str_match(int m,elem_t* a,int n,elem_t* b){ int match[MAXN+1][MAXN+1],i,j; memset(match,0,sizeof(match)); for (i=0;imatch[i+1][j]?match[i][j+1]:match[i+1][j]); last[i+1][j+1]=(match[i][j+1]>match[i+1][j]?3:1); if ((t=(match[i][j]+_value(a[i],b[i]))*_match(a[i],b[j]))>match[i+1][j+1]) match[i+1][j+1]=t,last[i+1][j+1]=2; } for (;match[i][j];i-=(last[t=i][j]>1),j-=(last[t][j]<3)) ret[match[i][j]-1]=(last[i][j]<3?a[i-1]:b[j-1]); return match[m][n]; } 13.12 最大子段和 //求最大子段和,复杂度O(n) //传入串长n和内容list[] //返回最大子段和,重载返回子段位置(maxsum=list[start]+...+list[end]) //可更改元素类型 typedef int elem_t; elem_t maxsum(int n,elem_t* list){ elem_t ret,sum=0; int i; for (ret=list[i=0];i0?sum:0)+list[i],ret=(sum>ret?sum:ret); return ret; } elem_t maxsum(int n,elem_t* list,int& start,int& end){ elem_t ret,sum=0; int s,i; for (ret=list[start=end=s=i=0];i0?s:i)) if ((sum=(sum>0?sum:0)+list[i])>ret) ret=sum,start=s,end=i; return ret; } 13.13 最大子阵和 //求最大子阵和,复杂度O(n^3) //传入阵的大小m,n和内容mat[][] //返回最大子阵和,重载返回子阵位置(maxsum=list[s1][s2]+...+list[e1][e2]) //可更改元素类型 #define MAXN 100 typedef int elem_t; elem_t maxsum(int m,int n,elem_t mat[][MAXN]){ elem_t matsum[MAXN][MAXN+1],ret,sum; int i,j,k; for (i=0;i0?sum:0)+matsum[i][k+1]-matsum[i][j],ret=(sum>ret?sum:ret); return ret; } elem_t maxsum(int m,int n,elem_t mat[][MAXN],int& s1,int& s2,int& e1,int& e2){ elem_t matsum[MAXN][MAXN+1],ret,sum; int i,j,k,s; for (i=0;i0?s:i)) if ((sum=(sum>0?sum:0)+matsum[i][k+1]-matsum[i][j])>ret) ret=sum,s1=s,s2=i,e1=j,e2=k; return ret; } 14.4 线性方程组 #define MAXN 100 #define fabs(x) ((x)>0?(x):-(x)) #define eps 1e-10 //列主元gauss消去求解a[][]x[]=b[] //返回是否有唯一解,若有解在b[]中 int gauss_cpivot(int n,double a[][MAXN],double b[]){ int i,j,k,row; double maxp,t; for (k=0;kfabs(maxp)) maxp=a[row=i][k]; if (fabs(maxp)=0;i--) for (j=i+1;jfabs(maxp)) maxp=a[row=i][col=j]; if (fabs(maxp)=0;i--) for (j=i+1;j #define MAXN 100 #define eps 1e-10 int linear_dependent(int m,int n,double vec[][MAXN]){ double ort[MAXN][MAXN],e; int i,j,k; if (m>n) return 1; for (i=0;i12) return 0; if (a.month==2) return a.day>0&&a.day<=28+leap(a.year); return a.day>0&&a.day<=days[a.month-1]; } //比较日期大小 inline int datecmp(date a,date b){ if (a.year!=b.year) return a.year-b.year; if (a.month!=b.month) return a.month-b.month; return a.day-b.day; } //返回指定日期是星期几 int weekday(date a){ int tm=a.month>=3?(a.month-2):(a.month+10); int ty=a.month>=3?a.year:(a.year-1); return (ty+ty/4-ty/100+ty/400+(int)(2.6*tm-0.2)+a.day)%7; } //日期转天数偏移 int date2int(date a){ int ret=a.year*365+(a.year-1)/4-(a.year-1)/100+(a.year-1)/400,i; days[1]+=leap(a.year); for (i=0;i=365+leap(ret.year);a-=365+leap(ret.year),ret.year++); days[1]+=leap(ret.year); for (ret.month=1;a>=days[ret.month-1];a-=days[ret.month-1],ret.month++); days[1]=28; ret.day=a+1; return ret; }

下载文档到电脑,查找使用更方便

文档的实际排版效果,会与网站的显示效果略有不同!!

需要 8 金币 [ 分享文档获得金币 ] 0 人已下载

下载文档