北京大学ACM暑期课讲义-计算几何教程共158页

Computational Geometry
计算几何教程
计算几何的恶心之处
代码长，难写。
需要讨论各种边界情况。后面所介绍的算
法，有些对于边界情况的处理很完美，不 需要再做讨论；有些则不然，需要自行处 理边界情况。
精度误差

二维矢量
2-Dimension Vector
矢量
既有大小又有方向的量。
算法演示
2 被弹出栈，3 进栈。
算法演示
没有点被弹出栈，4 进栈。
算法演示
4 被弹出栈，5 进栈。
算法演示
5 被弹出栈，6 进栈。
算法演示
6 被弹出栈，7 进栈。
算法演示
如此直到下凸壳被找到。
算法演示
倒序进行扫描找到上凸壳。
Graham 扫描的时间复杂度
这太傻缺了。算法的瓶颈在排序，所以时
总时间复杂度 O(N2 log N)。
算法扩展
SPOJ CIRUT
求恰好被覆盖 1 次、2 次、„、N 次的面积。
算法的本质在于对“裸露”部分进行处理。
这种方法也可以用于求很多圆与很多凸多 边形的面积并——但是情况复杂得多，要 讨论的情况也多得多。
三维矢量
3-Dimension Vector
总体来说讨论比较辛苦，写起来不算很难。
精度存在比较大的问题，而且半平面交的
题目一般都容许 O(N2) 算法，所以并不是 很常用。
半平面交练习题
POJ 1279 POJ 3525 POJ 2451 POJ 3384
求多边形的核。核中的点满足其到多边形边 界上任一点之间的线段只与多边形在那一点相交。 求凸多边形最大内切圆半径。 裸的半平面交。朱泽园专门为他的算法出的 题，要求 O(N log N) 算法。 求在凸多边形内部放两个半径为 r 的圆所能 覆盖的最大面积。两圆可以重叠，但不能与多边形相交。 需要用到旋转卡壳。

