I Zhejiang University ICPC Team Routine Library by WishingBone (Dec. 2002) Last Update (Nov. 2004) by Riveria II 1、 几何...........................................................................................................................1 1.1 注意...................................................................................................................1 1.2 几何公式...........................................................................................................1 1.3 多边形...............................................................................................................3 1.4 多边形切割.......................................................................................................6 1.5 浮点函数...........................................................................................................7 1.6 面积.................................................................................................................12 1.7 球面.................................................................................................................13 1.8 三角形.............................................................................................................14 1.9 三维几何.........................................................................................................16 1.10 凸包.................................................................................................................23 1.11 网格.................................................................................................................25 1.12 圆.....................................................................................................................25 1.13 整数函数.........................................................................................................27 2、 组合.................................................................................................................30 2.1 组合公式..................................................................................................................30 2.2 排列组合生成..........................................................................................................30 2.3 生成 gray 码.............................................................................................................32 2.4 置换(polya) ..............................................................................................................32 2.5 字典序全排列..........................................................................................................33 2.6 字典序组合..............................................................................................................33 3、............... 结构 34 3.1 并查集......................................................................................................................34 3.2 堆..............................................................................................................................35 3.3 线段树......................................................................................................................36 3.4 子段和......................................................................................................................41 3.5 子阵和......................................................................................................................41 4、.. 数论 42 4.1 阶乘最后非 0 位......................................................................................................42 4.2 模线性方程组..........................................................................................................43 4.3 素数..........................................................................................................................44 4.4 欧拉函数..................................................................................................................45 III 5、 数值计算.................................................................................................................46 5.1 定积分计算(Romberg) ............................................................................................46 5.2 多项式求根(牛顿法)...............................................................................................48 5.3 周期性方程(追赶法)...............................................................................................49 6、......... 图论—NP 搜索 50 6.1 最大团......................................................................................................................50 6.2 最大团(n<64)(faster) ...............................................................................................51 7、 图论—连通性 53 7.1 无向图关键点(dfs 邻接阵) .....................................................................................53 7.2 无向图关键边(dfs 邻接阵) .....................................................................................54 7.3 无向图的块(bfs 邻接阵) .........................................................................................55 7.4 无向图连通分支(dfs/bfs 邻接阵) ...........................................................................56 7.5 有向图强连通分支(dfs/bfs 邻接阵) .......................................................................57 7.6 有向图最小点基(邻接阵) .......................................................................................58 8、....... 图论—匹配 59 8.1 二分图最大匹配(hungary 邻接表).........................................................................59 8.2 二分图最大匹配(hungary 邻接阵).........................................................................60 8.3 二分图最大匹配(hungary 正向表).........................................................................60 8.4 二分图最佳匹配(kuhn_munkras 邻接阵) ...............................................................61 8.5 一般图匹配(邻接表)...............................................................................................62 8.6 一般图匹配(邻接阵)...............................................................................................63 8.7 一般图匹配(正向表)...............................................................................................63 9、....... 图论—网络流 64 9.1 最大流(邻接阵).......................................................................................................64 9.2 上下界最大流(邻接阵)...........................................................................................65 9.3 上下界最小流(邻接阵)...........................................................................................66 9.4 最大流无流量(邻接阵)...........................................................................................67 9.5 最小费用最大流(邻接阵) .......................................................................................67 10 、....... 图论—应用 68 10.1 欧拉回路(邻接阵).................................................................................................68 10.2 树的前序表转化....................................................................................................69 10.3 树的优化算法........................................................................................................70 10.4 拓扑排序(邻接阵).................................................................................................71 10.5 最佳边割集............................................................................................................72 10.6 最佳点割集............................................................................................................73 10.7 最小边割集............................................................................................................74 10.8 最小点割集............................................................................................................75 10.9 最小路径覆盖........................................................................................................77 11 、......... 图论—支撑树 77 11.1 最小生成树(kruskal 邻接表).................................................................................77 11.2 最小生成树(kruskal 正向表).................................................................................79 11.3 最小生成树(prim+binary_heap 邻接表)...............................................................80 11.4 最小生成树(prim+binary_heap 正向表)...............................................................81 11.5 最小生成树(prim+mapped_heap 邻接表) ............................................................82 IV 11.6 最小生成树(prim+mapped_heap 正向表) ............................................................84 11.7 最小生成树(prim 邻接阵).....................................................................................85 11.8 最小树形图(邻接阵) .............................................................................................85 12 、......... 图论—最短路径 87 12.1 最短路径(单源 bellman_ford 邻接阵)..................................................................87 12.2 最短路径(单源 dijkstra+bfs 邻接表)....................................................................87 12.3 最短路径(单源 dijkstra+bfs 正向表)....................................................................88 12.4 最短路径(单源 dijkstra+binary_heap 邻接表) .....................................................89 12.5 最短路径(单源 dijkstra+binary_heap 正向表) .....................................................90 12.6 最短路径(单源 dijkstra+mapped_heap 邻接表)...................................................91 12.7 最短路径(单源 dijkstra+mapped_heap 正向表)...................................................92 12.8 最短路径(单源 dijkstra 邻接阵) ...........................................................................93 12.9 最短路径(多源 floyd_warshall 邻接阵) ...............................................................94 13 、......... 应用 94 13.1 Joseph 问题.............................................................................................................94 13.2 N 皇后构造解.........................................................................................................95 13.3 布尔母函数............................................................................................................96 13.4 第 k 元素................................................................................................................96 13.5 幻方构造................................................................................................................97 13.6 模式匹配(kmp)......................................................................................................98 13.7 逆序对数................................................................................................................99 13.8 字符串最小表示....................................................................................................99 13.9 最长公共单调子序列..........................................................................................100 13.10 最长子序列........................................................................................................101 13.11 最大子串匹配....................................................................................................102 13.12 最大子段和........................................................................................................103 13.13 最大子阵和........................................................................................................103 14 、 .............. 其它 104 14.1 大数(只能处理正数)...........................................................................................104 14.2 分数......................................................................................................................110 14.3 矩阵......................................................................................................................112 14.4 线性方程组..........................................................................................................114 14.5 线性相关..............................................................................................................116 14.6 日期......................................................................................................................116 1 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) 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 3 圆柱: 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); 4 } // 判定凸多边形,顶点按顺时针或逆时针给出,允许相邻边共线 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 多边形切割多边形切割多边形切割多边形切割 // 多边形切割 // 可用于半平面交 #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 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(int& n,point* p,point l1,point l2,point side){ point pp[100]; int m=0,i; for (i=0;i #define eps 1e-8 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)); 10 } 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,u 2); } // 判两线段相交,不包括端点和部分重合 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; 11 } // 点到直线上的最近点 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; } 13 // 计算三角形面积,输入三边长 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))); } 14 // 计算球面距离,r 为球半径 inline double sphere_dist(double r,double lng1,double lat1,double lng2,double lat2){ return r*angle(lng1,lat1,lng2,lat2); } 1.8 三角形三角形三角形三角形 #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); 15 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; 16 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.9 三维几何三维几何三维几何三维几何 // 三维几何函数库 #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 19 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); } 26 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); 28 } 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 29 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,u 2); } // 判两线段相交,不包括端点和部分重合 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); } 30 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++); } 3、、、、 结构结构结构结构 3.1 并查集并查集并查集并查集 // 带路径压缩的并查集,用于动态维护查询等价类 // 图论算法中动态判点集连通常用 // 维护和查询复杂度略大于 O(1) // 集合元素取值 1..MAXN-1( 注意 0 不能用!), 默认不等价 #include #define MAXN 100000 #define _ufind_run(x) for(;p[t=x];x=p[x],p[t]=(p[x]?p[x]:x)) #define _run_both _ufind_run(i);_ufind_run(j) struct ufind{ int p[MAXN],t; void init(){memset(p,0,sizeof(p));} void set_friend(int i,int j){_run_both;p[i]=(i==j?0:j);} int is_friend(int i,int j){_run_both;return i==j&&i;} }; // 带路径压缩的并查集扩展形式 // 用于动态维护查询 friend-enemy 型等价类 // 维护和查询复杂度略大于 O(1) // 集合元素取值 1..MAXN-1( 注意 0 不能用!), 默认无关 #include #define MAXN 100000 #define sig(x) ((x)>0?1:-1) 35 #define abs(x) ((x)>0?(x):-(x)) #define _ufind_run(x) for(;p[t=abs(x)];x=sig(x)*p[abs(x)],p[t]=sig(p[t])*(p[abs(x)]?p[abs(x)]:abs(p[t]))) #define _run_both _ufind_run(i);_ufind_run(j) #define _set_side(x) p[abs(i)]=sig(i)*(abs(i)==abs(j)?0:(x)*j) #define _judge_side(x) (i==(x)*j&&i) struct ufind{ int p[MAXN],t; void init(){memset(p,0,sizeof(p));} int set_friend(int i,int j){_run_both;_set_side(1);return !_judge_side(-1);} int set_enemy(int i,int j){_run_both;_set_side(-1);return !_judge_side(1);} int is_friend(int i,int j){_run_both;return _judge_side(1);} int is_enemy(int i,int j){_run_both;return _judge_side(-1);} }; 3.2 堆堆堆堆 // 二分堆(binary) // 可插入,获取并删除最小(最大)元素,复杂度均 O(logn) // 可更改元素类型,修改比较符号或换成比较函数 #define MAXN 10000 #define _cp(a,b) ((a)<(b)) typedef int elem_t; struct heap{ elem_t h[MAXN]; int n,p,c; void init(){n=0;} void ins(elem_t e){ for (p=++n;p>1&&_cp(e,h[p>>1]);h[p]=h[p>>1],p>>=1); h[p]=e; } int del(elem_t& e){ if (!n) return 0; for (e=h[p=1],c=2;c1&&_cp(e,h[p>>1]);h[map[ind[p]=ind[p>>1]]=p]=h[p>>1],p>>=1); h[map[ind[p]=i]=p]=e; } int del(int i,elem_t& e){ i=map[i];if (i<1||i>n) return 0; for (e=h[p=i];p>1;h[map[ind[p]=ind[p>>1]]=p]=h[p>>1],p>>=1); for (c=2;c>1; if (lm0) inc_seg(t+t+1,m0,r0,m0>l?m0:l,r); if (cnt[t+t]&&cnt[t+t+1]){ cnt[t+t]--; update(t+t,l0,m0); cnt[t+t+1]--; update(t+t+1,m0,r0); cnt[t]++; } } update(t,l0,r0); } void segtree::dec_seg(int t,int l0,int r0,int l,int r){ if (l0==l&&r0==r) cnt[t]--; else if (cnt[t]){ cnt[t]--; if (l>l0) inc_seg(t,l0,r0,l0,l); if (r>1; if (lm0) dec_seg(t+t+1,m0,r0,m0>l?m0:l,r); } update(t,l0,r0); } int segtree::seg_len(int t,int l0,int r0,int l,int r){ if (cnt[t]||(l0==l&&r0==r)) return len[t]; else{ int m0=(l0+r0)>>1,ret=0; if (lm0) ret+=seg_len(t+t+1,m0,r0,m0>l?m0:l,r); return ret; 39 } } // 线段树扩展 // 可以计算长度和线段数 // 可以处理加入边和删除边不同的情况 //inc_seg 和 dec_seg 用于加入边 //seg_len 求长度,seg_cut 求线段数 //t 传根节点(一律为 1) //l0,r0 传树的节点范围(一律为 1..t) //l,r 传线段(端点) #define MAXN 10000 struct segtree{ int n,cnt[MAXN],len[MAXN],cut[MAXN],bl[MAXN],br[MAXN]; segtree(int t):n(t){ for (int i=1;i<=t;i++) cnt[i]=len[i]=cut[i]=bl[i]=br[i]=0; }; void update(int t,int l,int r); void inc_seg(int t,int l0,int r0,int l,int r); void dec_seg(int t,int l0,int r0,int l,int r); int seg_len(int t,int l0,int r0,int l,int r); int seg_cut(int t,int l0,int r0,int l,int r); }; int length(int l,int r){ return r-l; } void segtree::update(int t,int l,int r){ if (cnt[t]||r-l==1) len[t]=length(l,r),cut[t]=bl[t]=br[t]=1; else{ len[t]=len[t+t]+len[t+t+1]; cut[t]=cut[t+t]+cut[t+t+1]; if (br[t+t]&&bl[t+t+1]) cut[t]--; bl[t]=bl[t+t],br[t]=br[t+t+1]; } } void segtree::inc_seg(int t,int l0,int r0,int l,int r){ if (l0==l&&r0==r) cnt[t]++; 40 else{ int m0=(l0+r0)>>1; if (lm0) inc_seg(t+t+1,m0,r0,m0>l?m0:l,r); if (cnt[t+t]&&cnt[t+t+1]){ cnt[t+t]--; update(t+t,l0,m0); cnt[t+t+1]--; update(t+t+1,m0,r0); cnt[t]++; } } update(t,l0,r0); } void segtree::dec_seg(int t,int l0,int r0,int l,int r){ if (l0==l&&r0==r) cnt[t]--; else if (cnt[t]){ cnt[t]--; if (l>l0) inc_seg(t,l0,r0,l0,l); if (r>1; if (lm0) dec_seg(t+t+1,m0,r0,m0>l?m0:l,r); } update(t,l0,r0); } int segtree::seg_len(int t,int l0,int r0,int l,int r){ if (cnt[t]||(l0==l&&r0==r)) return len[t]; else{ int m0=(l0+r0)>>1,ret=0; if (lm0) ret+=seg_len(t+t+1,m0,r0,m0>l?m0:l,r); return ret; } } int segtree::seg_cut(int t,int l0,int r0,int l,int r){ if (cnt[t]) return 1; if (l0==l&&r0==r) return cut[t]; else{ int m0=(l0+r0)>>1,ret=0; if (lm0) ret+=seg_cut(t+t+1,m0,r0,m0>l?m0:l,r); if (lm0&&br[t+t]&&bl[t+t+1]) ret--; return ret; } } 3.4 子段和子段和子段和子段和 // 求 sum{[0..n-1]} // 维护和查询复杂度均为 O(logn) // 用于动态求子段和,数组内容保存在 sum.a[] 中 // 可以改成其他数据类型 #include #define lowbit(x) ((x)&((x)^((x)-1))) #define MAXN 10000 typedef int elem_t; struct sum{ elem_t a[MAXN],c[MAXN],ret; int n; void init(int i){memset(a,0,sizeof(a));memset(c,0,sizeof(c));n=i;} void update(int i,elem_t v){for (v-=a[i],a[i++]+=v;i<=n;c[i-1]+=v,i+=lowbit(i));} elem_t query(int i){for (ret=0;i;ret+=c[i-1],i^=lowbit(i));return ret;} }; 3.5 子阵和子阵和子阵和子阵和 42 // 求 sum{a[0..m-1][0..n-1]} // 维护和查询复杂度均为 O(logm*logn) // 用于动态求子阵和,数组内容保存在 sum.a[][] 中 // 可以改成其他数据类型 #include #define lowbit(x) ((x)&((x)^((x)-1))) #define MAXN 100 typedef int elem_t; struct sum{ elem_t a[MAXN][MAXN],c[MAXN][MAXN],ret; int m,n,t; void init(int i,int j){memset(a,0,sizeof(a));memset(c,0,sizeof(c));m=i,n=j;} void update(int i,int j,elem_t v){ for (v-=a[i][j],a[i++][j++]+=v,t=j;i<=m;i+=lowbit(i)) for (j=t;j<=n;c[i-1][j-1]+=v,j+=lowbit(j)); } elem_t query(int i,int j){ for (ret=0,t=j;i;i^=lowbit(i)) for (j=t;j;ret+=c[i-1][j-1],j^=lowbit(j)); return ret; } }; 4、、、、 数论数论数论数论 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--) 43 c=c*10+a[i],a[i]=c/5,c%=5; } return ret+ret%2*5; } 4.2 模线性方程组模线性方程组模线性方程组模线性方程组 #ifdef WIN32 typedef __int64 i64; #else typedef long long i64; #endif // 扩展 Euclid 求解 gcd(a,b)=ax+by int ext_gcd(int a,int b,int& x,int& y){ int t,ret; if (!b){ x=1,y=0; return a; } ret=ext_gcd(b,a%b,x,y); t=x,x=y,y=t-a/b*y; return ret; } // 计算 m^a, O(loga), 本身没什么用, 注意这个按位处理的方法 :-P int exponent(int m,int a){ int ret=1; for (;a;a>>=1,m*=m) if (a&1) ret*=m; return ret; } // 计算幂取模 a^b mod n, O(logb) int modular_exponent(int a,int b,int n){ //a^b mod n int ret=1; for (;b;b>>=1,a=(int)((i64)a)*a%n) if (b&1) ret=(int)((i64)ret)*a%n; return ret; } // 求解模线性方程 ax=b (mod n) // 返回解的个数,解保存在 sol[] 中 44 // 要求 n>0, 解的范围 0..n-1 int modular_linear(int a,int b,int n,int* sol){ int d,e,x,y,i; d=ext_gcd(a,n,x,y); if (b%d) return 0; e=(x*(b/d)%n+n)%n; for (i=0;i0,w[i] 与 w[j] 互质,解的范围 1..n,n=w[0]*w[1]*...*w[k-1] int modular_linear_system(int b[],int w[],int k){ int d,x,y,a=0,m,n=1,i; for (i=0;i1; } 45 void initprime(){ int i; for (plist[pcount++]=2,i=3;i<50000;i++) if (prime(i)) plist[pcount++]=i; } //miller rabin // 判断自然数 n 是否为素数 //time 越高失败概率越低,一般取 10 到 50 #include #ifdef WIN32 typedef __int64 i64; #else typedef long long i64; #endif int modular_exponent(int a,int b,int n){ //a^b mod n int ret; for (;b;b>>=1,a=(int)((i64)a)*a%n) if (b&1) ret=(int)((i64)ret)*a%n; return ret; } // Carmicheal number: 561,41041,825265,321197185 int miller_rabin(int n,int time=10){ if (n==1||(n!=2&&!(n%2))||(n!=3&&!(n%3))||(n!=5&&!(n%5))||(n!=7&&!(n%7))) return 0; while (time--) if (modular_exponent(((rand()&0x7fff<<16)+rand()&0x7fff+rand()&0x7fff)%(n-1)+1,n-1,n)!=1) return 0; return 1; } 4.4 欧拉函数欧拉函数欧拉函数欧拉函数 int gcd(int a,int b){ return b?gcd(b,a%b):a; } inline int lcm(int a,int b){ 46 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; } 5、、、、 数值计算数值计算数值计算数值计算 5.1 定积分计算定积分计算定积分计算定积分计算(Romberg) /* Romberg 求定积分 输入：积分区间[a,b] ，被积函数 f(x,y,z) 输出：积分结果 f(x,y,z) 示例： double f0( double x, double l, double t ) { return sqrt(1.0+l*l*t*t*cos(t*x)*cos(t*x)); } */ double Integral(double a, double b, double (*f)(double x, double y, double z), double eps, double l, double t) double Romberg (double a, double b, double (*f)(double x, double y, double z), double eps, double l, double t) { #define MAX_N 1000 int i, j, temp2, min; double h, R[2][MAX_N], temp4; 47 for (i=0; imin)) return R[1][i-1]; h *= 0.50; temp2 *= 2; for (j=0; j0; 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) 49 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)= 0; i --) { a[i] = -a[i] - c[i] * a[i + 1]; x[i] -= c[i] * x[i + 1]; } x[N - 1] -= (c[N - 1] * x[0] + a[N - 1] * x[N - 2]); x[N - 1] /= (c[N - 1] * a[0] + a[N - 1] * a[N - 2] + b[N - 1]); for (int i = N - 2; i >= 0; i --) x[i] += a[i] * x[N - 1]; } 6、、、、 图论图论图论图论————NP 搜索搜索搜索搜索 6.1 最大团最大团最大团最大团 // 最大团 // 返回最大团大小和一个方案,传入图的大小 n 和邻接阵 mat //mat[i][j] 为布尔量 #define MAXN 60 void clique(int n, int* u, int mat[][MAXN], int size, int& max, int& bb, int* res, int* rr, int* c) { int i, j, vn, v[MAXN]; if (n) { if (size + c[u[0]] <= max) return; for (i = 0; i < n + size - max && i < n; ++ i) { for (j = i + 1, vn = 0; j < n; ++ j) if (mat[u[i]][u[j]]) v[vn ++] = u[j]; rr[size] = u[i]; clique(vn, v, mat, size + 1, max, bb, res, rr, c); if (bb) return; } 51 } else if (size > 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; 52 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 53 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; } 7、、、、 图论图论图论图论————连通性连通性连通性连通性 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; 54 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]){ 56 for (st[sp]=-1,a[0]=now,m=1;st[sp]!=i;a[m++]=st[--sp]); dummy(m,a); } } else if (dfn[i] #define MAXN 310 #define _clr(x) memset(x,0xff,sizeof(int)*MAXN) struct edge_t{ int from,to; edge_t* next; }; int hungary(int m,int n,edge_t* list[],int* match1,int* match2){ int s[MAXN],t[MAXN],p,q,ret=0,i,j,k;edge_t* e; for (_clr(match1),_clr(match2),i=0;i=0)) for (_clr(t),s[p=q=0]=i;p<=q&&match1[i]<0;p++) for (e=list[k=s[p]];e&&match1[i]<0;e=e->next) if (t[j=e->to]<0){ s[++q]=match2[j],t[j]=k; if (s[q]<0) for (p=j;p>=0;j=p) match2[j]=k=t[j],p=match1[k],match1[k]=j; } return ret; } 60 8.2 二分图最大匹配二分图最大匹配二分图最大匹配二分图最大匹配(hungary 邻接阵邻接阵邻接阵邻接阵) // 二分图最大匹配,hungary 算法,邻接阵形式,复杂度 O(m*m*n) // 返回最大匹配数,传入二分图大小 m,n 和邻接阵 mat, 非零元素表示有边 //match1,match2 返回一个最大匹配,未匹配顶点 match 值为-1 #include #define MAXN 310 #define _clr(x) memset(x,0xff,sizeof(int)*MAXN) int hungary(int m,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; } 8.3 二分图最大匹配二分图最大匹配二分图最大匹配二分图最大匹配(hungary 正向表正向表正向表正向表) // 二分图最大匹配,hungary 算法,正向表形式,复杂度 O(m*e) // 返回最大匹配数,传入二分图大小 m,n 和正向表 list,buf( 只需一边) //match1,match2 返回一个最大匹配,未匹配顶点 match 值为-1 #include #define MAXN 310 #define _clr(x) memset(x,0xff,sizeof(int)*MAXN) int hungary(int m,int n,int* list,int* buf,int* match1,int* match2){ int s[MAXN],t[MAXN],p,q,ret=0,i,j,k,l; for (_clr(match1),_clr(match2),i=0;i=0)) for (_clr(t),s[p=q=0]=i;p<=q&&match1[i]<0;p++) for (l=list[k=s[p]];l=0;j=p) match2[j]=k=t[j],p=match1[k],match1[k]=j; } 61 return ret; } 8.4 二分图最佳匹配二分图最佳匹配二分图最佳匹配二分图最佳匹配(kuhn_munkras 邻接阵邻接阵邻接阵邻接阵) // 二分图最佳匹配,kuhn munkras 算法,邻接阵形式,复杂度 O(m*m*n) // 返回最佳匹配值,传入二分图大小 m,n 和邻接阵 mat, 表示权值 //match1,match2 返回一个最佳匹配,未匹配顶点 match 值为-1 // 一定注意 m<=n, 否则循环无法终止 // 最小权匹配可将权值取相反数 #include #define MAXN 310 #define inf 1000000000 #define _clr(x) memset(x,0xff,sizeof(int)*n) int kuhn_munkras(int m,int n,int mat[][MAXN],int* match1,int* match2){ int s[MAXN],t[MAXN],l1[MAXN],l2[MAXN],p,q,ret=0,i,j,k; for (i=0;il1[i]?mat[i][j]:l1[i]; for (i=0;i=0;j=p) match2[j]=k=t[j],p=match1[k],match1[k]=j; } if (match1[i]<0){ for (i--,p=inf,k=0;k<=q;k++) for (j=0;jnext) 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; } 63 8.6 一般图匹配一般图匹配一般图匹配一般图匹配(邻接阵邻接阵邻接阵邻接阵) // 一般图最大匹配,邻接阵形式,复杂度 O(n^3) // 返回匹配顶点对数,match 返回匹配,未匹配顶点 match 值为-1 // 传入图的顶点数 n 和邻接阵 mat #define MAXN 100 int aug(int n,int mat[][MAXN],int* match,int* v,int now){ int i,ret=0; v[now]=1; for (i=0;i=2;) if (match[i]<0&&aug(n,mat,match,v,i)) i=0,j-=2; else i++; for (i=j=0;i=0); return j/2; } 8.7 一般图匹配一般图匹配一般图匹配一般图匹配(正向表正向表正向表正向表) 64 // 一般图最大匹配,正向表形式,复杂度 O(n*e) // 返回匹配顶点对数,match 返回匹配,未匹配顶点 match 值为-1 // 传入图的顶点数 n 和正向表 list,buf #define MAXN 100 int aug(int n,int* list,int* buf,int* match,int* v,int now){ int i,t,ret=0; v[now]=1; for (i=list[now];i=2;) if (match[i]<0&&aug(n,list,buf,match,v,i)) i=0,j-=2; else i++; for (i=j=0;i=0); return j/2; } 9、、、、 图论图论图论图论————网络流网络流网络流网络流 9.1 最大最大最大最大流流流流(邻接阵邻接阵邻接阵邻接阵) 65 // 求网络最大流,邻接阵形式 // 返回最大流量,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); 70 } } 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); } 11 、、、、 图论图论图论图论————支撑树支撑树支撑树支撑树 11.1 最小生成树最小生成树最小生成树最小生成树(kruskal 邻接表邻接表邻接表邻接表) 78 // 无向图最小生成树,kruskal 算法,邻接表形式,复杂度 O(mlogm) // 返回最小生成树的长度,传入图的大小 n 和邻接表 list // 可更改边权的类型,edge[][2] 返回树的构造,用边集表示 // 如果图不连通,则对各连通分支构造最小生成树,返回总长度 #include #define MAXN 200 #define inf 1000000000 typedef double elem_t; struct edge_t{ int from,to; elem_t len; edge_t* next; }; #define _ufind_run(x) for(;p[t=x];x=p[x],p[t]=(p[x]?p[x]:x)) #define _run_both _ufind_run(i);_ufind_run(j) struct ufind{ int p[MAXN],t; void init(){memset(p,0,sizeof(p));} void set_friend(int i,int j){_run_both;p[i]=(i==j?0:j);} int is_friend(int i,int j){_run_both;return i==j&&i;} }; #define _cp(a,b) ((a).len<(b).len) struct heap_t{int a,b;elem_t len;}; struct minheap{ heap_t h[MAXN*MAXN]; int n,p,c; void init(){n=0;} void ins(heap_t e){ for (p=++n;p>1&&_cp(e,h[p>>1]);h[p]=h[p>>1],p>>=1); h[p]=e; } int del(heap_t& e){ if (!n) return 0; for (e=h[p=1],c=2;cnext) if (ito) e.a=i,e.b=t->to,e.len=t->len,h.ins(e); while (m #define MAXN 200 #define inf 1000000000 typedef double elem_t; struct edge_t{ int to; elem_t len; }; #define _ufind_run(x) for(;p[t=x];x=p[x],p[t]=(p[x]?p[x]:x)) #define _run_both _ufind_run(i);_ufind_run(j) struct ufind{ int p[MAXN],t; void init(){memset(p,0,sizeof(p));} void set_friend(int i,int j){_run_both;p[i]=(i==j?0:j);} int is_friend(int i,int j){_run_both;return i==j&&i;} }; #define _cp(a,b) ((a).len<(b).len) struct heap_t{int a,b;elem_t len;}; struct minheap{ heap_t h[MAXN*MAXN]; int n,p,c; void init(){n=0;} void ins(heap_t e){ for (p=++n;p>1&&_cp(e,h[p>>1]);h[p]=h[p>>1],p>>=1); 80 h[p]=e; } int del(heap_t& e){ if (!n) return 0; for (e=h[p=1],c=2;c1&&_cp(e,h[p>>1]);h[p]=h[p>>1],p>>=1); h[p]=e; } int del(heap_t& e){ if (!n) return 0; for (e=h[p=1],c=2;cnext) if (!v[t->to]&&t->lento]) pre[t->to]=t->from,min[e.v=t->to]=e.d=t->len,h.ins(e); return ret; } 11.4 最小生成树最小生成树最小生成树最小生成树(prim+binary_heap 正向表正向表正向表正向表) // 无向图最小生成树,prim 算法+二分堆,正向表形式,复杂度 O(mlogm) // 返回最小生成树的长度,传入图的大小 n 和正向表 list,buf // 可更改边权的类型,pre[] 返回树的构造,用父结点表示,根节点(第一个)pre 值为-1 // 必须保证图的连通的! #define MAXN 200 #define inf 1000000000 typedef double elem_t; struct edge_t{ int to; elem_t len; }; 82 #define _cp(a,b) ((a).d<(b).d) struct heap_t{elem_t d;int v;}; struct heap{ heap_t h[MAXN*MAXN]; int n,p,c; void init(){n=0;} void ins(heap_t e){ for (p=++n;p>1&&_cp(e,h[p>>1]);h[p]=h[p>>1],p>>=1); h[p]=e; } int del(heap_t& e){ if (!n) return 0; for (e=h[p=1],c=2;c1&&_cp(e,h[p>>1]);h[map[ind[p]=ind[p>>1]]=p]=h[p>>1],p>>=1); h[map[ind[p]=i]=p]=e; } int del(int i,elem_t& e){ i=map[i];if (i<1||i>n) return 0; for (e=h[p=i];p>1;h[map[ind[p]=ind[p>>1]]=p]=h[p>>1],p>>=1); for (c=2;cnext) if (!v[t->to]&&t->lento]) pre[t->to]=t->from,h.del(t->to,e),h.ins(t->to,min[t->to]=t->len); return ret; } 84 11.6 最小生成树最小生成树最小生成树最小生成树(prim+mapped_heap 正向表正向表正向表正向表) // 无向图最小生成树,prim 算法+映射二分堆,正向表形式,复杂度 O(mlogn) // 返回最小生成树的长度,传入图的大小 n 和正向表 list,buf // 可更改边权的类型,pre[] 返回树的构造,用父结点表示,根节点(第一个)pre 值为-1 // 必须保证图的连通的! #define MAXN 200 #define inf 1000000000 typedef double elem_t; struct edge_t{ int to; elem_t len; }; #define _cp(a,b) ((a)<(b)) struct heap{ elem_t h[MAXN+1]; int ind[MAXN+1],map[MAXN+1],n,p,c; void init(){n=0;} void ins(int i,elem_t e){ for (p=++n;p>1&&_cp(e,h[p>>1]);h[map[ind[p]=ind[p>>1]]=p]=h[p>>1],p>>=1); h[map[ind[p]=i]=p]=e; } int del(int i,elem_t& e){ i=map[i];if (i<1||i>n) return 0; for (e=h[p=i];p>1;h[map[ind[p]=ind[p>>1]]=p]=h[p>>1],p>>=1); for (c=2;c 86 #define MAXN 120 #define inf 1000000000 typedef int elem_t; elem_t edmonds(int n,elem_t mat[][MAXN*2],int* pre){ elem_t ret=0; int c[MAXN*2][MAXN*2],l[MAXN*2],p[MAXN*2],m=n,t,i,j,k; for (i=0;in;pre[k]=pre[m]) for (i=0;i=0&&min[k]+mat[k][i]next) if (min[t->to]==inf) min[que[p++]=t->to]=len*l,pre[t->to]=que[i]; } 12.3 最短路径最短路径最短路径最短路径(单源单源单源单源 dijkstra+bfs 正向表正向表正向表正向表) // 单源最短路径,用于路权相等的情况,dijkstra 优化为 bfs, 正向表形式,复杂度 O(m) // 求出源 s 到所有点的最短路经,传入图的大小 n 和正向表 list,buf, 边权值 len // 返回到各点最短距离 min[] 和路径 pre[],pre[i] 记录 s 到 i 路径上 i 的父结点,pre[s]=-1 // 可更改路权类型,但必须非负且相等! #define MAXN 200 #define inf 1000000000 typedef int elem_t; void dijkstra(int n,int* list,int* buf,elem_t len,int s,elem_t* min,int* pre){ int i,que[MAXN],f=0,r=0,p=1,l=1,t; for (i=0;i1&&_cp(e,h[p>>1]);h[p]=h[p>>1],p>>=1); h[p]=e; } int del(heap_t& e){ if (!n) return 0; for (e=h[p=1],c=2;cnext) if (!v[t->to]&&min[t->from]+t->lento]) pre[t->to]=t->from,min[e.v=t->to]=e.d=min[t->from]+t->len,h.ins(e); } 12.5 最短路径最短路径最短路径最短路径(单源单源单源单源 dijkstra+binary_heap 正向表正向表正向表正向表) // 单源最短路径,dijkstra 算法+二分堆,正向表形式,复杂度 O(mlogm) // 求出源 s 到所有点的最短路经,传入图的大小 n 和正向表 list,buf // 返回到各点最短距离 min[] 和路径 pre[],pre[i] 记录 s 到 i 路径上 i 的父结点,pre[s]=-1 // 可更改路权类型,但必须非负! #define MAXN 200 #define inf 1000000000 typedef int elem_t; struct edge_t{ int to; elem_t len; }; #define _cp(a,b) ((a).d<(b).d) struct heap_t{elem_t d;int v;}; struct heap{ heap_t h[MAXN*MAXN]; int n,p,c; void init(){n=0;} void ins(heap_t e){ for (p=++n;p>1&&_cp(e,h[p>>1]);h[p]=h[p>>1],p>>=1); h[p]=e; } int del(heap_t& e){ if (!n) return 0; for (e=h[p=1],c=2;c1&&_cp(e,h[p>>1]);h[map[ind[p]=ind[p>>1]]=p]=h[p>>1],p>>=1); h[map[ind[p]=i]=p]=e; } int del(int i,elem_t& e){ i=map[i];if (i<1||i>n) return 0; for (e=h[p=i];p>1;h[map[ind[p]=ind[p>>1]]=p]=h[p>>1],p>>=1); for (c=2;cnext) if (!v[t->to]&&min[i]+t->lento]) pre[t->to]=i,h.del(t->to,e),min[t->to]=e=min[i]+t->len,h.ins(t->to,e); } 12.7 最短路径最短路径最短路径最短路径(单源单源单源单源 dijkstra+mapped_heap 正向表正向表正向表正向表) // 单源最短路径,dijkstra 算法+映射二分堆,正向表形式,复杂度 O(mlogn) // 求出源 s 到所有点的最短路经,传入图的大小 n 和正向表 list,buf // 返回到各点最短距离 min[] 和路径 pre[],pre[i] 记录 s 到 i 路径上 i 的父结点,pre[s]=-1 // 可更改路权类型,但必须非负! #define MAXN 200 #define inf 1000000000 typedef int elem_t; struct edge_t{ int to; elem_t len; }; #define _cp(a,b) ((a)<(b)) struct heap{ elem_t h[MAXN+1]; int ind[MAXN+1],map[MAXN+1],n,p,c; void init(){n=0;} void ins(int i,elem_t e){ for (p=++n;p>1&&_cp(e,h[p>>1]);h[map[ind[p]=ind[p>>1]]=p]=h[p>>1],p>>=1); h[map[ind[p]=i]=p]=e; } int del(int i,elem_t& e){ i=map[i];if (i<1||i>n) return 0; for (e=h[p=i];p>1;h[map[ind[p]=ind[p>>1]]=p]=h[p>>1],p>>=1); for (c=2;c=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 96 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 ++; } 100 } 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]; 101 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; 102 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)) 103 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 、、、、 其它其它其它其它 14.1 大数大数大数大数(只能处理正数只能处理正数只能处理正数只能处理正数) #include #include #define DIGIT 4 #define DEPTH 10000 #define MAX 100 typedef int bignum_t[MAX+1]; int read(bignum_t a,istream& is=cin){ char buf[MAX*DIGIT+1],ch; int i,j; memset((void*)a,0,sizeof(bignum_t)); if (!(is>>buf)) return 0; for (a[0]=strlen(buf),i=a[0]/2-1;i>=0;i--) ch=buf[i],buf[i]=buf[a[0]-1-i],buf[a[0]-1-i]=ch; for (a[0]=(a[0]+DIGIT-1)/DIGIT,j=strlen(buf);j1;a[0]--); return 1; } void write(const bignum_t a,ostream& os=cout){ int i,j; for (os<=DEPTH;c[c[0]+1]=c[c[0]]/DEPTH,c[c[0]]%=DEPTH,c[0]++); return comp(a,c); } int comp(const bignum_t a,const int c,const int d,const bignum_t b){ int i,t=0,O=-DEPTH*2; if (b[0]-a[0]d;i--){ t=t*DEPTH+a[i-d]*c-b[i]; if (t>0) return 1; if (t0) return 1; if (t0; 106 } void add(bignum_t a,const bignum_t b){ int i; for (i=1;i<=b[0];i++) if ((a[i]+=b[i])>=DEPTH) a[i]-=DEPTH,a[i+1]++; if (b[0]>=a[0]) a[0]=b[0]; else for (;a[i]>=DEPTH&&i0); } void add(bignum_t a,const int b){ int i=1; for (a[1]+=b;a[i]>=DEPTH&&i=DEPTH;a[a[0]+1]=a[a[0]]/DEPTH,a[a[0]]%=DEPTH,a[0]++); } void sub(bignum_t a,const bignum_t b){ int i; for (i=1;i<=b[0];i++) if ((a[i]-=b[i])<0) a[i+1]--,a[i]+=DEPTH; for (;a[i]<0;a[i]+=DEPTH,i++,a[i]--); for (;!a[a[0]]&&a[0]>1;a[0]--); } void sub(bignum_t a,const int b){ int i=1; for (a[1]-=b;a[i]<0;a[i+1]+=(a[i]-DEPTH+1)/DEPTH,a[i]-=(a[i]-DEPTH+1)/DEPTH*DEPTH,i++); for (;!a[a[0]]&&a[0]>1;a[0]--); } void sub(bignum_t a,const bignum_t b,const int c,const int d){ int i,O=b[0]+d; for (i=1+d;i<=O;i++) if ((a[i]-=b[i-d]*c)<0) a[i+1]+=(a[i]-DEPTH+1)/DEPTH,a[i]-=(a[i]-DEPTH+1)/DEPTH*DEPTH; for (;a[i]<0;a[i+1]+=(a[i]-DEPTH+1)/DEPTH,a[i]-=(a[i]-DEPTH+1)/DEPTH*DEPTH,i++); for (;!a[a[0]]&&a[0]>1;a[0]--); } 107 void mul(bignum_t c,const bignum_t a,const bignum_t b){ int i,j; memset((void*)c,0,sizeof(bignum_t)); for (c[0]=a[0]+b[0]-1,i=1;i<=a[0];i++) for (j=1;j<=b[0];j++) if ((c[i+j-1]+=a[i]*b[j])>=DEPTH) c[i+j]+=c[i+j-1]/DEPTH,c[i+j-1]%=DEPTH; for (c[0]+=(c[c[0]+1]>0);!c[c[0]]&&c[0]>1;c[0]--); } void mul(bignum_t a,const int b){ int i; for (a[1]*=b,i=2;i<=a[0];i++){ a[i]*=b; if (a[i-1]>=DEPTH) a[i]+=a[i-1]/DEPTH,a[i-1]%=DEPTH; } for (;a[a[0]]>=DEPTH;a[a[0]+1]=a[a[0]]/DEPTH,a[a[0]]%=DEPTH,a[0]++); for (;!a[a[0]]&&a[0]>1;a[0]--); } void mul(bignum_t b,const bignum_t a,const int c,const int d){ int i; memset((void*)b,0,sizeof(bignum_t)); for (b[0]=a[0]+d,i=d+1;i<=b[0];i++) if ((b[i]+=a[i-d]*c)>=DEPTH) b[i+1]+=b[i]/DEPTH,b[i]%=DEPTH; for (;b[b[0]+1];b[0]++,b[b[0]+1]=b[b[0]]/DEPTH,b[b[0]]%=DEPTH); for (;!b[b[0]]&&b[0]>1;b[0]--); } void div(bignum_t c,bignum_t a,const bignum_t b){ int h,l,m,i; memset((void*)c,0,sizeof(bignum_t)); c[0]=(b[0]>1;h>l;m=(h+l+1)>>1) if (comp(b,m,i-1,a)) h=m-1; else l=m; for (;!c[c[0]]&&c[0]>1;c[0]--); c[0]=c[0]>1?c[0]:1; } 108 void div(bignum_t a,const int b,int& c){ int i; for (c=0,i=a[0];i;c=c*DEPTH+a[i],a[i]=c/b,c%=b,i--); for (;!a[a[0]]&&a[0]>1;a[0]--); } void sqrt(bignum_t b,bignum_t a){ int h,l,m,i; memset((void*)b,0,sizeof(bignum_t)); for (i=b[0]=(a[0]+1)>>1;i;sub(a,b,m,i-1),b[i]+=m,i--) for (h=DEPTH-1,l=0,b[i]=m=(h+l+1)>>1;h>l;b[i]=m=(h+l+1)>>1) if (comp(b,m,i-1,a)) h=m-1; else l=m; for (;!b[b[0]]&&b[0]>1;b[0]--); for (i=1;i<=b[0];b[i++]>>=1); } int length(const bignum_t a){ int t,ret; for (ret=(a[0]-1)*DIGIT,t=a[a[0]];t;t/=10,ret++); return ret>0?ret:1; } int digit(const bignum_t a,const int b){ int i,ret; for (ret=a[(b-1)/DIGIT+1],i=(b-1)%DIGIT;i;ret/=10,i--); return ret%10; } int zeronum(const bignum_t a){ int ret,t; for (ret=0;!a[ret+1];ret++); for (t=a[ret+1],ret*=DIGIT;!(t%10);t/=10,ret++); return ret; } void comp(int* a,const int l,const int h,const int d){ int i,j,t; for (i=l;i<=h;i++) for (t=i,j=2;t>1;j++) while (!(t%j)) a[j]+=d,t/=j; } 109 void convert(int* a,const int h,bignum_t b){ int i,j,t=1; memset(b,0,sizeof(bignum_t)); for (b[0]=b[1]=1,i=2;i<=h;i++) if (a[i]) for (j=a[i];j;t*=i,j--) if (t*i>DEPTH) mul(b,t),t=1; mul(b,t); } void combination(bignum_t a,int m,int n){ int* t=new int[m+1]; memset((void*)t,0,sizeof(int)*(m+1)); comp(t,n+1,m,1); comp(t,2,m-n,-1); convert(t,m,a); delete []t; } void permutation(bignum_t a,int m,int n){ int i,t=1; memset(a,0,sizeof(bignum_t)); a[0]=a[1]=1; for (i=m-n+1;i<=m;t*=i++) if (t*i>DEPTH) mul(a,t),t=1; mul(a,t); } #define SGN(x) ((x)>0?1:((x)<0?-1:0)) #define ABS(x) ((x)>0?(x):-(x)) int read(bignum_t a,int &sgn,istream& is=cin){ char str[MAX*DIGIT+2],ch,*buf; int i,j; memset((void*)a,0,sizeof(bignum_t)); if (!(is>>str)) return 0; buf=str,sgn=1; if (*buf=='-') sgn=-1,buf++; for (a[0]=strlen(buf),i=a[0]/2-1;i>=0;i--) ch=buf[i],buf[i]=buf[a[0]-1-i],buf[a[0]-1-i]=ch; for (a[0]=(a[0]+DIGIT-1)/DIGIT,j=strlen(buf);j1;a[0]--); if (a[0]==1&&!a[1]) sgn=0; return 1; } 14.2 分数分数分数分数 struct frac{ int num,den; }; double fabs(double x){ return x>0?x:-x; } int gcd(int a,int b){ int t; if (a<0) a=-a; if (b<0) b=-b; if (!b) return a; while (t=a%b) a=b,b=t; return b; } void simplify(frac& f){ int t; if (t=gcd(f.num,f.den)) f.num/=t,f.den/=t; else f.den=1; } frac f(int n,int d,int s=1){ frac ret; if (d<0) ret.num=-n,ret.den=-d; else ret.num=n,ret.den=d; 111 if (s) simplify(ret); return ret; } frac convert(double x){ frac ret; for (ret.den=1;fabs(x-int(x))>1e-10;ret.den*=10,x*=10); ret.num=(int)x; simplify(ret); return ret; } int fraqcmp(frac a,frac b){ int g1=gcd(a.den,b.den),g2=gcd(a.num,b.num); if (!g1||!g2) return 0; return b.den/g1*(a.num/g2)-a.den/g1*(b.num/g2); } frac add(frac a,frac b){ int g1=gcd(a.den,b.den),g2,t; if (!g1) return f(1,0,0); t=b.den/g1*a.num+a.den/g1*b.num; g2=gcd(g1,t); return f(t/g2,a.den/g1*(b.den/g2),0); } frac sub(frac a,frac b){ return add(a,f(-b.num,b.den,0)); } frac mul(frac a,frac b){ int t1=gcd(a.den,b.num),t2=gcd(a.num,b.den); if (!t1||!t2) return f(1,1,0); return f(a.num/t2*(b.num/t1),a.den/t1*(b.den/t2),0); } frac div(frac a,frac b){ return mul(a,f(b.den,b.num,0)); } 112 14.3 矩阵矩阵矩阵矩阵 define MAXN 100 #define fabs(x) ((x)>0?(x):-(x)) #define zero(x) (fabs(x)<1e-10) struct mat{ int n,m; double data[MAXN][MAXN]; }; int mul(mat& c,const mat& a,const mat& b){ int i,j,k; if (a.m!=b.n) return 0; c.n=a.n,c.m=b.m; for (i=0;it) t=fabs(a.data[is[k]=i][js[k]=j]); if (zero(t)) return 0; if (is[k]!=k) for (j=0;j=0;k--){ for (j=0;j0?(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++); 118 days[1]=28; ret.day=a+1; return ret; }