ACM几何基础篇

1 前言

1.1 生产生活

随着科学技术的飞速发展及计算机在国民经济各个领域中的普遍运用，计算机辅助设计，即CAD越来越为人们所重视。当前的CAD工作中，计算机远远不只是一种高效的计算工具，它已成为人们进行创造性设计活动的得力助手甚至参谋。计算几何作为CAD的基础理论之一，主要研究内容是几何形体的数学描述和计算机表述；它同计算机辅助几何设计，即CAGD有着十分密切的关系。而CAGD是由微分几何、代数几何、数值计算、逼近论、拓扑学以及数控技术等形成的一门新兴边缘学科，其主要研究对象和内容是对自由形曲线、曲面的数学描述、设计、分析及图形的显示、处理等。

1.2 ICPC竞赛

• 在竞赛中，计算几何相比于其他部分来说是比较独立的

• 曾出现过其与图论和动态规划相结合的题目，计算几何与其他的知识点很少有过多的结合

• 计算几何的题目难度不会很大，但也永远不会成为最水的题

• 计算几何题目具有代码量大、特殊情况多。精度问题难以控制等特点

2 准备知识

2.1 头文件及函数及常量

#include <cstdio>
#include <cmath>

using namespace std;

const double pi = acos(-1.0);
const double inf = 1e100;
const double eps = 1e-6;
int sgn(double d){
    if(fabs(d) < eps)
        return 0;
    if(d > 0)
        return 1;
    return -1;
}

int main() {
    double x = 1.49999;
    int fx = floor(x);//向下取整函数
    int cx = ceil(x);//向上取整函数
    int rx = round(x);//四舍五入函数
    printf("%f %d %d %d\n", x, fx, cx, rx);
    //输出结果 1.499990 1 2 1
    return  0 ;
}

2.2 浮点误差

计算几何中一般来说使用$double$类型比较频繁，该用实数是就用$double$，$float$容易失去精度

$double$类型读入时应用$\%lf$占位符

而输出时应用$\%f$占位符

2.2.1 计算误差

尽量少用三角函数、除法、开方、求幂、取对数运算

• 尽量把公式整理成使用的以上操作最少的形式

• $1.0/2.03.0/4.05.0 = (1.03.05.0)/(2.0*4.0)$

• 在不溢出的情况下将除式比较转化为乘式比较

• $\frac {a}{b} > c \Leftrightarrow a > bc$

在使用除法、开根号、和三角函数的时候，我们要考虑由浮点误差运算产生的代价

2.2.2 判等

不要直接用等号判断浮点数是否相等！
不要直接用等号判断浮点数是否相等！
不要直接用等号判断浮点数是否相等！

同样的一个数$1.5$如果从不同的途径计算得来，它们的储存方式可能会是$a = 1.4999999$ 与 $b = 1.5000001$，此时使用$a == b$判断将返回$false$

2.2.2.1 解决方案 1 误差判别法

借助$\varepsilon$也就是上面定义的常量eps

const double eps = 1e-9;
int dcmp(double x, double y){
    if(fabs(x - y) < eps)
        return 0;
    if(x > y)
        return 1;
    return -1;
}

2.2.2.2 解决方案 2 化浮为整

在不溢出整数范围的情况下，可以通过乘上$10 ^ k$转化为整数运算，最后再将结果转化为浮点数输出

2.2.3 负零

输出时一定要小心不要输出 $-0$，比如

double a = -0.000001;
printf("%.4f", a);

2.2.4 反三角函数

使用反三角函数时，要注意定义域的范围，比如，经过计算 x = 1.000001

double x = 1.000001;
double acx = acos(x);
//可能会返回runtime error

//此时我们可以加一句判断
double x = 1.000001;
if(fabs(x - 1.0) < eps || fabs(x + 1.0) < eps)
    x = round(x);
double acx = acos(x);

2.3 模板总结

#include <cstdio>
#include <cmath>

using namespace std;

const double pi = acos(-1.0);
const double inf = 1e100;
const double eps = 1e-6;
int sgn(double d){
    if(fabs(d) < eps)
        return 0;
    if(d > 0)
        return 1;
    return -1;
}
int dcmp(double x, double y){
    if(fabs(x - y) < eps)
        return 0;
    if(x > y)
        return 1;
    return -1;
}
int main() {
    double x = 1.49999;
    int fx = floor(x);//向下取整函数
    int cx = ceil(x);//向上取整函数
    int rx = round(x);//四舍五入函数
    printf("%f %d %d %d\n", x, fx, cx, rx);
    //输出结果 1.499990 1 2 1
    return  0 ;
}

