第7期：计算几何（持续更新中......）

2023-11-24 12:44:28

// 2022/01/22更新

更新内容主要为算法竞赛入门经典——训练指南（升级版）（刘汝佳、陈锋编著）第4章几何问题

1 二维几何基础

在平面坐标系下，向量和点一样，也用两个数x,y表示。第6章介绍齐次坐标的概念，从而在程序上区别开点和向量。以下点和向量都用两个数x,y表示。

不要把点和向量弄混。有如下点-点=向量，向量+向量=向量，向量-向量=向量，点+向量=点，而点+点是没有意义的。

struct Point{
	double x,y;
	Point(double x=0, double y=0):x(x),y(y) { } // 构造函数，方便代码编写 
};

typedef Point Vector; // 从程序实现上，Vector只是Point的别名

//向量+向量=向量，点+向量=点
Vector operator + (Vector A, Vector B) { return Vector(A.x+B.x, A.y+B.y); }
//点-点=向量
Vector operator - (Point A, Point B) { return Vector(A.x-B.x, A.y-B.y); }
//向量*数=向量
Vector operator * (Vector A, double p) { return Vector(A.x*p, A.y*p); }
//向量/数=向量
Vector operator / (Vector A, double p) { return Vector(A.x/p, A.y/p); }

bool operator < (const Point& a, const Point& b) {
	return a.x<b.x || ( a.x == b.x && a.y < b.y);
} 

const double eps = 1e-10;
int dcmp(double x){
	if(fabs(x) < eps) return 0; else return x < 0 ? -1 : 1;
}

bool operator == (const Point& a, const Point& b){
	return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y) == 0;
}

注意，上面的“相等”函数用到了“三态函数”dcmp，减少了精度问题。

极角：

向量的极角：从x轴正半轴旋转到该向量方向所需要的角度。
C标准库里的函数就是用来求极角的，如向量(x,y)的极角就是（单位：弧度）。

//极角函数 atan2()函数 的例子
#include<bits/stdc++.h>
using namespace std;
struct Point{
	double x,y;
	Point(double x=0, double y=0):x(x),y(y) { } // 构造函数，方便代码编写 
};
typedef Point Vector; // 从程序实现上，Vector只是Point的别名
int main(){
	ios::sync_with_stdio(false);
	double x,y,Polar_angle; Vector vec; 
	cin>>x>>y;
	vec=Vector(x,y);
	Polar_angle=atan2(y,x);
	cout<<Polar_angle<<endl;
	return 0;
}
/*
#1        #2        #3
0 1       -1 0      13 34
1.5708    3.14159   1.20559
*/

1.1 基本运算

点积

叉积

两个向量的位置关系

向量旋转

基于复数的几何计算

点积：两个向量v和w的点积等于二者长度的乘积再乘上它们夹角的余弦。其中的是指从v到w逆时针旋转的角，因此当夹角大于90度时点积为负。

余弦是偶函数，因此点积满足交换律。如果两向量垂直，点积等于0。不难推导出：在平面坐标系下，两个向量和的点积等于。这里给出点积计算方法，以及利用点积计算向量长度和夹角的函数。代码如下：

struct Point{
	double x,y;
	Point(double x=0, double y=0):x(x),y(y) { } // 构造函数，方便代码编写 
};

typedef Point Vector; // 从程序实现上，Vector只是Point的别名

double Dot(Vector A, Vector B) { return A.x*B.x + A.y*B.y; }
double Length(Vector A) { return sqrt(Dot(A, A)); }
double Angle(Vector A, Vector B) { return acos(Dot(A, B) / Length(A) / Length(B)); }

//2022/01/23

叉积：简单来说，两个向量v和w的叉积等于v和w组成的三角形的有向面积的两倍。不难发现，叉积不满足交换律。事实上，。在坐标系下，两个向量和的叉积等于。代码如下：

struct Point{
	double x,y;
	Point(double x=0, double y=0):x(x),y(y) { } // 构造函数，方便代码编写 
};

typedef Point Vector; // 从程序实现上，Vector只是Point的别名

//点-点=向量
Vector operator - (Point A, Point B) { return Vector(A.x-B.x, A.y-B.y); }

double Cross(Vector A, Vector B){ return  A.x*B.y-A.y*B.x; }
double Area2(Point A, Point B, Piont C){ return Cross(B-A, C-A); }

两个向量的位置关系：把叉积和点积组合在一起，我们可以更细致地判断两个向量的位置关系。