/*计算几何目录㈠点的基本运算1. 平面上两点之间距离12. 判断两点是否重合13. 矢量叉乘14. 矢量点乘25. 判断点是否在线段上26. 求一点饶某点旋转后的坐标27. 求矢量夹角2㈡线段及直线的基本运算1. 点与线段的关系32. 求点到线段所在直线垂线的垂足43. 点到线段的最近点44. 点到线段所在直线的距离45. 点到折线集的最近距离46. 判断圆是否在多边形内57. 求矢量夹角余弦58. 求线段之间的夹角59. 判断线段是否相交610.判断线段是否相交但不交在端点处611.求线段所在直线的方程612.求直线的斜率713.求直线的倾斜角714.求点关于某直线的对称点715.判断两条直线是否相交及求直线交点716.判断线段是否相交，如果相交返回交点7㈢多边形常用算法模块1. 判断多边形是否简单多边形82. 检查多边形顶点的凸凹性93. 判断多边形是否凸多边形94. 求多边形面积95. 判断多边形顶点的排列方向，方法一106. 判断多边形顶点的排列方向，方法二107. 射线法判断点是否在多边形内108. 判断点是否在凸多边形内119. 寻找点集的graham算法1210.寻找点集凸包的卷包裹法1311.判断线段是否在多边形内1412.求简单多边形的重心1513.求凸多边形的重心1714.求肯定在给定多边形内的一个点1715.求从多边形外一点出发到该多边形的切线1816.判断多边形的核是否存在19㈣圆的基本运算1 .点是否在圆内202 .求不共线的三点所确定的圆21㈤矩形的基本运算1.已知矩形三点坐标，求第4点坐标22㈥常用算法的描述22㈦补充1．两圆关系：242．判断圆是否在矩形内：243．点到平面的距离：254．点是否在直线同侧：255．镜面反射线：256．矩形包含：267．两圆交点：278．两圆公共面积：289. 圆和直线关系：2910. 内切圆：3011. 求切点：3112. 线段的左右旋：3113．公式：32*//* 需要包含的头文件*/#include <cmath >/* 常用的常量定义*/const double INF = 1E200const double EP = 1E-10const int MAXV = 300const double PI = 3.14159265/* 基本几何结构*/struct POINT{double x;double y;POINT(double a=0, double b=0) { x=a; y=b;} //constructor};struct LINESEG{POINT s;POINT e;LINESEG(POINT a, POINT b) { s=a; e=b;}LINESEG() { }};struct LINE // 直线的解析方程a*x+b*y+c=0 为统一表示，约定a >= 0{double a;double b;double c;LINE(double d1=1, double d2=-1, double d3=0) {a=d1; b=d2; c=d3;}};/*********************** ** 点的基本运算** ***********************/double dist(POINT p1,POINT p2) // 返回两点之间欧氏距离{return( sqrt( (p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y) ) );}bool equal_point(POINT p1,POINT p2) // 判断两个点是否重合{return ( (abs(p1.x-p2.x)<EP)&&(abs(p1.y-p2.y)<EP) );}/****************************************************************************** r=multiply(sp,ep,op),得到(sp-op)和(ep-op)的叉积r>0：ep在矢量opsp的逆时针方向；r=0：opspep三点共线；r<0：ep在矢量opsp的顺时针方向******************************************************************************* /double multiply(POINT sp,POINT ep,POINT op){return((sp.x-op.x)*(ep.y-op.y)-(ep.x-op.x)*(sp.y-op.y));}/*r=dotmultiply(p1,p2,op),得到矢量(p1-op)和(p2-op)的点积，如果两个矢量都非零矢量r<0：两矢量夹角为锐角；r=0：两矢量夹角为直角；r>0：两矢量夹角为钝角******************************************************************************* /double dotmultiply(POINT p1,POINT p2,POINT p0){return ((p1.x-p0.x)*(p2.x-p0.x)+(p1.y-p0.y)*(p2.y-p0.y));}/****************************************************************************** 判断点p是否在线段l上条件：(p在线段l所在的直线上) && (点p在以线段l为对角线的矩形内)******************************************************************************* /bool online(LINESEG l,POINT p){return( (multiply(l.e,p,l.s)==0) &&( ( (p.x-l.s.x)*(p.x-l.e.x)<=0 )&&( (p.y-l.s.y)*(p.y-l.e.y)<=0 ) ) ); }// 返回点p以点o为圆心逆时针旋转alpha(单位：弧度)后所在的位置POINT rotate(POINT o,double alpha,POINT p){POINT tp;p.x-=o.x;p.y-=o.y;tp.x=p.x*cos(alpha)-p.y*sin(alpha)+o.x;tp.y=p.y*cos(alpha)+p.x*sin(alpha)+o.y;return tp;}/* 返回顶角在o点，起始边为os，终止边为oe的夹角(单位：弧度)角度小于pi，返回正值角度大于pi，返回负值可以用于求线段之间的夹角原理：r = dotmultiply(s,e,o) / (dist(o,s)*dist(o,e))r'= multiply(s,e,o)r >= 1 angle = 0;r <= -1 angle = -PI-1<r<1 && r'>0 angle = arccos(r)-1<r<1 && r'<=0 angle = -arccos(r)*/double angle(POINT o,POINT s,POINT e){double cosfi,fi,norm;double dsx = s.x - o.x;double dsy = s.y - o.y;double dex = e.x - o.x;double dey = e.y - o.y;cosfi=dsx*dex+dsy*dey;norm=(dsx*dsx+dsy*dsy)*(dex*dex+dey*dey);cosfi /= sqrt( norm );if (cosfi >= 1.0 ) return 0;if (cosfi <= -1.0 ) return -3.1415926;fi=acos(cosfi);if (dsx*dey-dsy*dex>0) return fi; // 说明矢量os 在矢量oe的顺时针方向return -fi;}/*****************************\* ** 线段及直线的基本运算** *\*****************************//* 判断点与线段的关系,用途很广泛本函数是根据下面的公式写的，P是点C到线段AB所在直线的垂足AC dot ABr = ---------||AB||^2(Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay)= -------------------------------L^2r has the following meaning:r=0 P = Ar=1 P = Br<0 P is on the backward extension of ABr>1 P is on the forward extension of AB0<r<1 P is interior to AB*/double relation(POINT p,LINESEG l){LINESEG tl;tl.s=l.s;tl.e=p;return dotmultiply(tl.e,l.e,l.s)/(dist(l.s,l.e)*dist(l.s,l.e));}// 求点C到线段AB所在直线的垂足PPOINT perpendicular(POINT p,LINESEG l){double r=relation(p,l);POINT tp;tp.x=l.s.x+r*(l.e.x-l.s.x);tp.y=l.s.y+r*(l.e.y-l.s.y);return tp;}/* 求点p到线段l的最短距离,并返回线段上距该点最近的点np注意：np是线段l上到点p最近的点，不一定是垂足*/double ptolinesegdist(POINT p,LINESEG l,POINT &np){double r=relation(p,l);if(r<0){np=l.s;return dist(p,l.s);}if(r>1){np=l.e;return dist(p,l.e);}np=perpendicular(p,l);return dist(p,np);}// 求点p到线段l所在直线的距离,请注意本函数与上个函数的区别double ptoldist(POINT p,LINESEG l){return abs(multiply(p,l.e,l.s))/dist(l.s,l.e);}/* 计算点到折线集的最近距离,并返回最近点.注意：调用的是ptolineseg()函数*/double ptopointset(int vcount,POINT pointset[],POINT p,POINT &q) {int i;double cd=double(INF),td;LINESEG l;POINT tq,cq;for(i=0;i<vcount-1;i++)l.s=pointset[i];l.e=pointset[i+1];td=ptolinesegdist(p,l,tq);if(td<cd){cd=td;cq=tq;}}q=cq;return cd;}/* 判断圆是否在多边形内.ptolineseg()函数的应用2 */bool CircleInsidePolygon(int vcount,POINT center,double radius,POINT polygon[]){POINT q;double d;q.x=0;q.y=0;d=ptopointset(vcount,polygon,center,q);if(d<radius||fabs(d-radius)<EP)return true;elsereturn false;}/* 返回两个矢量l1和l2的夹角的余弦(-1 --- 1)注意：如果想从余弦求夹角的话，注意反余弦函数的定义域是从0到pi */double cosine(LINESEG l1,LINESEG l2){return (((l1.e.x-l1.s.x)*(l2.e.x-l2.s.x) +(l1.e.y-l1.s.y)*(l2.e.y-l2.s.y))/(dist(l1.e,l1.s)*dist(l2.e,l2.s))) );}// 返回线段l1与l2之间的夹角单位：弧度范围(-pi，pi)double lsangle(LINESEG l1,LINESEG l2){POINT o,s,e;o.x=o.y=0;s.x=l1.e.x-l1.s.x;s.y=l1.e.y-l1.s.y;e.x=l2.e.x-l2.s.x;e.y=l2.e.y-l2.s.y;return angle(o,s,e);// 如果线段u和v相交(包括相交在端点处)时，返回true////判断P1P2跨立Q1Q2的依据是：( P1 - Q1 ) ×( Q2 - Q1 ) * ( Q2 - Q1 ) ×( P2 - Q1 ) >= 0。

三计算模型
规定RAM的原始运算如下： 1）算术运算 2）两个实数之间的比较 3）间接寻址及求根运算，三角函数运算，指 数函数运算，对数函数运算。 执行每种运算耗费一个时间单位。
四常见计算几何基本问题
• 矢量的概念
• • • • • • • • • • • • • • • • • • • • • • • • • • • •
基本问题（续）
• 在实际编程中，没有必要计算所有的交点，首先应判断线段和多边形的边是否内 交，倘若线段和多边形的某条边内交则线段一定在多边形外；如果线段和多边形 的每一条边都不内交，则线段和多边形的交点一定是线段的端点或者多边形的顶 点，只要判断点是否在线段上就可以了。
得出算法如下： if 线端PQ的端点不都在多边形内 then return false; 点集pointSet初始化为空; for 多边形的每条边s do if 线段的某个端点在s上 then 将该端点加入pointSet; else if s的某个端点在线段PQ上 then 将该端点加入pointSet; else if s和线段PQ相交 // 这时候已经可以肯定是内交了 then return false; 将pointSet中的点按照X-Y坐标排序; for pointSet中每两个相邻点 pointSet[i] , pointSet[ i+1] do if pointSet[i] , pointSet[ i+1] 的中点不在多边形中 then return false; return true; 这个过程中的排序因为交点数目肯定远小于多边形的顶点数目n，所以最多是 常数级的复杂度，几乎可以忽略不计。因此算法的时间复杂度也是O(n)。
基本问题（续）
基本问题（续）