3 点与向量

3.2 相关定义

3.1.1 点的定义

二维平面下的点的表示只需要两个实数，即横坐标与纵坐标即可

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

3.1.2 向量的定义

• 既有大小又有方向的量叫做向量
• 在计算机中我们常用坐标表示

这样看来，向量这个结构体貌似与点没有任何区别，因此我们可以

typedef Point Vector;

3.2 运算

3.2.1 加减乘除

3.2.1.1 加法运算

• 点与点之间的加法运算没有意义

• 点与向量相加得到另一个点

• 向量与向量相加得到另外一个向量

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

3.2.1.2 减法运算

两个点之间作差将得到一个向量，$A - B$将得到由$B$指向$A$的向量$\vec{BA}$

Vector operator - (Point A, Point B){
    return Vector(A.x-B.x, A.y-B.y);
}

3.2.1.3 乘法运算

向量与实数相乘得到等比例缩放的向量

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

3.2.1.4 除法运算

向量与实数相除得到等比例缩放的向量

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

3.2.2 小于运算（Left Then Low 排序）

有时我们需要将点集按照$x$坐标升序排列，若$x$坐标相同，则按照$y$坐标升序排列

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

此比较器将在Andrew算法中用到

而Graham Scan算法用到的比较器基于极角排序

3.2.3 等于运算

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

3.2.4 内积运算

又称数量积，点积

$\alpha \cdot \beta = |\alpha||\beta|cos\theta$

对加法满足分配律

3.2.4.1 几何意义

向量$\alpha$在向量$\beta$的投影$\alpha ’$（带有方向性）与$\beta$的长度乘积

• 若$\alpha$与$\beta$的夹角为锐角，则其内积为正

• 若$\alpha$与$\beta$的夹角为钝角，则其内积为负

• 若$\alpha$与$\beta$的夹角为直角，则其内积为0

3.2.4.2 代码实现

常用的实现方法有重载*运算符，或是单独写成函数，下面给出后一种实现方式

double Dot(Vector A, Vector B){
    return A.x*B.x + A.y*B.y;
}

3.2.5 外积运算

又称向量积，叉积

$\alpha \times \beta = |\alpha||\beta|sin\theta$

$\theta$表示向量$\alpha$旋转到向量$\beta$所经过的夹角

对加法满足分配律

3.2.5.1 几何意义

向量$\alpha$与$\beta$所张成的平行四边形的有向面积

3.2.5.2 判断外积的符号

右手定则

$\alpha \times \beta$

若$\beta$在$\alpha$的逆时针方向，则为正值

顺时针则为负值

两向量共线则为0

3.2.5.3 代码实现

重载^运算符或者单独写成函数

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

3.3 常用函数

3.3.1 取模（长度）

double Length(Vector A){
    return sqrt(Dot(A, A));
}

3.3.2 计算两向量夹角

返回值为弧度制下的夹角

double Angle(Vector A, Vector B){
    return acos(Dot(A, B) / Length(A) / Length(B));
}

3.3.3 计算两向量构成的平行四边形有向面积

double Area2(Point A, Point B, Point C){
    return Cross(B-A, C-A);
}

3.3.4 计算向量逆时针旋转后的向量

Vector Rotate(Vector A, double rad){//rad为弧度 且为逆时针旋转的角
    return Vector(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad)+A.y*cos(rad));
}

3.3.5 计算向量逆时针旋转九十度的单位法向量

Vector Normal(Vector A){//向量A左转90°的单位法向量
    double L = Length(A);
    return Vector(-A.y/L, A.x/L);
}

3.3.6 ToLeftTest

判断折线$\vec{bc}$是不是向$\vec{ab}$的逆时针方向（左边）转向

凸包构造时将会频繁用到此公式

bool ToLeftTest(Point a, Point b, Point c){
    return Cross(b - a, c - b) > 0;
}

3.4 复数黑科技

利用复数黑科技实现平面点与向量

复数定义向量后，自动拥有构造函数、加减法和数量积

3.4.1 代码实现