如图，括号里的第一个数是点积的符号，第二个数是叉积的符号，第一个向量v总是水平向右，另一个向量w的各种情况都包含在了上图中。比如，当w的终点在左上方的第二象限时，点积为负但叉积均为正，用（-，+）表示。

向量旋转：向量可以绕起点旋转，公式为，其中为逆时针旋转的角。代码如下：

//rad是弧度
Vector Rotate(Vector A, double rad){
	return Vector(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad)+A.y*cos(rad));
}

特殊情况：下面函数用来计算向量的单位法线，即左转90度以后把长度归一化。

//调用前请确保A不是零向量
Vector Normal(Vector A){
	double L=Length(A);
	return Vector(-A.y/L, A.x/L);
}

基于复数的几何计算：使用C++里的复数实现点和向量的方法。如下定义后，我们自动拥有了构造函数、加减法和数量积。用real(p)和imag(p)访问实部和虚部，conj(p)返回共轭复数，即。相关代码如下：

#include<complex>
using namespace std;
typedef complex<double> Point;
typedef Point Vector;

double Dot(Vector A, Vector B){ return real(conj(A)*B); }
double Cross(Vector A, Vector B){ return imag(conj(A)*B); }
Vector Rotate(Vector A, double rad){ return A*exp(Point(0, rad)); }

上述函数的效率不是很高，但是相当方便、好记。

1.2 点和直线

直线的参数表示

直线交点

点到直线的距离

点到线段的距离

点在直线上的投影

线段相交判定

直线的参数表示：

公式：(P0是直线上一点，v是方向向量，t为参数)，如果已知直线上的两个不同点A和B，则方向向量为B-A，所以参数方程为。参数方程最方便的地方在于直线、射线和线段的方程形式是一样的，区别仅仅在于参数t。

类型	t范围
直线	没有范围限制
射线
线段

直线交点：设直线分别为和，设向量，设交点在第一条直线上的参数为，第二条直线上的参数为，则x和y坐标可以各列出一个方程，解得：

，。代码如下：

//调用前请确保两条直线 P+tv 和 Q+tw 有唯一交点，当且仅当 Cross(v,w)非0
Point GetLineIntersection(Point P, Vector v, Point Q, Vector w){
	Vector u = P-Q;
	double t = Cross(w, u) / Cross(v, w);
	return P+v*t; 
}

需要注意的是，从上述公式可以看出，如果P,v,Q,w的各个分量均为有理数，则交点坐标也是有理数。在精度要求极高的情况下，可以考虑自定义分数类。（？？？）

点到线段的距离：点到线段的距离有两种可能。设投影点为Q。

①如果Q在线段AB上，则所求距离就是P点到直线AB的距离。

②如果Q不在线段AB上，则分两种情况：

1）如果Q在射线BA上，则所求距离为PA距离；

2）如果Q在射线AB上，则所求距离为PB距离。

判断Q的位置可以用点积进行。代码如下：

double DistanceToSegment(Point P, Point A, Point B){
	if(A == B) return Length(P-A);
	Vector v1 = B - A, v2 = P - A, v3 = P - B;
	if(dcmp(Dot(v1, v2)) < 0) return Length(v2);
	else if(dcmp(Dot(v1, v3)) > 0) return Length(v3);
	else return fabs(Cross(v1, v2)) / Length(v1);
}

点在直线上的投影：为求上述的投影点Q，我们把直线AB写成参数式（v为向量）,并且设Q的参数为，那么。根据PQ垂直于AB，两个向量的点积应该为0，因此。根据分配律有，这样就可以解出，从而得到Q点。代码如下：

Point GetLineProjection(Point P, Point A, Point B){
	Vector v = B-A;
	return A+v*(Dot(v, P-A) / Dot(v, v) ); 
}

线段相交判定：即给定两条线段，判断是否相交。

我们定义“规范相交”为两线段恰好有一个公共点，且公共点不是任何一条线段的端点。（以可以理解成两条线段均不含端点）

线段规范相交的充要条件是：每条线段的两个端点都在另一条线段的两侧（这里的“两侧”是指叉积的符号不同）。代码如下：

bool SegmentProperIntersection(Point a1, Point a2, Point b1, Point b2){
	double c1 = Cross(a2-a1,b1-a1), c2 = Cross(a2-a1,b2-b1);
		   c2 = Cross(b2-b1,a1-b1), c4 = Cross(b2-b1,a2-b1);
	return dcmp(c1)*dcmp(c2)<0 && dcmp(c3)*dcmp(c4)<0;
}