E问题
OUTPUT FORMAT
输出应当有F+1行，每行一个整数，依次表示路径经 过的顶点号。注意数据可能有多组解，但是只有上面题目 要求的那一组解是认为正确的。
SAMPLE OUTPUT(fence.out) 1 2 3 4 2 5 4 6 5 7
E问题
很显然，这是一个欧拉路径问题，我们要做的就是读入栅栏的 构图后，找到图中的一条欧拉路径。
图的连通性
4. 有向图的强连通分支
在下面的几页中，我们可以看到求图的 强连通分支的实例。
首先，图(a)为有向图G，其中的阴影部 分是G的强连通分支，在对图G进行DFS 的过程中，我们对每个顶点都标出了其 开始搜索时刻preOrder与完成时刻 postOrder，黑色边为DFS搜索树树枝；
可以看到，图G共有 4个强连通分支：
{a,b,e}
{c,d}
{f,g}
{h}
图的连通性
4. 有向图的强连通分量
(b)图中G的转置图G*。图中说明了求 强连通分支算法第3部计算出的深度优 先树，其中黑色边是树枝。每个强连通
子图对应于一棵深度优先树。图中黑色 顶点b,c,g和h是强连通子图中每个顶点 的祖先，这些顶点也是对G*进行深度 优先搜索所产生的深度优先树的树根。
int map[MAXV][MAXV]; //map[i][j]记录顶点i和顶点j之 间的路径数目
int deg[MAXV]; int path[MAXE]; 径
//deg[i]记录顶点i的度数 //path数组用来存放找到的欧拉路
int fn,minv,maxv,pathnum=0; //minv为顶点最小编号，maxv为顶点最大编号
G有欧拉回路(G为欧拉图)：G连通，G中均为偶度顶点。