#include <complex>
using namespace std;
typedef complex<double> Point;
typedef Point Vector;
const double eps = 1e-9;
int dcmp(double x){
    if(fabs(x) < eps)
        return 0;
    if(x < 0)
        return -1;
    return 1;
}
double Length(Vector A){
    return abs(A);
}
double Dot(Vector A, Vector B){//conj(a+bi)返回共轭复数a-bi
    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));//exp(p)返回以e为底复数的指数
}

3.4.2 不足

复数运算会比自己写的向量运算慢，若题目时间要求比较苛刻，要谨慎使用

3.5 模板总结

3.5.1 手动实现

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){
    if(a.x == b.x)
        return a.y < b.y;
    return a.x < b.x;
}
const double eps = 1e-6;
int sgn(double x){
    if(fabs(x) < eps)
        return 0;
    if(x < 0)
        return -1;
    return 1;
}
bool operator == (const Point& a, const Point& b){
    if(sgn(a.x-b.x) == 0 && sgn(a.y-b.y) == 0)
        return true;
    return false;
}
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;
}
double Area2(Point A, Point B, Point C){
    return Cross(B-A, C-A);
}
Vector Rotate(Vector A, double rad){//rad为弧度 且为逆时针旋转的角
    return Vector(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad)+A.y*cos(rad));
}
Vector Normal(Vector A){//向量A左转90°的单位法向量
    double L = Length(A);
    return Vector(-A.y/L, A.x/L);
}
bool ToLeftTest(Point a, Point b, Point c){
    return Cross(b - a, c - b) > 0;
}

3.5.2 复数黑科技

#include <complex>
using namespace std;
typedef complex<double> Point;
typedef Point Vector;//复数定义向量后，自动拥有构造函数、加减法和数量积
const double eps = 1e-9;
int sgn(double x){
    if(fabs(x) < eps)
        return 0;
    if(x < 0)
        return -1;
    return 1;
}
double Length(Vector A){
    return abs(A);
}
double Dot(Vector A, Vector B){//conj(a+bi)返回共轭复数a-bi
    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));//exp(p)返回以e为底复数的指数
}

4 点与线

4.1 定义

4.1.1 直线定义

直线表示常用的有三种形式

• 一般式$ax + by + c = 0$

• 点向式$x_0 + y_0 + v_xt + v_yt = 0$

• 斜截式$y = kx + b$

计算机中常用点向式表示直线，即参数方程形式表示

4.1.2 点向式

直线可以用直线上的一个点$P_0$和方向向量$v$表示

$$P = P _ 0 + vt$$

其中$t$为参数

4.1.2.1 优点

• 可以表示所有直线

• 可以通过限制参数来表示线段和射线

4.1.3 线段与射线

利用带参数限制的直线点向式方程表示

4.2 实现

struct Line{//直线定义
    Point v, p;
    Line(Point v, Point p):v(v), p(p) {}
    Point point(double t){//返回点P = v + (p - v)*t
        return v + (p - v)*t;
    }
};

4.3 常用操作

4.3.1 判断点在直线上

• 利用三点共线的等价条件$\alpha \times \beta == 0$

• 直线上取两不同点与待测点构成向量求叉积是否为零

4.3.2 计算两直线交点

必须保证直线相交，否则将会出现除以零的情况

//调用前需保证 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;
}

4.3.3 计算点到直线的距离

//点P到直线AB距离公式
double DistanceToLine(Point P, Point A, Point B){
    Vector v1 = B-A, v2 = P-A;
    return fabs(Cross(v1, v2)/Length(v1));
}//不去绝对值，得到的是有向距离

4.3.4 计算点到线段的距离

//点P到线段AB距离公式
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);
    if(dcmp(Dot(v1, v3)) > 0)
        return Length(v3);
    return DistanceToLine(P, A, B);
}

4.3.5 求点在直线上的投影点

//点P在直线AB上的投影点
Point GetLineProjection(Point P, Point A, Point B){
    Vector v = B-A;
    return A+v*(Dot(v, P-A)/Dot(v, v));
}

4.3.6 判断点是否在线段上

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

4.3.7 判断两线段是否相交

通过两次跨立实验

不允许在顶点处相交

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

允许在端点处相交