如果允许在端点处相交，情况就比较复杂了：

如果c1和c2都是0，表示两线段共线，这时可能会有部分重叠的情况；

如果c1和c2不都是0，则只有一种相交方法，即某个端点在另外一条线段上。

为了判断上述情况是否发生，还需要判断一个点是否在一条线段上（不包含端点）。代码如下：

bool OnSegment(Point p, Point a1, Point a2){
	return dcmp(Cross(a1-p, a2-p)) == 0 && dcmp(Dot(a1-p, a2-p) < 0);
}

1.3 多边形

计算多边形的有向面积

①凸多边形，可以从第一个顶点出发把凸多边形分成n-2个三角形，然后把面积加起来。代码如下

double ConvexPolygonArea(Point* p, int n){
	double area = 0;
	for(int i = 1; i < n-1; i++)
		area += Cross(p[i]-p[0], p[i+1]-p[0]);
	return area/2;
}

②非凸多边形，上述方法也对非凸多边形适用。由于三角形面积是有向的，在外面的部分可以正负抵消掉。实际上，可以从任意点出发进行划分。

可以取p[0]点为划分顶点，一方面可以少算两个叉积（0和任意向量的叉积都等于0），另一方面也减少了乘法溢出的可能性，还不用特殊处理（i=n-1的时候，下一个顶点是p[0]而不是p[n]，因为p[n]不存在）。代码如下：

//多边形的有向面积
double PolygonArea(Point* p, int n){
	double area = 0;
	for(int i = 1; i < n-1; i++)
		area += Cross(p[i]-p[0], p[i+1]-p[0]);
	return area/2;
}

也可以取坐标原点为划分点，乘法次数减少，代码待自行编写。

// 2022/01/25

1.4 例题选讲

1. UVA11178 Morley's Theorem(Morley定理)

#include<bits/stdc++.h>
using namespace std;
struct Point{
	double x,y;
	Point(double x=0, double y=0):x(x),y(y) { }
};
typedef Point Vector;
Vector operator + (Vector A, Vector B) { return Vector(A.x+B.x, A.y+B.y); }
Vector operator - (Point A, Point B) { return Vector(A.x-B.x, A.y-B.y); }
Vector operator * (Vector A, double p) { return Vector(A.x*p, A.y*p); }
Vector operator / (Vector A, double p) { return Vector(A.x/p, A.y/p); }
double Dot(Vector A, Vector B) { return A.x*B.x + A.y*B.y; }
double Length(Vector A) { return sqrt(Dot(A, A)); }
double Angle(Vector A, Vector B) { return acos(Dot(A, B) / Length(A) / Length(B)); }
double Cross(Vector A, Vector B){ return  A.x*B.y-A.y*B.x; }
Vector Rotate(Vector A, double rad){
	return Vector(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad)+A.y*cos(rad));
} 
Point GetLineIntersection(Point P, Vector v, Point Q, Vector w){
	Vector u = P-Q;
	double t = Cross(w, u) / Cross(v, w);
	return P+v*t; 
} 
Point read_point(){
	Point p;
	scanf("%lf %lf",&p.x,&p.y);
	return p;
}
Point getD(Point A, Point B, Point C){
	Vector v1 = C-B;
	double a1 = Angle(A-B, v1);
	v1 = Rotate(v1, a1/3);
	
	Vector v2 = B-C;
	double a2 = Angle(A-C, v2);
	v2 = Rotate(v2, -a2/3); // 负号表示顺时针旋转
	
	return GetLineIntersection(B, v1, C, v2); 
}
int main(){
	int T;
	Point A, B, C, D, E, F;
	scanf("%d",&T);
	while(T--){
		A = read_point();
		B = read_point();
		C = read_point();
		D = getD(A, B, C);
		E = getD(B, C, A);
		F = getD(C, A, B);
		printf("%.6lf %.6lf %.6lf %.6lf %.6lf %.6lf\n", D.x, D.y, E.x, E.y, F.x, F.y); 
	}
	return 0;
}

2 好看的一笔画 That Nice Euler Circuit，上海 2004，LA 3263

分析：(还需要思考)

若是要直接找出所有区域，会非常麻烦，而且容易出错。但用欧拉定理可以将问题进行转化，使解法更容易。