Class 1: Expectations, variances, and basics of estimationBasics of matrix (1)I. Organizational Matters(1)Course requirements:1)Exercises: There will be seven (7) exercises, the last of which is optional. Eachexercise will be graded on a scale of 0-10. In addition to the graded exercise, ananswer handout will be given to you in lab sections.2)Examination: There will be one in-class, open-book examination.(2)Computer software: StataII. Teaching Strategies(1) Emphasis on conceptual understanding.Yes, we will deal with mathematical formulas, actually a lot of mathematical formulas. But, I do not want you to memorize them. What I hope you will do, is to understand the logic behind the mathematical formulas.(2) Emphasis on hands-on research experience.Yes, we will use computers for most of our work. But I do not want you to become a computer programmer. Many people think they know statistics once they know how to run a statistical package. This is wrong. Doing statistics is more than running computer programs. What I will emphasize is to use computer programs to your advantage in research settings. Computer programs are like automobiles. The best automobile is useless unless someone drives it. You will be the driver of statistical computer programs.(3) Emphasis on student-instructor communication.I happen to believe in students' judgment about their own education. Even though I will be ultimately responsible if the class should not go well, I hope that you will feel part of the class and contribute to the quality of the course. If you have questions, do not hesitate to ask in class. If you have suggestions, please come forward with them. The class is as much yours as mine.Now let us get to the real business.III(1). Expectation and VarianceRandom Variable: A random variable is a variable whose numerical value is determined by the outcome of a random trial.Two properties: random and variable.A random variable assigns numeric values to uncertain outcomes. In a common language, "give a number". For example, income can be a random variable. There are many ways to do it. You can use the actual dollar amounts.In this case, you have a continuous random variable. Or you can use levels of income, such as high, median, and low. In this case, you have an ordinal random variable [1=high,2=median, 3=low]. Or if you are interested in the issue of poverty, you can have a dichotomous variable: 1=in poverty, 0=not in poverty.In sum, the mapping of numeric values to outcomes of events in this way is theessence of a random variable.Probability Distribution: The probability distribution for a discrete random variable Xassociates with each of the distinct outcomes x i (i = 1, 2,..., k ) a probability P (X = x i ). Cumulative Probability Distribution: The cumulative probability distribution for a discreterandom variable X provides the cumulative probabilities P (X ≤ x ) for all values x .Expected Value of Random Variable: The expected value of a discrete random variable X isdenoted by E {X } and defined:E {X } = x i i k=∑1P (x i )where : P (x i ) denotes P (X = x i ). The notation E { } (read “expectation of”) is called theexpectation operator.In common language, expectation is the mean. But the difference is that expectation is a concept for the entire population that you never observe. It is the result of the infinitenumber of repetitions. For example, if you toss a coin, the proportion of tails should be .5 in the limit. Or the expectation is .5. Most of the times you do not get the exact .5, but a number close to it.Conditional ExpectationIt is the mean of a variable conditional on the value of another random variable.Note the notation: E(Y|X).In 1996, per-capita average wages in three Chinese cities were (in RMB):Shanghai: 3,778Wuhan: 1,709Xi ’an: 1,155Variance of Random Variable: The variance of a discrete random variable X is denoted by V {X } and defined:V {X } = i k =∑1(x i - E {X })2 P (x i )where : P (x i ) denotes P (X = x i ). The notation V { } (read “variance of”) is called the variance operator.Since the variance of a random variable X is a weighted average of the squared deviations,(X - E {X })2 , it may be defined equivalently as an expected value: V {X } = E {(X - E {X })2}. An algebraically identical expression is: V {X} = E {X 2} - (E {X })2.Standard Deviation of Random Variable: The positive square root of the variance of X is called the standard deviation of X and is denoted by σ{X }:σ {X } =V X {}The notation σ{ } (read “standard deviation of”) is called the standard deviation operator. Standardized Random Variables: If X is a random variable with expected value E {X } and standard deviation σ{X }, then:Y =}{}{X X E X σ-is known as the standardized form of random variable X .Covariance: The covariance of two discrete random variables X and Y is denoted by Cov {X,Y } and defined:Cov {X, Y } = ({})({})(,)xE X y E Y P x y i j i j j i --∑∑where: P (x i , y j ) denotes P Xx Y y i j (=⋂= ) The notation of Cov { , } (read “covariance of”) is called the covariance operator .When X and Y are independent, Cov {X, Y } = 0.Cov {X, Y } = E {(X - E {X })(Y - E {Y })}; Cov {X, Y } = E {XY } - E {X }E {Y }(Variance is a special case of covariance.)Coefficient of Correlation: The coefficient of correlation of two random variables X and Y isdenoted by ρ{X,Y } (Greek rho) and defined:ρσσ{,}{,}{}{}X Y Cov X Y X Y =where: σ{X } is the standard deviation of X; σ{Y } is the standard deviation of Y; Cov is the covariance of X and Y.Sum and Difference of Two Random Variables: If X and Y are two random variables, then the expected value and the variance of X + Y are as follows:Expected Value : E {X+Y } = E {X } + E {Y };Variance : V {X+Y } = V {X } + V {Y }+ 2 Cov (X,Y ).If X and Y are two random variables, then the expected value and the variance of X - Y are as follows:Expected Value : E {X - Y } = E {X } - E {Y };Variance : V {X - Y } = V {X } + V {Y } - 2 Cov (X,Y ).Sum of More Than Two Independent Random Variables: If T = X 1 + X 2 + ... + X s is the sumof s independent random variables, then the expected value and the variance of T are as follows:Expected Value: E T E X i i s {}{}==∑1; Variance: V T V X i i s {}{}==∑1III(2). Properties of Expectations and Covariances:(1) Properties of Expectations under Simple Algebraic Operations)()(x bE a bX a E +=+This says that a linear transformation is retained after taking an expectation.bX a X +=*is called rescaling: a is the location parameter, b is the scale parameter.Special cases are:For a constant: a a E =)(For a different scale: )()(X E b bX E =, e.g., transforming the scale of dollars into thescale of cents.(2) Properties of Variances under Simple Algebraic Operations)()(2X V b bX a V =+This says two things: (1) Adding a constant to a variable does not change the variance of the variable; reason: the definition of variance controls for the mean of the variable[graphics]. (2) Multiplying a constant to a variable changes the variance of the variable by a factor of the constant squared; this is to easy prove, and I will leave it to you. This is the reason why we often use standard deviation instead of variance2x x σσ=is of the same scale as x.(3) Properties of Covariance under Simple Algebraic OperationsCov(a + bX, c + dY) = bd Cov(X,Y).Again, only scale matters, location does not.(4) Properties of Correlation under Simple Algebraic OperationsI will leave this as part of your first exercise: ),(),(Y X dY c bX a ρρ=++That is, neither scale nor location affects correlation.IV: Basics of matrix.1. DefinitionsA. MatricesToday, I would like to introduce the basics of matrix algebra. A matrix is a rectangular array of elements arranged in rows and columns:11121211.......m n nm x x x x X x x ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦Index: row index, column index.Dimension: number of rows x number of columns (n x m)Elements: are denoted in small letters with subscripts.An example is the spreadsheet that records the grades for your home work in the following way:Name 1st 2nd ....6thA 7 10 (9)B 6 5 (8)... ... ......Z 8 9 (8)This is a matrix.Notation: I will use Capital Letters for Matrices.B. VectorsVectors are special cases of matrices:If the dimension of a matrix is n x 1, it is a column vector:⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=n x x x x (21)If the dimension is 1 x m, it is a row vector: y' = | 1y 2y .... m y |Notation: small underlined letters for column vectors (in lecture notes)C. TransposeThe transpose of a matrix is another matrix with positions of rows and columns beingexchanged symmetrically.For example: if⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⨯nm n m m n x x x x x x X 12111211)( (1121112)()1....'...n m n m nm x x x x X x x ⨯⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦It is easy to see that a row vector and a column vector are transposes of each other. 2. Matrix Addition and SubtractionAdditions and subtraction of two matrices are possible only when the matrices have the same dimension. In this case, addition or subtraction of matrices forms another matrix whoseelements consist of the sum, or difference, of the corresponding elements of the two matrices.⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡±±±±±=Y ±X mn nm n n m m y x y x y x y x y x (11)2121111111 Examples:⎥⎦⎤⎢⎣⎡=A ⨯4321)22(⎥⎦⎤⎢⎣⎡=B ⨯1111)22(⎥⎦⎤⎢⎣⎡=B +A =⨯5432)22(C 3. Matrix MultiplicationA. Multiplication of a scalar and a matrixMultiplying a scalar to a matrix is equivalent to multiplying the scalar to each of the elements of the matrix.11121211Χ...m n nm cx c cx cx ⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦ B. Multiplication of a Matrix by a Matrix (Inner Product)The inner product of matrix X (a x b) and matrix Y (c x d) exists if b is equal to c. The inner product is a new matrix with the dimension (a x d). The element of the new matrix Z is:c∑=kj ik ijy x zk=1Note that XY and YX are very different. Very often, only one of the inner products (XY and YX) exists.Example:⎥⎦⎤⎢⎣⎡=4321)22(x A⎥⎦⎤⎢⎣⎡=10)12(x BBA does not exist. AB has the dimension 2x1⎥⎦⎤⎢⎣⎡=42ABOther examples:If )53(x A , )35(x B , what is the dimension of AB? (3x3)If )53(x A , )35(x B , what is the dimension of BA? (5x5)If )51(x A , )15(x B , what is the dimension of AB? (1x1, scalar)If )53(x A , )15(x B , what is the dimension of BA? (nonexistent)4. Special MatricesA. Square Matrix)(n n A ⨯B. Symmetric MatrixA special case of square matrix.For )(n n A ⨯, ji ij a a =. All i, j .A' = AC. Diagonal MatrixA special case of symmetric matrix⎥⎥⎥⎥⎦⎢⎢⎢⎢⎣=X nn x x 0 (2211)D. Scalar Matrix0....0c c c c ⎡⎤⎢⎥⎢⎥=I ⎢⎥⎢⎥⎣⎦E. Identity MatrixA special case of scalar matrix⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=I 10 (101)Important: for r r A ⨯AI = IA = AF. Null (Zero) MatrixAnother special case of scalar matrix⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=O 00 (000)From A to E or F, cases are nested from being more general towards being more specific.G. Idempotent MatrixLet A be a square symmetric matrix. A is idempotent if ....32=A =A =AH. Vectors and Matrices with elements being oneA column vector with all elements being 1,⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⨯1 (111)r A matrix with all elements being 1, ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⨯1...1...111...11rr J Examples let 1 be a vector of n 1's: )1(1⨯n 1'1 = )11(⨯n11' = )(n n J ⨯I. Zero Vector A zero vector is ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⨯0 (001)r 5. Rank of a MatrixThe maximum number of linearly independent rows is equal to the maximum number oflinearly independent columns. This unique number is defined to be the rank of the matrix. For example,⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=B 542211014321 Because row 3 = row 1 + row 2, the 3rd row is linearly dependent on rows 1 and 2. The maximum number of independent rows is 2. Let us have a new matrix:⎥⎦⎤⎢⎣⎡=B 11014321* Singularity: if a square matrix A of dimension ()n n ⨯has rank n, the matrix is nonsingular. If the rank is less than n, the matrix is then singular.。

快速排斥试验引申
• 思考： • 简单的改造一下就可以用来判断一个矩形 是否在另一个矩形内部 • NKOJ 1177 rectangles
凸包
• convex hull • 是指对于平面上给定的一些点，选取一个 最小的凸多边形使得这些点或者在其内部， 或者在其边上
求凸包的Graham扫描法
• 对于一个有三个或以上点的点集Q • 令p0为Q中Y-X坐标排序下最小的点 设<p1,p2,...pm>为对其余点按以p0为中心的极角 逆时针排序所得的点集（如果有多个点有相同的 极角，除了距p0最远的点外全部移除） 压p0进栈S，压p1进栈S，压p2进栈S for(i=3;i<=m;++i){ while(由S的栈顶元素的下一个元素、S的栈顶 元素以及pi构成的折线段不拐向左侧){ 对S弹栈 压pi进栈S } } return S
计算几何
需要注意的细节
• 常用头文件#include<math.h> • 计算几何中一般来说使用double型比较频 繁，请注意数据类型的选择，该用实数的 时候就用double，而float容易失去精度。 • 判断double型的x是否为0，应当用x<eps && x>-eps（或者fabs(x)<eps），其中eps 代表某个精度，常常取eps=0.000001，还 有其他类似情况也要注意double类型的精 度问题
需要注意的细节
• 圆周率取3.141592654或者更精确，或者用 acos(-1) • 角度制和弧度制的转换，C/C++中的三角函 数均为弧度制 • 尽量少用除法，开方，三角函数，容易失 去精度。用除法时注意除数不为0
向量及其运算

– 如果当前是我的回合，无论怎么做，轮到对方时， 对方都存在一种必胜策略，那么当前我面对的局面 就是必败局面
• 例：取石子问题
– 有一堆n个石子 两人轮流取 每次可以取1..m个 没得取的判负 在游戏中，n不断变化，m是定值 一个局面可以用当前的石子数n来表示 n = 0 必败 n = 1..m都是必胜局面 n = m+1 是必败局面（为什么？） n = m+2.. 2m+1 是必胜局面（为什么？）
– 附加问题：先手必胜时，输出当前一步可选 步
• 最终局面都是P局面 • 对于一个局面,若至少有一种操作使它变为P局 面,则它是N局面
– 如果当前是我的回合，我有一种操作选择，使得下 回合对方面对的局面是个必输局面，那么我的局面 就是必胜局面 – 问题：如何寻找这一步操作选择？（输出解）
• 对于一个局面,无论如何操作都使它变为N局面, 则它是P局面
• 双方每步可以从并行存在的多个子游戏中选择 一个操作
– 维度可以进一步变得更多 – 可能存在一个子游戏，由同一个玩家连续两次对其 操作
• 对sg函数来说：怎么变都无所谓，它只关心所 有子局面的异或值
– Sg(x)的变量x永远是一维的。
竞赛中的博弈题策略
• 利用最简单的N-P图染色来判定 • 博弈搜索：用搜索方式来对图染色
• 同样可以用拓扑序对节点计算sg值
SG函数与NP局面的关系
• 初始节点（必须是P节点）没有出边， sg = 0 • 若v有边指向某sg(vi) = 0的vi，则sg(v)>0 • 若v没有边指向某sg(vi) = 0的vi，则sg(v) =0
• 所有P节点sg = 0 • 所有N节点sg > 0
– 可以看做α-β剪枝的特例

在边界上情况的特殊判断
◦ 一个避免精度问题的加速方法
Authored by Will, IIIS, THU 2020/1/21
8
预备知识：距离与交点的计算
点的距离计算
◦ 线段，直线 ◦ 圆，多边形
线段相关
◦ 线段距离？ ◦ 线段的交点？
圆相关
◦ 直线与圆的距离与交点 ◦ 点与圆的切点 ◦ 两圆的公切点
计算方法
◦ 类似Monotone Chain O(NlogN)
更一般的性质与问题
◦ 简单正交多边形的正交凸包 ◦ 口袋的面积 ◦ 最小的周长
Authored by Will, IIIS, THU 2020/1/21
16
凸包上的序——凸性
凸包的定向 直线与凸包 矩形覆盖问题 点与凸包 凸包与凸包 凸包上的动态计算问题 凸包的划分问题
◦ 线段的相交测试 ◦ 面积的计算
右手方向与左手方向
Authored by Will, IIIS, THU 2020/1/21
7
预备知识：另一个例子
判断点在多边形内部
◦ 凸多边形 ◦ 简单多边形
射线法——效率高！
◦ 共线的处理？
随机斜率 纵切变换
转角法——优美！
◦ 精度问题？
浅谈一些几何问题
清华大学 交叉信息院 吴翼 Will
Authored by Will, IIIS, THU 2020/1/21
1
前言：计算几何与计算机科学
机械自动控制 地理信息贪测 计算机辅助设计 模式识别
◦ 人脸识别 ◦ 相关图像处理
游戏领域
◦ 碰撞检测 ◦ 光线追踪 ◦ 三维场景渲染

acm-icpc暑期课_最短路北京⼤学暑期课《ACM/ICPC 竞赛训练》最短路算法本讲义为四处抄袭后改编⽽成，来源已不可考。

仅⽤于内部授课北京⼤学信息学院郭炜 guo_wei@/doc/d83c4f84bb68a98270fefa4a.html/doc/d83c4f84bb68a98270fefa4a.html /guoweiofpkuDijkstra 算法●解决⽆负权边的带权有向图或⽆向图的单源最短路问题●贪⼼思想，若离源点s 前k-1近的点已经被确定，构成点集P ，那么从s 到离s 第k 近的点t 的最短路径，{s,p 1,p 2…p i ,t}满⾜s,p 1,p 2…p i ∈P 。

●否则假设pi ?P ，则因为边权⾮负，pi 到t 的路径≥0，则d[pi]≤d[t],pi 才是第k 近。

将pi 看作t ，重复上⾯过程，最终⼀定会有找不到pi 的情况●d[i]=min(d[p i ]+cost(p i ,i)),i ?P,p i ∈P d[t]=min(d[i]) ,i ?P基本思想●初始令d[s]=0，d[i]=+∞，P=?●找到点i?P，且d[i]最⼩●把i添⼊P，对于任意j?P，若d[i]+cost(i,j)d[j]=d[i]+cost(i,j)。

●⽤邻接表，不优化，时间复杂度O(V2+E)●Dijkstra+堆的时间复杂度 o(ElgV)●⽤斐波那契堆可以做到O(VlogV+E)●若要输出路径，则设置prev数组记录每个节点的前趋点，在d[i]更新时更新prev[i]4 23 150200200 100080005000v Dist[v]0 01234 50100001000300025025082505250源点0加⼊P后：Dijkstra's AlgorithmDijkstra 算法也适⽤于⽆向图。

但不适⽤于有负权边的图。

231-234d[1,2] = 2但⽤Dijkstra 算法求得 d[1,2] = 3有N个孩⼦（N<=3000)分糖果。