bool SegmentProperIntersection(Point a1, Point a2, Point b1, Point b2){
    double c1 = Cross(a2-a1, b1-a1), c2 = Cross(a2-a1, b2-a1);
    double c3 = Cross(b2-b1, a1-b1), c4 = Cross(b2-b1, a2-b1);
    //if判断控制是否允许线段在端点处相交，根据需要添加
    if(!sgn(c1) || !sgn(c2) || !sgn(c3) || !sgn(c4)){
        bool f1 = OnSegment(b1, a1, a2);
        bool f2 = OnSegment(b2, a1, a2);
        bool f3 = OnSegment(a1, b1, b2);
        bool f4 = OnSegment(a2, b1, b2);
        bool f = (f1|f2|f3|f4);
        return f;
    }
    return (sgn(c1)*sgn(c2) < 0 && sgn(c3)*sgn(c4) < 0);
}

4.4 模板总结

struct Line{//直线定义
    Point v, p;
    Line(Point v, Point p):v(v), p(p) {}
    Point point(double t){//返回点P = v + (p - v)*t
        return v + (p - v)*t;
    }
};
//计算两直线交点
//调用前需保证 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到直线AB距离公式
double DistanceToLine(Point P, Point A, Point B){
    Vector v1 = B-A, v2 = P-A;
    return fabs(Cross(v1, v2)/Length(v1));
}//不去绝对值，得到的是有向距离
//点P到线段AB距离公式
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);
    if(dcmp(Dot(v1, v3)) > 0)
        return Length(v3);
    return DistanceToLine(P, A, B);
}
//点P在直线AB上的投影点
Point GetLineProjection(Point P, Point A, Point B){
    Vector v = B-A;
    return A+v*(Dot(v, P-A)/Dot(v, v));
}
//判断p点是否在线段a1a2上
bool OnSegment(Point p, Point a1, Point a2){
    return dcmp(Cross(a1-p, a2-p)) == 0 && dcmp(Dot(a1-p, a2-p)) < 0;
}
//判断两线段是否相交
bool SegmentProperIntersection(Point a1, Point a2, Point b1, Point b2){
    double c1 = Cross(a2-a1, b1-a1), c2 = Cross(a2-a1, b2-a1);
    double c3 = Cross(b2-b1, a1-b1), c4 = Cross(b2-b1, a2-b1);
    //if判断控制是否允许线段在端点处相交，根据需要添加
    if(!sgn(c1) || !sgn(c2) || !sgn(c3) || !sgn(c4)){
        bool f1 = OnSegment(b1, a1, a2);
        bool f2 = OnSegment(b2, a1, a2);
        bool f3 = OnSegment(a1, b1, b2);
        bool f4 = OnSegment(a2, b1, b2);
        bool f = (f1|f2|f3|f4);
        return f;
    }
    return (sgn(c1)*sgn(c2) < 0 && sgn(c3)*sgn(c4) < 0);
}

5 多边形

5.1 三角形

5.1.1 三角形面积

• 利用两条边叉积除以二取绝对值

• 海伦公式

$$S = \sqrt{p(p - a)(p - b)(p - c)}, p = \frac{(a + b + c)}{2}$$

• $S = \frac{absinC}{2}$

5.1.2 三角形四心

5.1.2.1 外心

三边中垂线交点，到三角形三个顶点距离相同

5.1.2.2 内心

角平分线的交点，到三角形三边的距离相同

5.1.2.3 垂心

三条高线的交点

5.1.2.4 重心

三条中线的交点，到三角形三顶点距离的平方和最小的点，三角形内到三边距离之积最大的点

5.2 普通多边形

通常按照逆时针储存所有顶点

5.2.1 定义

5.2.1.1 多边形

由在同一平面且不再同一直线上的多条线段首位顺次连接且不相交所组成的图形交多边形

5.2.1.2 简单多边形

简单多边形是除相邻边外其它边不相交的多边形

5.2.1.3 凸多边形

过多边形的任意一边做一条直线，如果其他各个顶点都在这条直线的同侧，则把这个多边形叫做凸多边形

任意凸多边形外角和均为$360°$

任意凸多边形内角和为$(n - 2)180°$

5.2.2 常用函数

5.2.2.1 求多边形面积

我们可以从第一个顶点除法把凸多边形分成$n - 2$个三角形，然后把面积加起来

最后返回值说为有向面积更贴近本质

double PolygonArea(Point* p, int n){//p为端点集合，n为端点个数
    double s = 0;
    for(int i = 1; i < n-1; ++i)
        s += Cross(p[i]-p[0], p[i+1]-p[0]);
    return s;
}