欧拉定理：设平面图的顶点数、边数和面数分别为V，E和F，则。

这样，只需求出顶点数V和边数E，就可以求出。

该平面图的结点由两部分组成，即原来的结点和新增的结点。由于可能出现三线共点，需要删除重复的点。代码如下：

/*********************************************************************************
欧拉定理：设平面图的顶点数、边数和面数分别为V、E和F，则V+F-E=2
so...F=E+2-V;
该平面图的结点有原来的和新增结点构成，由于可能出现三线共点，需要删除重复点
*********************************************************************************/
#include<bits/stdc++.h>
using namespace std;

struct Point{
    double x,y;
    Point(double x=0,double y=0):x(x),y(y){}
};

typedef Point Vector;

//向量+向量=向量;    向量+点=点
Vector operator + (Vector A,Vector B){return Vector(A.x+B.x,A.y+B.y);}

//点-点=向量
Vector operator - (Point A,Point B){return Vector(A.x-B.x,A.y-B.y);}

//向量*数=向量
Vector operator * (Vector A,double p){return Vector(A.x*p,A.y*p);}

//向量/数=向量
Vector operator / (Vector A,double p){return Vector(A.x/p,A.y/p);}


bool operator < (const Point& a,const Point& b){return a.x<b.x||(a.x==b.x && a.y<b.y);}


const double eps = 1e-10;

int dcmp(double x){if(fabs(x)<eps)return 0;else return x < 0 ? -1 : 1;}

bool operator == (const Point& a,const Point& b){return dcmp(a.x-b.x)==0 && dcmp(a.y-b.y)==0;}


double Dot(Vector A,Vector B){return A.x*B.x+A.y*B.y;}
double length(Vector A){return sqrt(Dot(A,A));}
double Angle(Vector A,Vector B){return acos(Dot(A,B)/length(A)/length(B));}

double Cross(Vector A,Vector B){return A.x*B.y-B.x*A.y;}
double Area2(Point A,Point B,Point C){return Cross(B-A,C-A);}

/*******两直线交点*******/
//调用前确保两条直线P+tv和Q+tv有唯一交点，当且仅当Cross(v,w)非0;
Point GetLineIntersection(Point P,Vector v,Point Q,Vector w){
    Vector u=P-Q;
    if(Cross(v,w))
    {
        double t=Cross(w,u)/Cross(v,w);//精度高的时候，考虑自定义分数类
        return P+v*t;
    }
//    else
 //       return ;
}

/************************
线段相交判定（规范相交）
************************/
bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2){
    double c1=Cross(a2-a1,b1-a1);
    double c2=Cross(a2-a1,b2-a1);
    double c3=Cross(b2-b1,a1-b1);
    double c4=Cross(b2-b1,a2-b1);
    return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;
}
/**如果允许在端点处相交：如果c1和c2都是0，表示共线，如果c1和c2不都是0，则表示某个端点在另一条直线上**/
bool OnSegment(Point p,Point a1,Point a2){
    return dcmp(Cross(a1-p,a2-p))==0&&dcmp(Dot(a1-p,a2-p))<0;
}

const int mmax=310;
Point P[mmax],V[mmax*mmax];

Point read_point(Point &P){
    scanf("%lf%lf",&P.x,&P.y);
    return P;
}

int main(){
    int n,kase=0;
    while(scanf("%d",&n)==1 && n){
        for(int i=0;i<n;i++){
            P[i]=read_point(P[i]);
            V[i]=P[i];
        }
        n--;
        int c=n,e=n;
        for(int i=0;i<n;i++){
            for(int j=i+1;j<n;j++){
                if(SegmentProperIntersection(P[i],P[i+1],P[j],P[j+1])){//严格相交
                    V[c++]=GetLineIntersection(P[i],P[i+1]-P[i],P[j],P[j+1]-P[j]);//交点
                }
            }
        }
 //       printf("c=%d\n",c);
        sort(V,V+c);
        c=unique(V,V+c)-V;
  //      printf("%d=%d-%d\n",c,unique(V,V+c),V);
        for(int i=0;i<c;i++){
            for(int j=0;j<n;j++){
                if(OnSegment(V[i],P[j],P[j+1])) e++;//边数
            }
        }
        printf("Case %d: There are %d pieces.\n",++kase,e+2-c);
    }
    return 0;
}

/*
5
0 0 0 1 1 1 1 0 0 0
7
1 1 1 5 2 1 2 5 5 1 3 5 1 1
7
1 1 1 5 2 1 2 5 5 1 3 9 1 1
0
*/

