第7期：计算几何（持续更新中......）_如何利用delta函数和eplison函数证明三叉积-CSDN博客

本文链接：https://blog.csdn.net/qq_57351574/article/details/122644183

// 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标准库里的 $atan2$ 函数就是用来求极角的，如向量(x,y)的极角就是 $atan2(y,x)$ （单位：弧度）。

//极角函数 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的点积等于二者长度的乘积再乘上它们夹角 $\theta$ 的余弦。其中的 $\theta$ 是指从v到w逆时针旋转的角，因此当夹角大于90度时点积为负。

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

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组成的三角形的有向面积的两倍。不难发现，叉积不满足交换律。事实上， $cross(v,w)=-cross(w,v)$ 。在坐标系下，两个向量 $\overrightarrow{OA}$ 和 $\overrightarrow{OB}$ 的叉积等于 $x_{A}y_{B}-x_{B}y_{A}$ 。代码如下：

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的终点在左上方的第二象限时，点积为负但叉积均为正，用（-，+）表示。

向量旋转：向量可以绕起点旋转，公式为 $x^{'}=xcos\alpha -ysin\alpha,y^{'}=xsin\alpha+ycos\alpha$ ，其中 $\alpha$ 为逆时针旋转的角。代码如下：

//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)返回共轭复数，即 $conj(a+bi)=a-bi$ 。相关代码如下：

#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 点和直线

直线的参数表示

直线交点

点到直线的距离

点到线段的距离

点在直线上的投影

线段相交判定

直线的参数表示：

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

类型	t范围
直线	没有范围限制
射线	$t>0$
线段	$t\epsilon [0,1]$

直线交点：设直线分别为 $P+t\mathbf{v}$ 和 $Q+t\textbf{w}$ ，设向量 $\textbf{u}=\overrightarrow{QP}$ ，设交点在第一条直线上的参数为 $t_{1}$ ，第二条直线上的参数为 $t_{2}$ ，则x和y坐标可以各列出一个方程，解得：

$t_{1}=\frac{cross(w,u)}{cross(v,w)}$ ， $t_{2}=\frac{cross(v,u)}{cross(v,w)}$ 。代码如下：

//调用前请确保两条直线 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写成参数式 $A+t\textbf{v}$ （v为向量 $\overrightarrow{AB}$ ）,并且设Q的参数为 $t_{0}$ ，那么 $Q=A+t_{0}\mathbf{v}$ 。根据PQ垂直于AB，两个向量的点积应该为0，因此 $Dot(\textbf{v},P-(A+t_{0}\mathbf{v}))=0$ 。根据分配律有 $Dot(\mathbf{v},P-A)-t_{0}*Dot(\mathbf{v},\mathbf{v})=0$ ，这样就可以解出 $t_{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+F-E=2$ 。

这样，只需求出顶点数V和边数E，就可以求出 $F=E+2-V$ 。

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