北京大学暑期课《ACM/ICPC竞赛训练》北京大学信息学院郭炜guo_wei@/guoweiofpku课程网页：/summerschool/pku_acm_train.htm信息科学技术学院《程序设计与算法》ACM/ICPC中的数学郭炜/林舒( a + b ) mod c =((a mod c) + (b mod c)) mod c( a * b ) mod c =((a mod c) * (b mod c)) mod c消去律：若 gcd(c,p) = 1,则ac ≡ bc mod p => a ≡ b mod p消去律证明：ac ≡ bc mod p => ac = bc + kp => c(a-b) = kpc和p没有除1以外的公因子 =>1) c能整除k 或 2) a = b如果2不成立，则c|kpc和p没有公因子 => c|k，所以k = ck'=> c(a-b)=kp可以表示为c(a-b) =ck'p因此a-b = k'p，得出a ≡ b mod pint Pow(int a,int b){ //快速求a^b ，复杂度 log(b)if(b == 0)return 1;if(b & 1) { //b是奇数return a * Pow(a,b-1);}else {int t = Pow(a,b/2);return t * t;}}int Pow(int a,int b){ //快速求a^b ，复杂度 log(b) int result = 1;int base = a;while(b) {if( b & 1)result *= base;base *= base;b >>= 1;}return result;}快速幂取模int PowMod(int a,int b,int c){//快速求 a^b % c ，要避免计算中间结果溢出int result = 1;int base = a%c;while(b) {if( b & 1)result = (result * base)%c;base = (base * base) % c;b >>= 1;}return result;}Sn= a+a2+...+a n要求 Snmod p如果用公式算，可能溢出，因此用二分法求1) 若 n是偶数Sn= a+...+a n/2 + a n/2+1 + a n/2+2 +...+ a n/2+n/2 =(a+...+a n/2) + a n/2(a+...+a n/2)=Sn/2+ a n/2Sn/2=(1+a n/2)Sn/22) 若n是奇数= a+...+a(n-1)/2 + a(n-1)/2+1 +...Sn+ a(n-1)/2+(n-1)/2 + a(n-1)/2+(n-1)/2 + 1 =S+ a(n-1)/2(a+...+a(n-1)/2)+a n (n-1)/2+a n=(1+a(n-1)/2)S(n-1)/2等比数列二分求和取模int PowSumMod(int a,int n,int p){// return (a+ a^2 + ... + a^n) Mod p;if( n == 1)return a%p;if( n %2 == 0)return (1+PowMod(a,n/2,p))*PowSumMod(a,n/2,p) % p;elsereturn ((1+PowMod(a,(n-1)/2,p)) * PowSumMod(a,(n-1)/2,p) + PowMod(a,n,p)) % p;}POJ3233 Matrix Power Series矩阵快速幂+等比数列二分求和取模给一个n×n的整数矩阵A和正整数k,m, 令S = A + A2 + A3+ … +A k求 S mod m (S的每个都 mod m)n (n≤ 30),k (k≤ 109) and m (m < 104).struct Matrix{T a[32][32];int r; //行列数Mat(int rr):r(rr) { }void MakeI() { //变为单位矩阵memset(a,0,sizeof(a));for(int i = 0;i < r; ++i)a[i][i] = 1;}};Matrix Pow(const Matrix & m,int k,int p){ //求m k mod pint r = m.r;Matrix result(r);result.MakeI(); //MakeI是将result变为单位矩阵Matrix base = m;while(k) {if( k & 1)result = Product(result,base,p); //result*base mod pk >>= 1;base = Product(base,base,p);}return result;} 14欧几里得算法求最大公约数int Gcd(int a,int b){if( b == 0)return a;return Gcd(b,a%b);}最小公倍数：lcm(a,b) = a*b/gcd(a,b)ax+by=gcd(a,b) 有整数解(x,y)<=>bx+(a%b)y = gcd(a,b) 有整数解(x1,y1)，且x = y1, y = x1-[a/b]*y1因此，可以在求gcd(a,b)的同时，对 ax+by=gcd(a,b)求解int GcdEx(int a,int b,int &x ,int & y) //求 ax+by=gcd(a,b)的整数解,返回gcd(a,b) {if( b == 0) {x = 1;y = 0;return a;}int x1,y1;int gcd = GcdEx(b,a%b,x1,y1);x = y1;y = x1 - a/b * y1;return gcd;}• ax+by=c 有解的充要条件是 gcd(a,b)|c•设 d = gcd(a,b), k = c/d，ax+by = d的解是 (x1,y1) 则ax+by = c的解集是：x = k*x1 + t*(b/d)y = k*y1 – t*(a/d) t为任意整数18•求ax≡c(mod b)等价于求 ax+by = c •设 d = gcd(a,b), k = c/d，ax+by = d的解是 (x1,y1) 则•ax≡c(mod b) 的解集是：x = k*x1 + t*(b/d) t为任意整数其中模b不同的解共有 d 个，为：x = k*x1 + t*(b/d) t=0,1,..d-1•求ax≡c(mod b) 的最小非负整数解：令 ans = k * x1;s = b/d;则 x = ans + t*s t为任意整数则最小非负整数解是：(ans%s + s)%s•对下面的循环：for (variable = A; variable != B; variable += C)statement;A,B,C和variable都是K(K<=32)比特的无符号数,+=运算的结果也是K比特的无符号数，给定 A,B,C和K，问循环执行多少次。