3 UVA11796 狗的距离 Dog Distance

分析：

先来看一种简单的情况：甲狗和乙狗的路线都是一条线段。因为运动是相对的，因此也可以认为甲狗静止不动，乙狗自己沿着直线走，因此问题转化为求点到线段的最小或最大距离。

有了简化版的分析，只需模拟整个过程。设现在甲狗的位置在，刚经过编号为的拐点；乙的位置是，刚经过编号为的拐点，则我们只需要计算两狗谁先到拐点，那么在这个时间点之前的问题就是我们刚才讨论过的“简化版”、求解完毕之后，需要更新甲狗和乙狗的位置，如果正好到达下一个拐点，还要更新和，然后继续模拟。因为每次至少有一条狗到达拐点，所以子问题的求解次数不超过A+B。代码如下：

#include<bits/stdc++.h>
using namespace std;
struct Point{
	double x,y;
	Point(double x=0, double y=0):x(x),y(y) { }
};
typedef Point Vector;
Vector operator + (Vector A, Vector B) { return Vector(A.x+B.x, A.y+B.y); }
Vector operator - (Point A, Point B) { return Vector(A.x-B.x, A.y-B.y); }
Vector operator * (Vector A, double p) { return Vector(A.x*p, A.y*p); }
Vector operator / (Vector A, double p) { return Vector(A.x/p, A.y/p); }
bool operator < (const Point& a, const Point& b) {
	return a.x<b.x || ( a.x == b.x && a.y < b.y);
} 
const double eps = 1e-10;
int dcmp(double x){
	if(fabs(x) < eps) return 0; else return x < 0 ? -1 : 1;
}
bool operator == (const Point& a, const Point& b){
	return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y) == 0;
}
double Dot(Vector A, Vector B) { return A.x*B.x + A.y*B.y; }
double Length(Vector A) { return sqrt(Dot(A, A)); }
double Angle(Vector A, Vector B) { return acos(Dot(A, B) / Length(A) / Length(B)); }
double Cross(Vector A, Vector B){ return  A.x*B.y-A.y*B.x; }
Vector Rotate(Vector A, double rad){
	return Vector(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad)+A.y*cos(rad));
} 
Point GetLineIntersection(Point P, Vector v, Point Q, Vector w){
	Vector u = P-Q;
	double t = Cross(w, u) / Cross(v, w);
	return P+v*t; 
} 
Point read_point(){
	Point p;
	scanf("%lf %lf",&p.x,&p.y);
	return p;
}
double DistanceToSegment(Point P, Point A, Point B){
	if(A == B) return Length(P-A);
	Vector v1 = B - A, v2 = P - A, v3 = P - B;
	if(dcmp(Dot(v1, v2) < 0)) return Length(v2);
	else if(dcmp(Dot(v1, v3)) > 0) return Length(v3);
	else return fabs(Cross(v1, v2)) / Length(v1); 
}
const int maxn = 60;
int T, A, B;
Point P[maxn], Q[maxn];
double Min, Max;
void update(Point P, Point A, Point B){
	Min = min(Min, DistanceToSegment(P, A, B));
	Max = max(Max, Length(P-A));
	Max = max(Max, Length(P-B));
}
int main(){
	scanf("%d",&T);
	for(int kase = 1; kase <= T; kase++){
		scanf("%d%d",&A,&B);
		for(int i=0;i<A;i++) P[i]=read_point();
		for(int i=0;i<B;i++) Q[i]=read_point();
		
		double LenA=0,LenB=0;
		for(int i=0;i<A-1;i++) LenA+=Length(P[i+1]-P[i]);
		for(int i=0;i<B-1;i++) LenB+=Length(Q[i+1]-Q[i]);
		
		int Sa=0,Sb=0;
		Point Pa=P[0],Pb=Q[0];
		Min=1e9,Max=-1e9;
		while(Sa<A-1&&Sb<B-1){
			double La=Length(P[Sa+1]-Pa); //甲到下一拐点的距离 
			double Lb=Length(Q[Sb+1]-Pb); //乙到下一拐点的距离 
			double T=min(La/LenA,Lb/LenB);
			//取合适的单位，可以让甲和乙的速度分别是LenA和LenB 
			Vector Va=(P[Sa+1]-Pa)/La*T*LenA; //甲的位移向量 
			Vector Vb=(Q[Sb+1]-Pb)/Lb*T*LenB; //乙的位移向量 
			update(Pa, Pb, Pb+Vb-Va); //求解“简化版”，更新最小最大距离 
			Pa=Pa+Va;
			Pb=Pb+Vb;
			if(Pa==P[Sa+1]) Sa++;
			if(Pb==Q[Sb+1]) Sb++;
		}
		printf("Case %d: %.0lf\n",kase,Max-Min); 
	}
	return 0;
}