/*********************************************************************************
欧拉定理：设平面图的顶点数、边数和面数分别为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

分析：

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

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

①第一种方法是解方程组。设交点为 $P=A+t(B-A)$ ，代入圆方程后整理得到 $(at+b)^{2}+(ct+d)^{2}=r^{2}$ ，进一步整理后得到一元二次方程 $et^{2}+ft+g=0$ 。根据判别式的值可以分成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，再求向量 $\overrightarrow{AB}$ 对应的单位向量 $\mathbf{v}$ ，则两个交点分别为 $P-L\mathbf{v}$ 和 $P+L\mathbf{v}$ ，其中L为P到交点的距离（P与两个交点等距），可以由勾股定理算出。代码略，代补。

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

// 计算向量极角
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;
}

过定点作圆的切线：先求出向量 $\overrightarrow{PQ}$ 的距离和向量 $\overrightarrow{PC}$ 的夹角ang，则向量 $\overrightarrow{PC}$ 的极角加减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度后即可知道切线的方向向量。这样，问题的关键是求出情况四、五中的外公切线和情况六中的内外公切线。

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

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

外公切线类似。假定 $r_{1}\geqslant r_{2}$ ，不管两圆是相离、相切还是相交， $cos\alpha$ 都是 $(r_{1}-r_{2})/d$ 。剩下的过程又和前面一样了。代码如下：

//返回切线的条数，-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 球面相关问题

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

$x=rcos\theta cos\phi$

$y=rcos\theta sin\phi,0\leqslant \theta \leqslant 2\pi,-\pi/2\leqslant \phi \leqslant \pi/2$

$z=rsin\theta$

代码如下：

//角度转换成弧度
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的圆弧长度。圆心角为 $2arcsin(d/2r)$ ，因此弧长为 $2arcsin(d/2r)r$ 。

2.3 例题选讲

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

2 UVA1308 Viva Confetti 圆盘问题

3 二维几何常用算法

点在多边形内的判定

凸包

半平面交

平面区域

4 三维几何基础

二维几何的内容也适用于三维几何，比如点+向量=点，向量+向量=向量，点+点没有定义等。代码如下：

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

typedef Point3 Vector3;

Vector3 operator + (Vector3 A, Vector3 B){
	return Vector3(A.x+B.x, A.y+B.y, A.z+B.z);
}

Vector3 operator - (Vector3 A, Vector3 B){
	return Vector3(A.x-B.x, A.y-B.y, A.z-B.z);
}

Vector3 operator * (Vector3 A, double p){
	return Vector3(A.x*p, A.y*p, A.z*p);
}

Vector3 operator / (Vector3 A, double p){
	return Vector3(A.x/p, A.y/p, A.z/p);
}

直线的表示：直线仍然可以用用参数方程（点和向量）来表示，射手和线段仍然可以看成“参数有取值范围限制”的直线，而且点到直线的投影也和二维情形一样。

平面的表示：通常用点法式 $(p_{0},\overrightarrow{n})$ 来描述一个平面。其中，点 $p_{0}$ 是平面上的一个点，向量 $\overrightarrow{n}$ 是平面的法向量。

每个平面把空间分为了两个部分，我们用点法式表示其中一个半空间。具体是哪一个呢？是这个法向量所背离的那一个（即法向量指向远离半空间的方向）。平面上的任意点p满足 $Dot(\overrightarrow{n},p-p_{0})=0$ 。设点p的坐标为 $(x_{0},y_{0},z_{0})$ ，向量 $\overrightarrow{n}$ 的坐标表示为 $(A,B,C)$ ，上述等式等价于 $A(x-x_{0})+B(y-y_{0})+C(z-z_{0})=0$ ，

整理得 $Ax+By+Cz-(Ax_{0}+By_{0}+Cz_{0})=0$ 。

如果令 $D=-(Ax_{0}+By_{0}+Cz_{0})$ ，我们就得到了

平面的一般式： $Ax+By+Cz+D=0$ 。注意，当 $Ax+By+Cz+D>0$ 时，上述点积大于0，即点 $(x,y,z)$ 在半空间 $(p_{0},\overrightarrow{n})$ 外。换句话说， $Ax+By+Cz+D>0$ 表示的是一个半空间（half space）。

4.1 三维点积

三维点积的定义和二维非常类似，而且也能用点积计算向量的长度和夹角。代码如下：

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

用三维点积可以解决很多基本问题。

过定点垂直于定直线的平面：平面的法向量就是这条直线，所以可以直接写出所求平面的点法式。

直线和平面的夹角、两平面的夹角、两直线的夹角：注意到，平面与平面的夹角可以转化为两者法向量之间的夹角，这些问题都不难解决。

点到平面的距离：把向量 $\overrightarrow{P_{0}P}$ 投影到法向量 $\overrightarrow{n}$ 上，可得点p到平面的有向距离为 $\frac{Dot(p-p_{0},\overrightarrow{n})}{Length(\overrightarrow{n})}$ 。这是一个相当简洁的结论，如果 $\overrightarrow{n}$ 是单位向量，甚至会更简单。

//点p到平面p0-n的距离。n必须为单位向量
double DistanceToPlane(const Point3& p, const Point3& p0, const Vector3& n){
	return fabs(Dot(p-p0,n)); //如果不取绝对值，得到的是有向距离 
}

点到平面上的投影：有了距离，投影点本身就不难求了。设点p到平面 $(p_{0},\overrightarrow{n})$ 上的投影为 $p'$ ，则 $p'-p$ 平行于 $\overrightarrow{n}$ ，则 $p'-p=d\overrightarrow{n}$ ，其中d就是p到平面的有向距离。代码如下：

//点p到平面p0-n上的投影。n必须为单位向量
Point3 GetPlaneProjection(const Point3& p, const Point3& p0, const Vector3& n){
	return p-n*Dot(p-p0,n);
}

直线和平面的交点：可以简单地通过解方程得到。设平面方程为 $Dot(\overrightarrow{n},p-p_{0})=0$ ，过点p1和p2的直线的参数方程为 $p=p_{1}+t(p_{2}-p_{1})$ ，则与平面方程联立解得： $t=\frac{Dot(\overrightarrow{n},p_{0}-p_{1})}{Dot(\overrightarrow{n},p_{2}-p_{1})}$ ，其中，分母为0的情况对应于直线和平面平行，或者直线在平面上。

//直线p1-p2到平面p0-n的交点，假定交点唯一存在
Point3 LinePlaneIntersection(Point3 p1, Point3 p2, Point3 p0, Vector3 n){
	Vector3 v = p2-p1;
	double t = (Dot(n,p0-p1) / Dot(n,p2-p1)); //判断分母是否为0
	return p1+v*t; //如果是线段，判断t是不是在0和1之间 
}

另，如果平面用一般式 $Ax+By+Cz+D=0$ ，则联立解出的表达式为： $t=\frac{Ax_{1}+By_{1}+Cz_{1}+D}{A(x_{1}-x_{2})+B(y_{1}-y_{2})+C(z_{1}-z_{2})}$

4.2 三维叉积

三维空间里也有叉积的概念，但形式和二维叉积不太一样。它是一个向量。

$v_{1}\ x\ v_{2}=\begin{bmatrix} y_{1}z_{2}-y_{2}z_{1}\\ z_{1}x_{2}-z_{2}x_{1}\\ x_{1}y_{2}-x_{2}y_{1} \end{bmatrix}$ ，虽然在学术上并不严谨，但在算法竞赛中可以认为叉积同时垂直于v1和v2，方向遵循右手定则。当且仅当v1和v2平行时，叉积为 $\overrightarrow{0}$ 。两空间向量叉积程序如下：

Vector3 Cross(Vector3 A, Vector3 B){
	return Vector3(A.y*B.z - A.z*B.y, A.z*B.x - A.x*B.z, A.x*B.y - A.y*B.x);
}
double Area2(Point3 A, Point3 B, Point3 C){ return Length(Cross(B-A, C-A)); }

有了叉积这个工具，可以解决更多的基础问题。

过不共线三点的平面：法向量为 $Cross(p_{2}-p_{0},p_{1}-p_{0})$ ，任取一个点即可得到平面的点法式。

三角形的有向面积：和二维情形一样，注意求出叉积之后要取长度。

判断点是否在三角形内：先判断点是否在三角形所在平面上，然后利用简单的面积关系即可判定，代码如下（假定点P在P0,P1,P2确定的平面上）：

//点P在三角形P0P1P2中
bool PointInTri(Point3 P, Point3 P0, Point3 P1,Point3 P2){
	double area1 = Area2(P, P0, P1);
	double area2 = Area2(P, P1, P2);
	double area3 = Area2(P, P2, P0);
	return dcmp(area1 + area2 + area3 - Area2(P0, P1, P2)) == 0;
}

判断线段和三角形是否相交：利用上面的函数不难判断，代码如下（为简单起见，这里没有考虑线段在平面上的情况）。

//三角形P0P1P2是否和线段AB相交
bool TriSegIntersection(Point3 P0, Point3 P1, Point3 P2, Point3 A, Point3 B, Point3& P){
	Vector3 n = Cross(P1-P0, P2-P0);
	if(dcmp(Dot(n, B-A)) == 0) return false; //线段AB和平面P0P1P2平行或共面
	else {                                  //平面A和直线P1-P2有唯一交点 
		double t = Dot(n, P0-A) / Dot(n, B-A);
		if(dcmp(t) < 0 || dcmp(t-1) > 0) return false; //交点不在线段AB上(在线段AB上是0~1)
		P = A + (B-A)*t; //计算交点
		return PointInTri(P, P0, P1, P2); //判断交点是否在三角形P0P1P2内 
	} 
}

点到直线/线段的距离：仍然可以用面积法（注意三维叉积是向量，要用Length函数而不是fabs）代码如下：

//点P到直线AB的距离
double DistanceToLine(Point3 P, Point3 A, Point3 B){
	Vector3 v1 = B - A, v2 = P - A;
	return Length(Cross(v1, v2)) / Length(v1);
} 

//点P到线段AB的距离
double DistanceToSegment(Point3 P, Point3 A, Point3 B){
	if(A == B) return Length(P-A);
	Vector3 v = 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 Length(Cross(v1, v2)) / Length(v1); 
}

四面体的体积：已知四面体的4个顶点A,B,C,D。根据叉积和点积的定义不难得出四面体的带符号体积为：

$V=\frac{1}{3}S*h=\frac{1}{6}(\overrightarrow{AB}\ x\ \overrightarrow{AC})*h=\frac{1}{6}((\overrightarrow{AB}\ x\ \overrightarrow{AC})*\overrightarrow{AD})$

其中， $\overrightarrow{AB},\overrightarrow{AC},\overrightarrow{AD}$ 呈右手系时为正。括号内部分也称为混合积。代码如下：

//返回AB,AC,AD的混合积，等于四面体ABCD的有向面积的6倍
double Volume6(Point3 A, Point3 B, Point3 C, Point3 D){
	return Dot(D-A, Cross(B-A, C-A));
}

多面体的体积：平面多边形的面积等于三角形的有向面积之和，空间多面体也类似。不过首先需要规定好多面体的存储方式。

一种简单的表示法是：点-面，即一个顶点数组V和面数组F。其中，V数组保存着各个顶点的空间坐标，而F数组保存着各个面的3个顶点在V数组中的索引。为了简单起见，假设各个面都是三角形，且3个点由右手定则确定的方向指向多边形的外部（即从外部看，3个顶点呈逆时针排列），所以这些面上3个点的排列顺序并不是任意的。