国际大学生程序设计竞赛辅导教程郭嵩山崔昊吴汉荣陈明睿著北京大学出版社前言ACM 国际大学生程序设计竞赛（ACM International Collegiate Programming Contest，简称ACM/ICPC）是由国际计算机界历史最悠久、颇具权威性的组织ACM 学会（Association for Computer Machinery）主办，是世界上公认的规模最大、水平最高的国际大学生程序设计竞赛，其目的旨在使大学生运用计算机来充分展示自已分析问题和解决问题的能力。

该项竞赛从 1970 年举办至今已历 25 届，因历届竞赛都荟萃了世界各大洲的精英，云集了计算机界的“希望之星”，而受到国际各知名大学的重视，并受到全世界各著名计算机公司的高度关注，成为世界各国大学生最具影响力的国际级计算机类的赛事，ACM 所颁发的获奖证书也为世界各著名计算机公司、各知名大学所认可。

该项竞赛分区域预赛和世界决赛两个阶段进行，各预赛区第一名自动获得参加世界决赛的资格，世界决赛安排在每年的 3~4 月举行，而区域预赛安排在上一年的 9 月~12 月在各大洲举行。

ACM/ICPC 的区域预赛是规模很大，范围很广的赛事，近几年，全世界有 1000 多所大学，近 2000 支参赛队在六大洲的 28~30 个赛站中争夺世界决赛的 60 个名额，其激烈程度可想而知。