2 与圆和球有关的计算问题

2.1 圆的相关计算

圆上任意一点拥有唯一的圆心角，所以在定义圆的时候，可以加一个通过圆心角求坐标的函数。代码如下：

struct Point{
	double x,y;
	Point(double x, double y):x(x), y(y) { }
};
struct Circle {
	Point c; //圆心
	double r; //半径
	Circle(Point c, double r):c(c), r(r) { }
	Point point(double a){
		return Point(c.x + cos(a)*r, c.y + sin(a)*r);
	} 
};

直线和圆的交点：假设直线为AB，圆的圆心为C，半径为r。

①第一种方法是解方程组。设交点为，代入圆方程后整理得到，进一步整理后得到一元二次方程。根据判别式的值可以分成3种情况，即无交点（相离）、一个交点（相切）和两个交点（相交）。代码如下（变量a,b,c,d,e,f,g对应于上述方程中的字母）：

int getLineCircleIntersection(Line L, Circle C, double& t1, double& t2, Vector<Point>& sol){
	double a = L.v.x, b = L.p.x-C.c.x, c = L.v.y, d = L.p.y-C.c.y;
	double e = a*a+c*c, f = 2*(a*b+c*d), g = b*b+d*d-C.r*C.r;
	double delta = f*f-4*e*g; //判别式
	if(dcmp(delta) < 0) return 0; //相离
	if(dcmp(delta) == 0){ //相切 
		t1=t2= -f / (2*e); sol.push_back(L.point(t1));
		return 1;
	}
	//相交
	t1 = (-f - sqrt(delta)) / (2*e); sol.push_back(L.point(t1)); 
	t2 = (-f + sqrt(delta)) / (2*e); sol.push_back(L.point(t2));
	return 2;
}

函数返回的是交点的个数，参数sol存放的是交点本身。注意，上述代码并没有清空sol，这给很多题目带来方便：可以反复调用这个函数，把所有交点放在一个sol里。

②第二种方法是几何法。先求圆心C在AB上的投影P，再求向量对应的单位向量，则两个交点分别为和，其中L为P到交点的距离（P与两个交点等距），可以由勾股定理算出。代码略，代补。

两圆相交：假设圆心分别为和，半径分别为和，圆心距为，根据余弦定理可以算出到的角，由向量的极角，加减就可以得到和的极角。有了极角，就可以很方便地计算出和的坐标了（和是两圆交点），代码如下：

// 计算向量极角
double angle(Vector v){ return atan2(v.y, v.x); }
// 两圆相交 
int getCircleCircleIntersection(Circle C1, Circle C2, vector<Point>& sol){
	double d = Length(C1.c - C2.c); //圆心距
	if(dcmp(d) == 0){
		if(dcmp(C1.r - C2.r) == 0) return -1; //两圆重合
		return 0; //两个同心圆 
	} 
	if(dcmp(C1.r + C2.r - d) < 0) return 0; //两圆相离
	if(dcmp(fabs(C1.r-C2.r)-d) > 0) return 0; 
	
	double a = angle(C2.c - C1.c); //向量C1C2的极角
	double da = acos((C1.r*C1.r + d*d - C2.r*C2.r) / (2*C1.r*d));
	//C1C2到C1P1的角
	Point p1 = C1.point(a-da), p2 = C1.point(a+da);
	
	sol.push_back(p1); 
	if(p1 == p2) return 1;
	sol.push_back(p2);
	return 2;
}

过定点作圆的切线：先求出向量的距离和向量的夹角ang，则向量的极角加减ang就是两条切线的极角，注意切线不存在和只有一条的情况。代码如下：

//过点p到圆C的切线，v[i]是第i条切线的向量，返回切线条数
int getTangents(Point p, Circle C, Vector* v){
	Vector u = C.c - p;
	double dist = Length(u);
	if(dist < C.r) return 0;
	else if(dcmp(dist-C.r) == 0){ //p在圆上，只有一条切线
		v[0] = Rotate(u, PI/2);
		return 1; 
	}else{
		double ang = asin(C.r / dist);
		v[0] = Rotate(u, -ang);
		v[1] = Rotate(u, +ang);
		return 2;
	}
}