5.2.2.2 判断点在多边形内

有射线法与转角法。

转角法的基本思想是看多边形相对于这个点转了多少度

• 如果是三百六十度，说明点在多边形内

• 如果是零度，说明点在多边形外

• 如果是一百八十度，说明点在多边形边界上

如果直接按照定义来算，则需要计算大量反三角函数，不仅速度慢，而且容易产生精度问

我们采用winding number绕数来计算

//判断点是否在多边形内，若点在多边形内返回1，在多边形外部返回0，在多边形上返回-1
int isPointInPolygon(Point p, vector<Point> poly){
    int wn = 0;
    int n = poly.size();
    for(int i = 0; i < n; ++i){
        if(OnSegment(p, poly[i], poly[(i+1)%n])) return -1;
        int k = sgn(Cross(poly[(i+1)%n] - poly[i], p - poly[i]));
        int d1 = sgn(poly[i].y - p.y);
        int d2 = sgn(poly[(i+1)%n].y - p.y);
        if(k > 0 && d1 <= 0 && d2 > 0) wn++;
        if(k < 0 && d2 <= 0 && d1 > 0) wn--;
    }
    if(wn != 0)
        return 1;
    return 0;
}

5.2.2.4 判断点在凸多边形内

只需要判断点是否在所有边的左边（按逆时针顺序排列的顶点集）ToLeftTest

5.3 Pick定理

5.3.1 内容

皮克定理是指一个计算点阵中顶点在格点上的多边形面积公式该公式可以表示为

$$2S = 2a + b - 2$$

其中$a$表示多边形内部的点数，$b$表示多边形边界上的点数，$S$表示多边形的面积。

常用形式

$$S = a + \frac{b}{2} - 1$$

5.3.2 常用计算

给你多边形的顶点，问多边形内部有多少点

$a = S - \frac{b}{2} + 1$

5.4 模板总结

//多边形有向面积
double PolygonArea(Point* p, int n){//p为端点集合，n为端点个数
    double s = 0;
    for(int i = 1; i < n-1; ++i)
        s += Cross(p[i]-p[0], p[i+1]-p[0]);
    return s;
}
//判断点是否在多边形内，若点在多边形内返回1，在多边形外部返回0，在多边形上返回-1
int isPointInPolygon(Point p, vector<Point> poly){
    int wn = 0;
    int n = poly.size();
    for(int i = 0; i < n; ++i){
        if(OnSegment(p, poly[i], poly[(i+1)%n])) return -1;
        int k = sgn(Cross(poly[(i+1)%n] - poly[i], p - poly[i]));
        int d1 = sgn(poly[i].y - p.y);
        int d2 = sgn(poly[(i+1)%n].y - p.y);
        if(k > 0 && d1 <= 0 && d2 > 0) wn++;
        if(k < 0 && d2 <= 0 && d1 > 0) wn--;
    }
    if(wn != 0)
        return 1;
    return 0;
}

6 圆

计算机中储存圆通常记录圆心坐标与半径即可

6.1 定义

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);
    }
};

6.2 常用函数

6.2.1 圆与直线交点

//求圆与直线交点
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(sgn(delta) < 0)//相离
        return 0;
    if(sgn(delta) == 0){//相切
        t1 = -f /(2*e);
        t2 = -f /(2*e);
        sol.push_back(L.point(t1));//sol存放交点本身
        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;
}

6.2.2 求两圆交点

6.2.3 点到圆的切线

6.2.4 两圆的公切线

6.2.5 两圆相交面积

通过计算两个圆相交所构成的两个扇形面积和减去其构成的筝形的面积

double AreaOfOverlap(Point c1, double r1, Point c2, double r2){
    double d = Length(c1 - c2);
    if(r1 + r2 < d + eps)
        return 0.0;
    if(d < fabs(r1 - r2) + eps){
        double r = min(r1, r2);
        return pi*r*r;
    }
    double x = (d*d + r1*r1 - r2*r2)/(2.0*d);
    double p = (r1 + r2 + d)/2.0;
    double t1 = acos(x/r1);
    double t2 = acos((d - x)/r2);
    double s1 = r1*r1*t1;
    double s2 = r2*r2*t2;
    double s3 = 2*sqrt(p*(p - r1)*(p - r2)*(p - d));
    return s1 + s2 - s3;
}