与其他编程竞赛相比，ACM/ICPC 题目难度更大，更强调算法的高效性，不仅要解决一个指定的命题，而且必需要以最佳的方式解决指定的命题；它涉及知识面广，与大学计算机系本科以及研究生如程序设计、离散数学、数据结构、人工智能、算法分析与设计等相关课程直接关联，对数学要求更高，由于采用英文命题，对英语要求较高，ACM/ICPC 采用3 人合作、共用一台电脑，所以它更强调团队协作精神；由于许多题目并无现成的算法，需要具备创新的精神，ACM/ICPC 不仅强调学科的基础，更强调全面素质和能力的培养；由于ACM/ICPC 是采用5 小时全封闭式竞赛，参赛队员与外界完全隔离，完全独立完成，没有任何水份，是其实际能力的真实表露，其成绩可信度甚高；但ACM/ICPC 又是一种“开卷考试”，可以带任何书籍、资料甚至源程序代码清单（但不能带软盘），不需要去死背算法，而强调的是算法的灵活运用；与其它计算机竞赛（如软件设计，网站设计等）相比，ACM/ICPC 有严谨而客观的评判规则（严格的数据测试），排除了因评委的主观因素而造成评审不公平的现象，所以，ACM/ICPC 对成绩的争议较少，大家比较心服口服。

等价类
• Ek(等价类)：设G是1…n的置换群。若K是1…n中 某个元素，K在G作用下的轨迹，记作Ek。即K在 G的作用下所能变化成的所有元素的集合。
方案1在四个置换作用 下都是方案1，所以 E1={1}；方案3，在a1 下是3，在a2下变成6， 在a3下变成5，在a4下 变成4，所以 E3={3,4,5,6}；方案11， 在a1、a3下是11，在 a2、a4下变成12，所 以E11={11,12}。
1 x p % p 1 x p % p 显然 : p为质数时，
( 1 x )n ( 1 x ) pn / p ( 1 x )a0 (1 x p )n / p ( 1 x )a0 n / p n / p a0 a0 pi j x x i j i 0 j 0 上式左边的 x m系数一定是右边 x m / p p的系数和 x b0 的系数相乘得到 n n/p a0 于是我们就有了 m m / p b 0
组合数恒等式
n n n n 2 0 1 0 n n n n n n 0 2 4 1 3 5 n n n n 1 1 2 n n 2 1 2 n 0 1 n n 1 k k k k 1
n n (n 1) ... (n m 1) m m!
• 直接求 O(m)
大组合数取模(*)
n,m为大数,p<=10^5,p为质数
n % p m

常用公式 (2)划分问题： (3)Stirling公式 (3)皮克定理 (3)卡特兰数 (4)错排公式 (4)等比数列 (5)等差数列 (5)二次函数 (6)二次方程 (7)约瑟夫环 (7)多边形面积 (7)均值不等式的简介 (8)均值不等式的变形 (8)Lucas 定理 (9)斐波那契数列 (10)欧拉函数 (11)蚂蚁爬绳 (12)(a/b)%m (13)泰勒公式 (13)乘法与因式分解公式 (14)三角不等式 (14)某些数列的前n项和 (15)二项式展开公式 (15)三角函数公式 (16)常用公式划分问题：1、n个点最多把直线分成C(n,0)+C(n,1)份；2、n条直线最多把平面分成C(n,0)+C(n,1)+C(n,2)份；3、n个平面最多把空间分成C(n,0)+C(n,1)+C(n,2)+C(n,3)=(n³+5n+6)/6份；4、n个空间最多把“时空”分成C(n,0)+C(n,1)+C(n,2)+C(n,3)+C(n,4)份. Stirling公式lim(n→∞) √(2πn) * (n/e)^n = n!也就是说当n很大的时候,n!与√(2πn) * (n/e) ^ n的值十分接近这就是Stirling公式.(2πn) ^0.5×n^ n ×e^(-n) =n!;皮克定理一个多边形的顶点如果全是格点，这多边形就叫做格点多边形。

