ACM模板

  1 /*
  2  * text.cpp
  3  *
  4  *  Created on: 2014-8-20
  5  *      Author: Administrator
  6  */
  7 
  8 #include <stdio.h>
  9 #include <iostream>
 10 #include <algorithm>
 11 #include <math.h>
 12 #include <string.h>
 13 #include <vector>
 14 #include <queue>
 15 using namespace std;
 16 #define N 1010
 17 #define LL __int64
 18 #define INF (int)1e7
 19 
 20 /* 2.1   基本计数方法：
 21  *
 22  * 1.容斥原理：集合之间有重复的部分，加加减减：大意为：A、B集合的相加减去A、B集合重叠的部分。
 23  * 2.组合性质1：组合数的递推公式：经常用于预处理：C(n+1,k+1) = C(n,k+1) +C(n,k);  （补充不漏的递推）
 24  * 3.组合性质2：求n一定时，k从1-n的所有组合数，O(n)的递推：C(n,k+1) = C(n,k)*(n-k)/(k+1); （注意不要溢出）
 25  * 4.组合性质3：可重复选择的组合，n个不同的数，每个元素可以选多多次，一共选k个：问题转化：x1 + x2 + x3.. xn = k 的解的个数，
 26  *                 答案: C(k+n-1,n-1) = C(k+n-1,k);
 27  * 4.排列性质1：有重复数字的全排列个数：设x为结果，则：n1! * n2! *... nk! * x = n! （n[i]为每种元素重复的个数）
 28  * 5.单色三角形的反向考虑：枚举点，每个点上连接a条红边和n-1-a条黑边，这些边属于a(n-1-a)个三角形！每个考虑了两次，最后除以2，累加即可。
 29  * 6.swap(a,b); 可以交换两个数字
 30  * 7.定义long long ，之后使用cout输出即可。
 31  * 在题目中要注意，加法原理，乘法原理的应用，将问题分解，逐部分去求解。
 32  * 有的时候，看了题目，明显发现时间复杂度过高，一般都是需要进行一些数学分析。
 33  */
 34 
 35 /* 2.2   递推关系
 36  *
 37  * 一般在n范围较大，又设计树、图，或者说给定一个很大的范围，要计算其中每一个数字是否满足某种情况，考虑递推降低复杂度。
 38  */
 39 
 40 /* 2.3   数论
 41  *
 42  * 基本概念：
 43  * 1.素数筛选法
 44  * 2.欧几里得算法：计算最大公约数 && ax + by = d的 |x|+|y|最小
 45  * 3.取模算法：乘积取模与幂取模
 46  * 4.欧拉phi函数：返回不超过x且和x互素的正整数个数（一个计算单个，一个筛选法计算数组）
 47  * 5.素数筛选 + 欧拉函数 的应用：最大公约数之和。
 48  *             --------------技巧：递推 + 转换 + 逆转思路筛选(枚举n的倍数可以很大程度降低复杂度，例如暑期训练第九期的1002题)
 49  *
 50  * 模方程：
 51  * 1.线性模方程
 52  * 2.中国剩余定理等。    （还没看）
 53  *
 54  */
 55 
 56 //======素数筛选法======
 57 
 58 
 59 //======欧几里得算法======
 60 //原始算法只能求得两个整数a、b的最大公约数
 61 LL gcd(LL a,LL b){
 62     return b==0 ? a : gcd(b,a%b);
 63 }
 64 //扩展的还可以求出x，y，使得 ax + by = d，|x|+|y| 最小。
 65 //注意，即使a,b在int范围内，x，y也有可能超出int范围
 66 void gcd_two(LL a,LL b,LL& d,LL& x,LL& y){
 67     if(!b){ d = a; x = 1; y = 0;}
 68     else{ gcd_two(b,a%b,d,y,x); y -= x*(a/b);}
 69 }
 70 
 71 //======取模运算======
 72 //返回ab mod n 要求： 0 <= a,b <n
 73 LL mul_mod(LL a ,LL b,int n){
 74     return a*b%n;
 75 }
 76 //返回a^p mod n,要求：0 <= a < n  (进行了降幂处理 log(p))
 77 LL pow_mod(LL a,LL p,LL n){
 78     if(p == 0) return 1;
 79     LL ans = pow_mod(a,p/2,n);
 80     ans = ans*ans%n;
 81     if(p%2 == 1) ans = ans*a%n;
 82     return ans;
 83 }
 84 //======欧拉phi函数======
 85 //计算某个x的欧拉phi函数，不超过n且与n互素的正整数个数
 86 int euler_phi(int n){
 87     int m = (int)sqrt(n+0.5);
 88     int ans = n;
 89     for(int i=2;i<=m;i++) if(n%i==0){
 90         ans = ans/i*(i-1);
 91         while(n%i==0) n/=i;
 92     }
 93     if(n>1) ans=ans/n*(n-1);
 94     return ans;
 95 }
 96 //使用类似筛选法的方法计算phi(1),phi(2)...
 97 int phi[200];
 98 void phi_table(int n){
 99     for(int i=2;i<=n;i++) phi[i]=0;
100     phi[1]=1;
101     for(int i=2;i<=n;i++) if(!phi[i])
102         for(int j=i;j<=n;j+=i){
103             if(!phi[j]) phi[j]=j;
104             phi[j] = phi[j]/i*(i-1);
105         }
106 }
107 
108 /* 2.4   组合游戏
109  *
110  * 这里大多指：博弈游戏
111  * 规则：1.一个状态是必败态当且仅当它所有的后继都是必胜状态；
112  *          2.一个状态是必胜态当且仅当它至少有一个后继是必败状态。
113  * --- 有了这个方法，我们可以用递推的方法，判断整个状态图的每一个节点是必胜态还是必败态。（此状态图无环，按照拓扑排序的逆序进行递推判断）
114  * 1.Nim游戏：取火柴，一共三堆，判断先手胜还是后手胜。
115  * 2.Ferguson游戏：两个盒子：清空一个，将另一个的分到两个盒子中，每个盒子至少有一个，不能继续的败！-- 后有代码
116  * 3.Chomp!游戏：有一个m*n的棋盘，每次可以取走一个方格并拿掉它右边和上边的所有盒子，拿到(1,1)的算输。
117  *         --- 除了(1,1)是先手必败之外，其余都是先手必胜；考虑到先手在第一次取方格的时候，都可以抢先一步后手进入必胜态。
118  * 4.约数游戏：有1~n个数字，两人轮流选，并把它和它所有的约数擦去。擦去最后一个数的人赢。
119  *         --- 补充Nim游戏的定理： 状态(x1,x2...xn)为必败态，当且仅当所有xi的异或结果为0.（适用于任意堆的情况）
120  * 5.组合游戏的和：使用SG函数和SG定理
121  *             --- 137页，回头看
122  *
123  */
124 
125 //Ferguson游戏：输出所有的k<t的必败状态
126 const int maxn = 100;
127 int winning[maxn][maxn];
128 void Ferguson(int t){
129     winning[1][1]=false;
130     for(int k=3;k<20;k++)
131         for(int n=1;n<k;n++){
132             int m = k-n;
133             winning[n][m] = false;
134             for(int i=1;i<n;i++)  //后继至少有一个必败态的时候， 此时为必胜态
135                 if(!winning[i][n-i]) winning[n][m]=true;
136             for(int i=1;i<m;i++)  //后继至少有一个必败态的时候， 此时为必胜态
137                 if(!winning[i][m-i]) winning[n][m]=true;
138             if(n<=m && !winning[n][m]) printf("%d %d\n",n,m);
139         }
140 }
141 
142 
143 
144 
145 /* 2.5   概率和数学期望 --- 一般需要开设数组进行计算，包括一维 + 二维 ；有一个特点，如果建图，一般都是拓扑结构
146  *
147  * 1.全概率公式（相当于平均概率）：关键是“划分样本空间”，只有把所有情况不重复不遗漏地进行分类，并计算出每个分类下事件
148  * 2.数学期望：所有可能值按照概率加权之后的结果
149  *         --- 过程： 一般先求出所有的可能值，之后求出相应值得概率，之后加权求和。
150  *         --- 两个求解工具：
151  *             (1)期望的线性性质： E(x+y) = E(x) + E(y)
152  *             (2)全期望公式：类似全概率公式，将所有情况不重复、不遗漏地分成若干类，每类计算数学期望，之后将每类的概率加权求和。
153  * 3.例题：麻球繁衍，用到求期望的独立性质，每个麻球的死亡是独立的，k天后，所有都死亡的概率则为其的k次方。
154  * 4.例题：和朋友会面，相当于区间概率的问题，建模改为二维平面围成的面积。（但是需要分类讨论，强调二维建模）
155  * 5.例题：玩纸牌（142页）（强调独立性）
156  *         --- (1)每晚输赢是相互独立的，先计算出单独一天垂头丧气去睡觉的概率 Q
157  *        --- (2)设d(i,j)表示前i局每局结束后的获胜比例不超过p，且前i局一共获胜j局的概率，
158  *                     根据全概率公式：j/i <= p 时，d(i,j) = d(i-1,j)*(1-p) + d(i-1,j-1)*p;
159  *        --- (3)则大Q = d(n,0) + d(n,1) +... 。
160  *        --- (4)得到每天垂头丧气去睡觉的概率之后，这里还有另外一个结论：游戏总天数的期望即为： 1/Q
161  * 6.例题：得到1： 本题可以看做一个随机状态的状态机（这里强调一个特点随机性）
162  *        --- 这个使用的记忆化搜索，但还不太会。
163  *
164  */
165 
166 
167 
168 /* 2.6   置换及其应用：这里指一一映射的一种关系。
169  *
170  * 1.置换满足结合律，但是不满足交换律，其置换式遵循一定顺序的，不能交换。
171  * 2.置换分解显示在图中，会表示成多个有向圈，其中循环的个数我们称为该置换的循环节
172  * 3.一个经典问题：等价类计数问题：即题目中会在已有情况的基础上，会给我们定义一种等价的情况。
173  *         --- 一般这类问题，都需要统计等价类的个数，再延伸出来一些应用。
174  *        --- Burnside引理：用来统计等价类个数：
175  *                 对于一个置换f，如果有一个方案s经过置换后不变，则s称为f的一个不动点，将f的不动点数目记为C(f)
176  *         --- 如何求C(f)? ,
177  * 4.例题：项链和手镯：输出使用t种颜色的n颗珠子能做成的项链或者手镯的个数。（手镯可以旋转，隐藏有等价类）
178  *         --- 。。。
179  * 5.例题：Leonardo笔记本：给出26个大写字母的置换B，问是否存在一个置换A，使得A^2 = B
180  *         。。。还有两个例题
181  *
182  */
183 
184 
185 /* 2.7   矩阵和线性方程组
186  *
187  * 1.矩阵运算满足结合律!多校合练中有一道(A*B)^n的矩阵想成，就用到了此性质，将矩阵行列数目变得很小，再使用矩阵快速幂解决问题。
188  * 2.逆矩阵：几乎不使用，因为矩阵求逆大多数情况下是可以避免的；
189  * 3.高斯消元法计算举例，解方程组 --- 用到增广矩阵（后有代码）
190  * 4.例题：递推关系
191  * 5.例题：细胞自动机
192  *        ...
193  *
194  */
195 
196 //求解多项式的每个未知变量的唯一解
197 //要求系数矩阵可逆 && A是增广矩阵， A[i][n] 是第i个方程右边的常数bi
198 //运行结束后 A[i][n]是第i个未知数的解
199 #define maxn 1010
200 void gauss_elimination(double A[maxn][maxn],int n){
201     int i,j,k,r;
202     for(i=0;i<n;i++){
203         //选一行r并与第i行交换
204         r = i;
205         for(j=i+1;j<n;j++)
206             if(fabs(A[j][i]) > fabs(A[r][i])) r=j;
207         //与第i+1 ~ n行进行消元e
208         for(j=n;j>=i;j--)
209             for(k=i+1;k<n;++k){
210                 A[k][j]-=A[k][i]/A[i][i]*A[i][j]; //这里消元使用的高精度要求的方法。
211             }
212     }
213     //回带过程
214     for(i =n-1;i>=0;i--){
215         for(j=i+1;j<n;j++)
216             A[i][n] -= A[j][n]*A[i][j];
217         A[i][n] /= A[i][i];
218     }
219 }
220 
221 //整数快速幂模板
222 LL quickpow(LL m ,LL n ,LL k){
223     LL ans = 1;
224     while(n){
225         if(n&1)//如果n是奇数
226             ans = (ans * m) % k;
227         n = n >> 1;//位运算“右移1类似除2”
228         m = (m * m) % k;
229     }
230     return ans;
231 }
232 
233 //矩阵快速幂模板（迭代法，还有一种使用递归实现）
234 
235 
236 
237 /* 2.8   数值方法简介
238  *
239  * 1.大意是说使用数值的方法算出近似解。常用：非线性方程求根 + 凸函数求极值 + 数值积分
240  * 2.例题：解方程
241  * 3.例题：误差曲线
242  * 4.例题：桥上的绳索：使用辛普森公式求解积分。（二维图上，图像与坐标轴围成的面积）
243  *
244  */
245 
246 
247 
248 /* 3.1   基础数据结构回顾
249  *
250  * 1.优先队列
251  * 2.并查集
252  */
253 
254 /* 3.2   区间信息的维护与查询：连续和查询问题
255  * 1.二叉搜索树(树状数组)：动态和查询问题
256  *         --- 标准用法：我们需要动态地修改单个元素值并求前缀和。
257  * 2.RMQ问题：范围最小值问题：预处理n*log(n),每次查询O(1)
258  *         --- 标准用法：给定一个区间(L,R);求min(A[L]...A[R]);
259  *         --- 使用数据结构： d(i,j)二维数组,表示从i开始，长度为2^j的一段元素中的最小值,可以使用递推的方法来计算d(i,j),d(i,j)=min}{d(i,j-1),d(i+2^j-1,j-1)};
260  * 3.线段树：点修改 + 区间修改 --- 动态范围的最小值问题。（RMQ不提供动态操作，只能是预处理不变的数组）
261  *
262  */
263 
264 
265 //树状数组：C[]为主要数据结构
266 //初始化使用方法： 先预处理a[]和c[]数组为0，之后执行n此add操作即可。
267 //a[]为长度为n的数组，add时将a[]中的数字依次add进去。（记得0号位不使用!）
268 int n; //表示节点个数
269 int c[105];
270 int lowbit(int x){
271     return x&(-x);
272 }
273 int sum(int x){
274     int ret =0;
275     while(x>0) {ret += c[x]; x-= lowbit(x);}
276     return ret;
277 }
278 void add(int x,int d){
279     while(x<=n){
280         c[x] += d;x+=lowbit(x);
281     }
282 }
283 void init(int A[],int n){
284     memset(c,0,sizeof(c));
285     for(int i=0;i<n;i++)
286         add(i+1,A[i]);
287 }
288 
289 //RMQ问题：查询L~R之间的最小值
290 //A表示一个n个元素的数组
291 int d[101][101];
292 void RMQ_init(const vector<int>& A){
293     int n = A.size();
294     for(int i=0;i<n;i++) d[i][0]=A[i];
295     for(int j=1;(1<<j) <=n ;j++)
296         for(int i=0;i+(1<<j)-1<n;i++)
297             d[i][j]=min(d[i][j-1],d[i+(1<<(j-1))][j-1]);
298 }
299 //查询操作:返回值即为最小值的数值
300 int RMQ(int L,int R){
301     int k=0;
302     while((1<<(k+1)) <= R-L+1) k++; //如果2^(k+1) <= R-L+1,那么k还可以加一;
303     return min(d[L][k],d[R-(1<<k)+1][k]);
304 }//这样就在O(n*log(n))的预处理之后，做到了O(1)的查询
305 
306 //线段树
307 
308 
309 
310 /* 3.3   字符串(1)
311  *
312  * 1.Trie(前缀树,也叫字典树)：是一种用于快速检索的多叉树结构，如英文字母的字典树是一个26叉树，数字的字典树是一个10叉树。
313  * 2.Trie的应用：例题：背单词，讲一个长字符串分解为已有的字典中的单词(这些单词建成了字典树) + 递推解决。
314  * 3.Trie的应用：例题：strcmp(char *s,char *t)：将所有的字符串构建成Trie树，之后分析情况，按照节点计算需要进行比较的次数。
315  * 4.KMP算法：字符串进行匹配。文本长度为n，匹配串长度为m，在n中寻找子串m。（满足T[i]=p[0] && T[i+1]=p[1] &&.. 连续的相等）
316  * 5.Aho-Corasick自动机：当待匹配串（短的）有多个的时候KMP算法就不是很合适了，因为每次匹配都需要去遍历整个文本串。
317  *         这里使用KMP + Trie树结合，完成AC自动机的算法。
318  * 6.AC自动机的应用：例题：出现次数最多的子串
319  *
320  */
321 
322 //Trie树，用于多字符串(或者多字符数串)，快速查找。
323 const int maxnode=105; //每个字符串的最大长度，解释在Trie树种，表示最大的深度
324 const int sigma_size=26; //表示几叉树，如26个字母组成的字符串，则其为26
325 struct Trie{
326     int ch[maxnode][sigma_size];
327     int val[maxnode];
328     int sz;  //节点总数
329     Trie(){sz=1;memset(ch[0],0,sizeof(ch[0]));}  //初始时只有一个根节点
330     int idx(char c){return c-'a';} //字符c的编号
331 
332     //插入字符串s，附加信息为v ---- 注意v必须非0，因为0代表"本节点不是单词节点"
333     void insert(char *s,int v){
334         int u=0,n=strlen(s);
335         for(int i=0;i<n;i++){
336             int c = idx(s[i]);
337             if(!ch[u][c]){
338                 memset(ch[sz],0,sizeof(ch[sz]));
339                 val[sz]=0;
340                 ch[u][c]=sz++;
341             }
342             u=ch[u][c];
343         }
344         val[u]=v;
345     }
346     //查询和插入类似
347     int serch(){
348         //待完善
349 
350         return 0;
351     }
352 };
353 
354 //KMP算法：先调用getFail(主串,t); 再调用find(主串,模式串,t);
355 //当有多个匹配的时候，因为没有break，会打印多个匹配项。
356 void getFail(char *P,int *f){
357     int m = strlen(P);
358     f[0]=0;f[1]=0;
359     for(int i=1;i<m;i++){
360         int j=f[i];
361         while(j&&P[i]!=P[j]) j = f[j];
362         f[i+1] = P[i] == P[j] ? j+1 : 0;
363     }
364 }
365 void find(char *T,char *P,int *f){
366     int n =strlen(T),m = strlen(P);
367     getFail(P,f);
368     int j=0;
369     for(int i=0;i<n;i++){
370         while(j&&P[j]!=T[i]) j=f[j];
371         if(P[j]==T[i]) j++;
372         if(j==m) printf("%d\n",i-m+1); //找到了，打印下坐标，这里当做return返回值也可
373     }
374 }
375 
376 /* 3.4   字符串(2)
377  *
378  * 1.后缀数组：AC自动机中，需要提前知道所有的模板串（短的），但是实际中常常不能满足，所以这里延伸出后缀Trie，预先处理文本串T，而不是模板。
379  *         --- 构造方法：将文本串的所有后缀都插入到一个后缀Trie树种，查找一个长度为m的模板的时候只需要进行一次普通的Trie查找，时间复杂度为O(m)
380  *        --- 但是实际中，后缀Trie较为难懂且易错，常常使用后缀数组代替!
381  *        --- 后缀数组中，字符串从0开始，且返回匹配下标；在实际中为叙述方便，直接把"以下标k开头的后缀"叫做"后缀k"。
382  *        --- 得到后缀数组之后，可以进行在线的多模板匹配问题。直接在后缀数组里进行二分查找即可 （221页代码）
383  * 2.最长公共前缀（LCP）
384  *
385  */
386 
387 /* 3.5   排序二叉树：每个节点可以保存着一个可以比较大小的东西，如：整数、字符串等。
388  *
389  * 1.排序二叉树用来实现STL中的set等，不常用。（有时候直接使用STL中的县城库更好一些）
390  * 2.用Treap实现名次数 --- 不会
391  *         。。。
392  *
393  */
394 
395 
396 
397 /* 4.1二维几何基础
398  *
399  * 1.二维几何的数据结构定义，及常用运算的重载。包括+-*\/ < ==
400  *         --- 基本运算：点积,等于两个向量v和w的的长度乘积*它们夹角的余弦。
401  *                点积计算 + 利用点积计算向量长度和夹角的函数
402  *         --- 基本运算：叉积：以及利用叉积计算有向面积
403  *         --- 两个向量的位置关系：使用点积与叉积的符号放在一起进行判断。（256页的图4-4）
404  *         --- 向量旋转：向量可以绕起点旋转： x' = xcosa - ysina, y' = xsina + ycosa; a为逆时针旋转的角度，弧度制
405  *         --- 基于复数的集合计算：略
406  * 2.点和直线：
407  *         --- 直线的参数表示：直线可以使用直线上一点P0和方向向量v表示； P = P0 + tv ;对于直线、射线、线段的不同，仅限于t的取值范围依次减小
408  *         --- 直线交点：
409  *         --- 点到线的距离： 点到直线 + 点到线段
410  *         --- 点在直线上的投影点
411  *         --- 线段相交判定
412  *         --- 判断点是否在线段上（不包括端点）
413  * 3.多边形：
414  *         --- 计算多边形的有向面积：凹凸均可，因为这里计算有向面积，正负会相互抵消。（将多边形划分成多个三角形，可以从任意顶点出发进行划分计算!）
415  *        --- 欧拉定理： 设平面图的定点数、边数和面数分别为： V,E和F，则V+F-E = 2;
416  */
417 /* 4.2   与圆和球有关的计算问题：
418  *
419  * 1.定义圆的结构体 + 通过圆心角求点在圆上坐标的函数
420  * 2.直线与圆的交点：求圆心到直线的距离，之后与半径r做比较来判断；
421  * 3.两圆相交代码： 增加一个计算向量极角的方法： 传入参数：两个圆 + 一个Vector<Point> 记录交点坐标
422  * 4.过定点求圆的切线：
423  * 5.两圆的公切线：
424  * 6.外公切线：
425  */
426 //向量和点的定义，相同，但是起了不同的别名
427 struct Point{
428     double x,y;
429     Point(double x=0,double y=0):x(x),y(y){} //构造函数，方便代码的编写
430 };
431 typedef Point Vector; //从程序上，Vector只是Point的别名
432 
433 //向量 + 向量 = 向量 , 点 + 向量 = 点
434 Vector operator + (Vector A, Vector B){return Vector(A.x+B.x,A.y+B.y);}
435 //点 - 点 = 向量
436 Vector operator - (Point A, Point B){return Vector(A.x-B.x,A.y-B.y);}
437 //向量 * 数 = 向量
438 Vector operator * (Vector A, double p){return Vector(A.x*p,A.y*p);}
439 //向量 / 数 = 向量
440 Vector operator / (Vector A, double p){return Vector(A.x/p,A.y/p);}
441 bool operator < (const Point& a,const Point& b){
442     return a.x < b.x || (a.x == b.x && a.y<b.y);
443 }
444 const double eps = 1e-10;
445 int dcmp(double x){
446     if(fabs(x) <eps) return 0;
447     else return x < 0 ? -1 : 1;
448 }
449 bool operator == (const Point& a,const Point& b){
450     return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y)==0 ;
451 }
452 
453 
454 //点积：以及计算向量长度和其夹角
455 double Dot(Vector A,Vector B){  //返回两向量的点积
456     return A.x*B.x + A.y*B.y;
457 }
458 double Length(Vector A){
459     return sqrt(Dot(A,A));
460 }
461 double Angle(Vector A,Vector B){ //返回两向量夹角
462     return acos(Dot(A,B)/Length(A)/Length(B));
463 }
464 
465 //叉积：以及计算有向面积
466 double Cross(Vector A,Vector B){
467     return A.x*B.y - A.y*B.x;
468 }
469 double Area2(Point A,Point B,Point C){  //返回两向量的有向面积，第二个向量在第一个向量左边返回正值，否则返回负值
470     return Cross(B-A,C-A);  // 是三角形有向面积的两倍
471 }
472 
473 //向量旋转：
474 Vector Rotate(Vector A,double rad){
475     return Vector(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad) + A.y*cos(rad));
476 }
477 //作为特殊情况，下面的函数计算向量的单位法线，即左转90°之后再把长度归一化的 一个单位法向量
478 Vector Normal(Vector A){
479     double L = Length(A);
480     return Vector(-A.y/L,A.x/L);
481 }
482 
483 //==============点和直线================
484 //直线交点：直线分别为 P + tv 和 Q +tw ，
485 Point GetLineIntersection(Point P,Vector v,Point Q, Vector w){
486     Vector u = P-Q;
487     double t = Cross(w,u)/Cross(v,w);
488     return P+v*t; //返回交点
489 }
490 //点到直线的距离：
491 double DistanceToLine(Point P,Point A,Point B){
492     Vector v1 = B - A,v2 = P - A;
493     return fabs(Cross(v1,v2))/Length(v1);  //如果不取绝对值，得到的是有向距离
494 }
495 //点到线段的距离：
496 double DistanceToSegment(Point P,Point A,Point B){
497     if(A==B) return Length(P-A);
498     Vector v1 = B-A,v2=P-A,v3=P-B;
499     if(dcmp(Dot(v1,v2))<0) return Length(v2);
500     else if(dcmp(Dot(v1,v3))>0) return Length(v3);
501     else return fabs(Cross(v1,v2))/Length(v1);
502 }
503 //点在直线上的投影点
504 Point GetLineProjection(Point P,Point A,Point B){
505     Vector v = B-A;
506     return A+v*(Dot(v,P-A)/Dot(v,v));
507 }
508 //线段判定相交：
509 //规范相交：两个线段有且仅有一个交点，且不在任何一条线段的端点
510 bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2){
511     double c1 = Cross(a2-a1,b1-a1),c2 = Cross(a2-a1,b2-a1),
512             c3 = Cross(b2-b1,a1-b1),c4 = Cross(b2-b1,a2-b1);
513     return dcmp(c1)*dcmp(c2)<0 && dcmp(c3)*dcmp(c4)<0 ;
514 }
515 //判断一个点是否在一条线段（不包含端点）上：点p，线段 a1-a2;
516 bool OnSegment(Point p,Point a1,Point a2){
517     return dcmp(Cross(a1-p,a2-p))==0 && dcmp(Dot(a1-p,a2-p))<0 ;
518 }
519 
520 //===============多边形=================
521 //计算多边形面积：
522 double ConvexPolygonArea(Point* p,int n){
523     //传进来顶点的数组，与顶点的个数
524     double area = 0.0;
525     for(int i=1;i<n-1;i++)
526         area += Cross(p[i]-p[0],p[i+1]-p[0]); //这里讲p[0]当做了分割点
527     return area/2;
528 }
529 struct Circle{
530     Point c;
531     double r;
532     Circle(Point c,double r):c(c),r(r){}
533     Point point(double a){ //圆上任意一点拥有唯一的圆心角，这里加一个通过圆心角求点在圆上坐标的函数
534         return Point(c.x + cos(a)*r, c.y + sin(a)*r);
535     }
536 };
537 //两圆相交 + 传入两圆，得交点坐标
538 //计算向量极角
539 double angle(Vector v){
540     return atan2(v.y,v.x);
541 }
542 int getCircleCircleIntersection(Circle c1,Circle c2,vector<Point>& sol){
543     double d = Length(c1.c - c2.c);
544     if(dcmp(d)==0){
545         if(dcmp(c1.r - c2.r) == 0) return -1; //两圆重合
546         return 0;
547     }
548     if(dcmp(c1.r + c2.r - d)<0) return 0;
549     if(dcmp(fabs(c1.r - c2.r)-d)>0) return 0;
550 
551     double a = angle(c2.c - c1.c); //向量c1c2的极角
552     double da = acos((c1.r*c1.r + d*d -c2.r*c2.r)/(2*c1.r*d));
553     //c1c2到c1p1的角
554     Point p1 = c1.point(a-da),p2 = c1.point(a+da);
555 
556     sol.push_back(p1);
557     if(p1 == p2) return 1;
558     sol.push_back(p2);
559     return 2;
560 }
561 
562 /* 4.3   二维几何常用算法：
563  *
564  * 1.点在多边形内的判定：这里使用转角法：0°在外、180°在边上、360°在内（这个多边形甚至可以自交）--- 271页
565  * 2.获得凸包：就是把给定点包围在内部的，面积最小的凸多边形。  N*log(N)
566  * 3.半平面交问题：相交结果一般是一个凸多边形 or 无界多边形 or 一条直线、线段、一个点、
567  * 4.
568  *
569  */
570 
571 //(1)先将坐标排序：x从小到大，相同时，y从小到大，并且删除重复点。
572 //计算凸包，输入点数组p，个数为n，输出点个数为ch。 函数返回凸包顶点个数；
573 //输入不能有重复点
574 //如果不希望在凸包的边上有点，把两个<= 换成 <
575 //在京都要求高的时候，建议使用dcmp()
576 int cmd(Point a,Point b){
577     if(fabs(a.x-b.x)<= 1e-9) return a.y<b.y;
578     else return a.x<b.x;
579 }
580 int ConvexHull(Point *p,int n,Point* ch){
581     sort(p,p+n,cmp); //定义排序函数
582     int m = 0;
583     for(int i=0;i<n;i++){
584         while(m>1 && Cross(ch[m-1] - ch[m-2],p[i] - ch[m-2]) <= 0) m--;
585         ch[m++] = p[i];
586     }
587     int k =m;
588     for(int i=n-2;i>=0;i--){
589         while(m>k && Cross(ch[m-1] - ch[m-2],p[i] - ch[m-2]) <= 0) m--;
590         ch[m++] = p[i];
591     }
592     if(n>1) m--;
593     return m; //返回凸包顶点个数
594 }
595 
596 //每个半平面使用一条有向直线表示： 276页
597 struct Line{
598     Point p;
599     Vector v;
600     double ang;
601     Line(){}
602     Line(Point p,Vector v):p(p),v(v){ang = atan2(v.y,v.x);}
603     bool operator < (const Line& L) const{
604         return ang < L.ang;
605     }
606 };
607 
608 
609 //=====================第5章：图论算法与模型========================
610 //无根树转换为有根数 + 表达式求值 + 最小生成树 + 最短路 + 最大流 + 最小割 + 最小费用流等 --- 充分展示 建模和分析 的技巧
611 
612 /* 5.1   基础题目选讲：图的宽度优先遍历 + 欧拉回路 + 拓扑排序
613  *
614  * 1.例题：大火蔓延的迷宫：宽度优先遍历 + 最短路
615  *         --- 输出走出迷宫的最短时间； 加上火只需要预处理每个格子着火的时间，在判断的时候，将此点作为无法通过的点即可。
616  *         --- 如何求出每个格子什么时间会起火，是一个最短路的问题（只不过一点不同是起点变成了多个，在初始化队列的时候将所有着火点都加入进去即可）
617  * 2.例题：独轮车：
618  *         --- 此题多了两个附加因素：一个车轮颜色、一个是到达格子时的朝向
619  * 3.例题：项链：彩色珠子，每颗珠子左右两部分颜色不同，相连时必须在接触的地方颜色相同，问是否可以复原成完整的项链
620  *         --- 这题充分体现了建模的重要性；
621  *         --- 将每个颜色当做一个点，左右两边两种颜色，也即两点之间加上一条边，之后判断全图是否存在一条欧拉回路即可！
622  *         --- 若图G中存在这样一条路径，使得它恰通过G中每条边一次,则称该路径为欧拉路径。若该路径是一个圈，则称为欧拉(Euler)回路。
623  * 4.例题：猜序列：连续和转换为前缀和之差。
624  *         --- 本题可以转换为： 0,b1,b2...bn(前缀和)的一些大小关系，求它们的值，只需要一组解即可。
625  *         --- 此题可以通过拓扑排序完成
626  *
627  */
628 
629 /* 5.2   深度优先遍历：  dfs + 连通分量 + 二分图判定
630  *
631  * 1.dfs框架：依次递归访问当前结点的所有相邻结点
632  * 2.连通分量：在图中，我们将相互可达的结点称为一个连通分量  ---  如果只需要计算连通分量，还可以使用并查集，效果更好，甚至图都可以不需要保存
633  * 3.二分图判定：非连通的图是二分图当且仅当每个连通分量都是二分图，且这里只考虑无向连通图
634  * 4.无向图的割顶和桥：
635  *         --- 割顶：对于非连通图，如果删除某个点后连通分量数目增加，则称u为图的关节点，即割顶；对于连通图，割顶就是删除之后使图不再连通的点。
636  *         --- 桥：删除一条边使图非连通的边。
637  *         --- 不知道可以应用在什么地方，这里没总结。
638  * 5.无向图的双连通分量：315页
639  *         --- 内部无割点：此图是点-双连通的，任意两条边都在一个简单的环中；
640  *         --- 类似的，内部无桥：每条边都至少在一个简单环中。
641  *         --- 存在点-双通分量 && 边-双通分量
642  *         --- 计算点-连通分量： 315页
643  * 6.例题：圆桌骑士
644  *         --- 以骑士为节点建立无向图G，如果两个骑士可以相邻，即不相互憎恨，则他们之间连接一条无向边，则题目转化为求不在任何一个简单奇圈上的结点个数。
645  * 7.有向图的强连通分量：
646  *         。。。
647  * 8.
648  *
649  */
650 
651 //DFS框架：
652 vector <int> G[10010];  //图使用vector式的邻接表，其中G[u][i]表示结点u的第i个子节点
653 int vis[10010]; //初始化为 0
654 void PreVisit(int x){}
655 void PostVist(int x){}  //可以在实际中完善其细节
656 void dfs(int u){ //u为起始结点
657     vis[u]=1;
658     PreVisit(u);
659     int d = G[u].size();
660     for(int i=0;i<d;i++){
661         int v = G[u][i];
662         if(!vis[v]) dfs(v);
663     }
664     PostVist(u);
665 }
666 
667 //连通分量：下面的代码为每一个结点计算出了该节点所属的连通分量编号
668 int current_cc;
669 int cc[10010]; //之后每个结点属于哪个连通分量的编号存储入数组 cc[] 中
670 int nn = 100; //nn表示图中点的个数
671 void find_cc(){
672     current_cc = 0;
673     memset(vis,0,sizeof(vis));
674     for(int u=0;u<nn;u++){
675         if(!vis[u]){
676             current_cc++;
677             dfs(u);  //在dfs的代码中，有Previsit()函数，在其中，将cc[u] = current_cc; 代表当前结点属于某个连通分量
678         }
679     }
680 }
681 
682 int color[10010];
683 //判断结点u所在的连通分量是否为二分图
684 bool bipartite(int u){
685     for(int i=0;i< (int)G[u].size();i++){
686         int v = G[u][i];
687         if(color[v] == color[u]) return false;
688         if(!color[v]){
689             color[v] = 3 - color[u];
690             if(!bipartite(v)) return false;
691         }
692     }
693     return true;
694 }
695 
696 
697 /* 5.3   最短路问题
698  *
699  * 1.单源最短路：可以计算任意点到其它所有点的最短路
700  * 2.定义的结构体： 保存加权图中的有向边，包括权值
701  * 3.例题：机场快线：
702  *         --- 快线分为商业线和经济线，只能做一站商业线，问最短的时间从起点到终点
703  *         --- 预处理单源最短时间，之后总时间为： f(a) + T(a,b) + f(b); T(a,b)枚举，取最小值
704  *         --- 思路特别巧妙： 使用单源最短路径进行预处理，之后按照d[B]<d[A]的条件，重新建图，题目的目标就是求出所有起点到终点的路径条数。（因为是DAG有向图，所以这里可以使用动态规划求解）
705  * 4.再谈Bellman-Ford算法：一个重要应用：判负圈
706  *         。。。
707  * 5.
708  *
709  */
710 
711 //定义有向图的结构体
712 struct Edge{
713     int from,to,dist;  //其中的dist会根据需求不同进行改变，如网络流中的容量和费用等。
714 };
715 struct HeapNode{
716     int d,u;
717     bool operator < (const HeapNode& rhs) const{
718         return d> rhs.d;
719     }
720 };
721 struct Dijkstra{
722     int n,m;   //点数和边数
723     vector<Edge> edges;  //边列表
724     vector<int> G[N];  //每个节点出发的边编号（从0开始编号）
725     bool done[N]; //是否已永久标号
726     int d[N];         //s到各个点的距离
727     int p[N];         //最短路中的上一条边
728 
729     void init(int n){
730         this->n = n;
731         for(int i=0;i<n;i++) G[i].clear(); //清空邻接表
732         edges.clear(); //清空边表
733     }
734 
735     void AddEdge(int from,int to,int dist){
736         //如果是无向图，每条无向边需要调用两次AddEdge()
737         edges.push_back((Edge){from,to,dist});
738         m = edges.size();
739         G[from].push_back(m-1);
740     }
741 
742     void dijkstra(int s){   //求s到所有点的距离 ---   调用之后，d[]数组之中就是存的距离，应该是表示到第i个结点的距离
743         priority_queue<HeapNode> Q;  //优先队列
744         for(int i=0;i<n;i++) d[i] = INF;
745         d[s] = 0;
746         memset(done,0,sizeof(done));
747         Q.push((HeapNode){0,s});
748         while(Q.empty()){
749             HeapNode x = Q.top();Q.pop();
750             int u = x.u;
751             if(done[u]) continue;
752             done[u] = true;
753             for(int i=0;i<(int)G[u].size();i++){
754                 Edge& e = edges[G[u][i]];
755                 if(d[e.to] > d[u] + e.dist){
756                     d[e.to] = d[u] + e.dist;
757                     p[e.to] = G[u][i];
758                     Q.push((HeapNode){d[e.to],e.to});
759                 }
760             }
761         }
762 
763     }
764 };
765 
766 
767 /* 5.4   生成树相关问题：最小生成树等  343页
768  *
769  * 1.最小生成树的两个性质：
770  *         --- 切割性质 && 回路性质
771  * 2.增量最小生成树：从包含n个点的空图开始，依次加入m条带权边。每加入一条边，输出当前图中的最小生成树权值。（如果当前图不连通，输出无解）O(n*m)
772  * 3.最小瓶颈生成树：给出加权无向图，求一棵生成树，使得最大边权值尽量小
773  * 4.最小瓶颈路给定加权无向图的两个节点u和v，求得从u到v的一条路径，使得路径上的最长边尽量短。
774  * 5.次小生成树：把所有生成树按照权值之和从大到小的顺序排列，求排在第二位的生成树。注意：如果最小生成树不唯一，次小生成树的权值和最小生成树相同。
775  * 6.最小有向生成树：指一个类似树的有向图
776  *
777  */
778 
779 /* 5.5   二分图
780  *
781  * 1.二分图最大匹配：找一个边数最大的匹配，使得任意两条选中的边均没有公共点。
782  * 2.二分图最佳完美匹配：假定有一个完全二分图G，每条边有一个权值（可以为负数）。如何求出权值和最大的完美匹配？（KM算法，也称匈牙利算法）
783  * 3.例题：蚂蚁匹配
784  *         --- 这里用到了一个特点，二分图的最佳匹配中不会出现线段相交的情况（最佳指所有的长度之和最小）
785  *         --- 换句话说：最佳匹配中不会出现线段相交的情况 ???
786  * 4.二分图有一些常见的应用：
787  *
788  */
789 
790 /* 5.6   网络流问题
791  *
792  * 1.
793  *
794  *
795  */
796 
797 
798 
799 
800 int main()
801 {
802     freopen("in.txt","r",stdin);
803     //freopen("out.txt","w",stdout);
804 
805 
806     return 0;
807 }