两圆的公切线：根据两圆的圆心距，从小到大排列一共有6种情况：

①情况一：两圆完全重合，有无数条公切线。

②情况二：两圆内含，没有公共点，没有公切线。

③情况三：两圆内切，有1条外公切线。

④情况四：两圆相交，有2条外公切线。

⑤情况五：两圆外切，有3条公切线（1条内公切线，2条外公切线）。

⑥情况六：两圆外离，有4条公切线（2条内公切线，2条外公切线）。

可以根据圆心距和半径的关系辨别出这6种情况，然后逐一求解。情况一和情况二没什么需要求解的；情况三和情况五中的内公切线都对应于“过圆上一点求圆的切线”，只需连接圆心和切点，旋转90度后即可知道切线的方向向量。这样，问题的关键是求出情况四、五中的外公切线和情况六中的内外公切线。

先考虑情况六中的内公切线，它对应于两圆相离的情况。

根据三角函数定义不难求出角度，然后和前面一样通过极角进行计算即可。

外公切线类似。假定，不管两圆是相离、相切还是相交，都是。剩下的过程又和前面一样了。代码如下：

//返回切线的条数，-1表示有无穷多条切线
//a[i]和b[i]分别是第i条切线在圆A和圆B上的切点
int getTangents(Circle A, Circle B, Point* a, Point* b){
	int cnt = 0;
	if(A.r < B.r) { swap(A, B); swap(a, b); }
	int d2 = (A.x-B.x)*(A.x-B.x) + (A.y-B.y)*(A.y-B.y);
	int rdiff = A.r-B.r;
	int rsum = A.r+B.r;
	if(d2 < rdiff*rdiff) return 0; //内含
	
	double base = atan2(B.y-A.y, B.x-A.x);
	if(d2 == 0 && A.r == B.r) return -1; //无限多条切线
	if(d2 == rdiff*rdiff){ //内切，1条切线 
		a[cnt] = A.getPoint(base); b[cnt] = B.getPoint(base); cnt++;
		return 1; 
	} 
	//有外公切线
	double ang = acos((A.r-B.r) / sqrt(d2));
	a[cnt] = A.getPoint(base+ang); b[cnt] = B.getPoint(base+ang); cnt++;
	a[cnt] = A.getPoint(base-ang); b[cnt] = B.getPoint(base-ang); cnt++;
	if(d2 == rsum*rsum){ //一条内公切线
		a[cnt] = A.getPoint(base); b[cnt] = B.getPoint(PI+base); cnt++; 
	}else if(d2 > rsum*rsum){ //两条内公切线
		double ang = acos((A.r+B.r) / sqrt(d2));
		a[cnt] = A.getPoint(base+ang); b[cnt] = B.getPoint(PI+base+ang); cnt++;
		a[cnt] = A.getPoint(base-ang); b[cnt] = B.getPoint(PI+base-ang); cnt++; 
	}
	return cnt;
}

2.2 球面相关问题

经纬度转换为空间坐标：公式如下：

代码如下：

//角度转换成弧度
double torad(double deg){
	return deg/180 *acos(-1); //acos(-1)就是PI 
} 

//经纬度（角度）转化为空间坐标
void get_coord(double R, double lat, double lng, double& x, double& y, double& z){
	lat = torad(lat);
	lng = torad(lng);
	x = R*cos(lat)*cos(lng);
	y = R*cos(lat)*sin(lng);
	z = R*sin(lat);
}

球面距离：问题：已知球面两点，如何求出它们的最短路?注意，只能沿着球面走，不能穿过球的内部。从表面走的话，最近的路径是走圆弧，具体来说是走大圆(Great Circle)圆弧。用一个穿过球心的平面截这个球，截面就是一个大圆。

怎么求大圆弧长呢？你无须想象那个大圆的空间位置，而可以把它们想象成一个平面问题：求半径为r，弦长为d的圆弧长度。圆心角为，因此弧长为。

2.3 例题选讲

1 UVA12304 2D Geometry 110 in 1! 二维几何110合一！

2 UVA1308 Viva Confetti 圆盘问题

3 二维几何常用算法

点在多边形内的判定

凸包

半平面交

平面区域

码农公寓