有趣的是，这种格点多边形的面积计算起来很方便，只要数一下图形边线上的点的数目及图内的点的数目，就可用公式算出。

这个公式是皮克(Pick)在1899年给出的，被称为“皮克定理”，这是一个实用而有趣的定理。

给定顶点坐标均是整点（或正方形格点）的简单多边形，皮克定理说明了其面积S和内部格点数目a、边上格点数目b的关系：S=a+ b/2 - 1。

(其中a表示多边形内部的点数,b表示多边形边界上的点数,S表示多边形的面积) 卡特兰数原理:令h(1)＝1,h(0)=1，catalan数满足递归式：h(n)= h(0)*h(n-1)+h(1)*h(n-2) + ... + h(n-1)h(0) (其中n>=2)另类递归式：h(n) = h(n-1)*(4*n-2)/(n+1);该递推关系的解为：h(n)=C(2n,n)/(n+1) (n=1,2,3,...)卡特兰数的应用（实质上都是递归等式的应用）错排公式当n个编号元素放在n个编号位置,元素编号与位置编号各不对应的方法数用M(n)表示,那么M(n-1)就表示n-1个编号元素放在n-1个编号位置,各不对应的方法数,其它类推.第一步,把第n个元素放在一个位置,比如位置k,一共有n-1种方法;第二步,放编号为k的元素,这时有两种情况.1,把它放到位置n,那么,对于剩下的n -2个元素,就有M(n-2)种方法;2,不把它放到位置n,这时,对于这n-1个元素,有M(n-1)种方法;综上得到递推公式：M(n)=(n-1)[M(n-2)+M(n-1)] 特殊地，M(1)=0,M(2)=1通项公式：M(n)=n![(-1)^2/2!+…+(-1)^(n-1)/(n-1)!+(-1)^n/n!]优美的式子：Dn=[n!/e+0.5],[x]为取整函数．公式证明较简单．观察一般书上的公式，可以发现e-1的前项与之相同，然后作比较可得/Dn-n!e-1/<1/(n+1)<0.5,于是就得到这个简单而优美的公式（此仅供参考）等比数列(1) 等比数列：a (n+1)/an=q (n∈N)。

常见置换的循环数
计算置换的循环数,是这一算法的瓶颈.如果能 够快速计算出各置换的循环数,就可以大大提 高程序的运行效率 旋转:n个点顺时针(或逆时针)旋转i个位置的置 换,循环数为gcd(n,i) 翻转:
n为偶数时,
对称轴不过顶点:循环数为n/2 对称轴过顶点:循环数为n/2+1
置换群
以置换为元素的群 置换群G={a1,a2,...,a|G|} 例子中G内共有6个置换 123456 123456 123456 123456 612345 561234 123456 123456 123456 456123 345612 234561
循环
在一个置换下,x1->x2,x2->x3,...,xn->x1,这 样x1,x2,...,xn就构成了一个循环 定义ck为在置换ak下的循环总数 例子中: c1=6,c2=1,c3=2,c4=3,c5=2,c6=1
Polya定理
设G={a1,a2,...,a|G|}是N={1,2,...,N}上的置 换群,现用m种颜色对这N个点染色,则不 同的染色方案数为 S=(mc1+mc2+...+mc|G|)/|G| 证明比较复杂,略
利用Polya定理 解决组合计数问题的步骤
写出置换群 123456 123456 123456 123456 612345 561234 123456 123456 123456 456123 345612 234561 求出每个置换的循环数 c1=6,c2=1,c3=2,c4=3,c5=2,c6=1 计算染色方案 S=(26+21+22+23+22+21)/6=14
解二元模线性方程

22
for (w=G->first(v); w<G->n(); w = G>next(v,w))
if (G->getMark(w) == UNVISITED) { G->setMark(w, VISITED); Q->enqueue(w); }
PostVisit(G, v); // Take action }}
Q Q
31
() ((1,1)) ((1,1) (2,3)) ((1,1) (2,4)) ((1,1) (2,4) (3.2))
Q Q
Q
32
() ((1,1)) ((1,1) (2,3)) ((1,1) (2,4)) ((1,1) (2,4) (3.2))
Q Q
33
Q () ((1,1)) ((1,1) (2,3)) ((1,1) (2,4)) ((1,1) (2,4) (3.2))
求解算法
穷举法 贪心法 深度/广度优先搜索
图和树的搜索 回溯法（深度优先）
递归
分治
动态规划法
55
谢谢！
56
49
动态规划法的思想
自底向上构造 在原问题的小子集中计算 每一步列出局部最优解，并用一张
表保留这些局部最优解，以避免重 复计算 子集越来越大，最终在问题的全集 上计算，所得出的就是整体最优 解。
50
Fib(1)=Fib(2)=1
Fibonacci数列 Fib(n)=Fib(n-1) + Fib (n-2)
23
八数码 深度优先搜索
1
23 184 765
2
c
3
283 164 75
283 14 765

57 2020/4/15
58 2020/4/15
59 2020/4/15
60 2020/4/15
61 2020/4/15
62 2020/4/15
63 2020/4/15
64 2020/4/15
特别提醒：
以上介绍的线段的三个属性， 是计算几何的基础，在很多方 面都有应用，比如求凸包等等， 请务必掌握！
15 2020/4/15
第二单元
多边形面积 和重心
16 2020/4/15
基本问题（1）：
给定一个简单多边形，求其 面积。
输入：多边形（顶点按逆时 针顺序排列）
输出：面积S
17 2020/4/15
A=sigma(Ai) (i=1…N-2)
P1
A1 P2
P6 A4
P5 A3
A2 P4
P3
25 2020/4/15
凹多边形的面积？
P3
P2 P4
P1
26 2020/4/15
依然成立！！！
多边形面积公式：A=sigma(Ai) (i=1…N-2)
结论： “有向面积”A比“面积”S其实更本
质！
27 2020/4/15
思考如下图形：
18 2020/4/15
Any good idea?
19 2020/4/15
先讨论最简单的多边形——三角形
20 2020/4/15
三角形的面积：
在解析几何里， △ABC的面积可以通过 如下方法求得：
点坐标 => 边长 => 海伦公式 => 面积
21 2020/4/15
思考：此方法的缺点：
C=sigma((↑Pi +↑Pi+1)(↑Pi ×↑Pi+1) ) /(6A)