计算几何模板 | Kyons'Blog 


 




计算几何模板

发表于 | 分类于 | | 阅读次数:

计算几何模板

正如不知何方大佬所言,计算几何精妙之处,就是不用解析几何的方法去做
为了方便查找,防止自己迷路,我把函数名都写成了拼音
绝对不是因为我英语不好!!!

基本数据结构

点和向量:

  • 点和向量都可以用一个坐标$(x,y)$来表示.
  • 故向量$Vector$可以写为
  • 1
    typedef struct point vec;

完整定义如下

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
typedef struct point vec;  //向量vec
struct point {    //点的基本数据结构
    double x, y;
    point(double _x=0, double _y=0):
    x(_x),y(_y)
    {
    }
    point operator*(const point& i_T) const
    {
        return point(x * i_T.x, y * i_T.y);
    }
    point operator*(double u) const
    {
        return point(x * u, y * u);
    }
    bool operator==(const point& i_T) const
    {
        return x == i_T.x && y == i_T.y;
    }
    point operator/(double u) const
    {
        return point(x / u, y / u);
    }
    point operator+(const point& i_T)
    {
        return point(x + i_T.x, y + i_T.y);
    }
    point operator-(const point& i_T)
    {
        return point(x - i_T.x, y - i_T.y);
    }
    friend bool operator<(point a, point b)
    {
        return a.y == b.y ? a.x < b.x : a.y < b.y;
    }
    friend ostream& operator<<(ostream& out, point& a)
    {
        //cout << a.x << ' ' << a.y;
        printf("%.8f %.8f", a.x, a.y);
        return out;
    }
    friend istream& operator>>(istream& in, point& a)
    {
        in >> a.x >> a.y;
        return in;
    }
};

直线和线段:

  • 直线和线段都可以用两个点的坐标来表示
  • 故线段$Segment$可以写为
  • 1
    typedef struct Line Segment;

完整定义如下

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
typedef struct Line Segment;   //线段Segment
struct Line {       //直线
    vec a, b;
    Line(point _a = point(), point _b = point())
        : a(_a)
        , b(_b)
    {
    }
    friend istream& operator>>(istream& in, Line& a)
    {
        cin >> a.a >> a.b;
        return in;
    }
    friend ostream& operator<<(ostream& out, Line& a)
    {
        out << a.a << ' ' << a.b;
        return out;
    }
};

圆:

  • 圆的表示有一个点圆心,以及其半径组成

完整定义如下

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
struct cirles {
    point o;
    double r;
    cirles(point _o = point(), double _r = 0.0)
        : r(_r)
        , o(_o)
    {
    }
    point Point(double t) //圆上任意一点
    {
        return point(o.x + r * cos(t), o.y + r * sin(t));
    }
    friend istream& operator>>(istream& in, cirles& a)
    {
        in >> a.o >> a.r;
        return in;
    }
    friend ostream& operator<<(ostream& out, cirles& a)
    {
        out << a.o << ' ' << a.r;
        return out;
    }
};

基本函数以及常量

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
const double pi = acos(-1.0);
const double eps = 1e-12;
int zhengfu(double d)   //判断正负,即sign/dcmp
{
    if (fabs(d) < eps)
        return 0;
    if (d > 0)
        return 1;
    return -1;
}
int bijiao(double x, double y)  //判断x和y的大小关系
{
    if (fabs(x - y) < eps)
        return 0;
    if (x > y)
        return 1;
    return -1;
}
int dayu_dengyu(double x, double y) //x>=y否?
{
    if (fabs(x - y) < eps || x > y)
        return 1;
    return 0;
}
double chaji(vec A, vec B)  //求向量叉积,即cross
{
    return A.x * B.y - A.y * B.x; // 正为A->B左旋
}
double xuanzhuan(point a, point b, point c)   //求三点叉积,即side
{
    return chaji(b - a, c - a);
}
bool cmp(vec a, vec b) //极角排序,运用叉积的极角排序,相比于atan2慢,但是精度高
{
    vec c(0, 0);
    if (chaji(a - c, b - c) == 0)
        return a.x < b.x;
    return chaji(a - c, b - c) > 0;
}
double dianji(vec A, vec B) //点积
{
    return A.x * B.x + A.y * B.y;
}
double changdu(vec a) //向量长度
{
    return sqrt(a.x * a.x + a.y * a.y);
}

AOJ相关习题

CGL_1_A:Projection

求一个点在向量$\overrightarrow{ab}$上的投影坐标
设点$c$,投影在$\overrightarrow{ab}$上为$c’$,则$c’$的坐标就是:$cos<\overrightarrow{ac},\overrightarrow{ab}>\times |\overrightarrow{ac}|+a$

代码内容
1
2
3
4
5
6
7
8
point touying(Line l, point c) //c投影在直线ab上的位置
{
    vec A = l.b - l.a;
    vec B = c - l.a;
    double La = changdu(A);
    double Lc = dianji(A, B) / (La * La);
    return A * Lc + l.a;
}

CGL_1_B:Reflection

求一个点$c$关于向量$\overrightarrow{ab}$的对称点$c’’$
先求出$c$在$ab$上的投影,那么$c’’=2\times \overrightarrow{cc’}+c$

代码内容
1
2
3
4
5
point fanshe(Line l, point c) //求c关于直线ab的对称点c'
{
    point A = touying(l, c);
    return (A - c) * 2.0 + c;
}

CGL_1_C:Counter-Clockwise

就是..根据图中的判断就是了

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
void Counter_Clockwise(point p,point p1,point p2)
{
    if (zhengfu(chaji(p2 - p1, p - p1)) == 1)
            cout << "COUNTER_CLOCKWISE" << endl;
        else if (zhengfu(chaji(p2 - p1, p - p1)) == -1)
            cout << "CLOCKWISE" << endl;
        else {
            if (zhengfu(dianji(p2 - p1, p - p1)) == -1)
                cout << "ONLINE_BACK" << endl;
            else {
                double j = changdu(p2 - p1);
                double k = changdu(p - p1);
                if (bijiao(j, k) >= 0)
                    cout << "ON_SEGMENT" << endl;
                else
                    cout << "ONLINE_FRONT" << endl;
            }
        }
}

CGL_2_A:Parallel/Orthogonal

先判平行,再用点积判垂直

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
bool pingxing(vec a,vec b)
{
    return bijiao(a.x*b.y,a.y*b.x)==0;
}
void Parallel/Orthogonal(Line l1,Line l2)
{
    if (pingxing(l1.b-l1.a,l2.b-l2.a))
            cout << "2" << endl;
    else {
        if (zhengfu(dianji(l1.b-l1.a,l2.b-l2.a))==0)
            cout << "1" << endl;
        else {
            cout << "0" << endl;
        }
    }
}

CGL_2_B:Intersection

线段相交要考虑蛮多的,首先,先按照x后y从小到大排一下.
最简单的情况,$ab$穿过$cd$,那么必定有交点.
第二种,$a$在$cd$上或者$b$在$cd$上
第三种,共线时,$a$在$cd$之间或$b$在$cd$之间.
处理好以上问题,就解决了

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
bool bijiao3(vec a, vec b, vec c)
{
    if (a.x <= c.x && c.x <= b.x && ((a.y <= c.y && c.y <= b.y || b.y <= c.y && c.y <= a.y)))
        return 1;
    return 0;
}
bool xianduan_xiangjiao(Segment l1,Segment l2) //两线段是否有交点
{
    if (l1.a.x == l1.b.x) {
        if (l1.a.y > l1.b.y)
            swap(l1.a, l1.b);
    } else if (l1.a.x > l1.b.x)
        swap(l1.a, l1.b);
    if (l2.a.x == l2.b.x) {
        if (l2.a.y > l2.b.y)
            swap(l2.a, l2.b);
    } else if (l2.a.x > l2.b.x)
        swap(l2.a, l2.b);

    double c1 = chaji(l1.b - l1.a, l2.a - l1.a), d1 = chaji(l1.b - l1.a, l2.b - l1.a);
    double c2 = chaji(l2.b - l2.a, l1.a - l2.a), d2 = chaji(l2.b - l2.a, l1.b - l2.a);
    
    if (zhengfu(c1 * d1) < 0 && zhengfu(c2 * d2) < 0) //ab横穿cd
        return 1;
    else if (zhengfu(c1 * d1) != 0 && zhengfu(c2 * d2) == 0) { //ab不穿过cd
        if (zhengfu(c2) == 0) {
            if (bijiao3(l2.a, l2.b, l1.a))
                return 1;
        }
        if (zhengfu(d2) == 0)
            if (bijiao3(l2.a, l2.b, l1.b))
                return 1;
    } else if (zhengfu(c1 * d1) == 0 && zhengfu(c2 * d2) != 0) { //cd不穿过ab
        if (c1 == 0)
            if (bijiao3(l1.a, l1.b, l2.a))
                return 1;
        if (d1 == 0)
            if (bijiao3(l1.a, l1.b, l2.b))
                return 1;
    } else if (zhengfu(c1 * d1) == 0 && zhengfu(c2 * d2) == 0) { //平行
        if (l1.a == l2.a || l1.a == l2.b || l1.b == l2.a || l1.b == l2.b)
            return 1;
        if (bijiao3(l1.a, l1.b, l2.a) == 1)
            return 1;
        if (bijiao3(l2.a, l2.b, l1.a) == 1)
            return 1;
    }
    return 0;
}

CGL_2_C:Cross Point

给你两个必定相交的线段,求交点

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
point xianduan_jiaodian(Segment l1,Segment l2)//两线段交点
{ 
    double tmpLeft, tmpRight, x = inf, y = inf;
    if (xianduan_xiangjiao(l1,l2)) {
        tmpLeft = (l2.b.x - l2.a.x) * (l1.a.y - l1.b.y) - (l1.b.x - l1.a.x) * (l2.a.y - l2.b.y);
        tmpRight = (l1.a.y - l2.a.y) * (l1.b.x - l1.a.x) * (l2.b.x - l2.a.x) + l2.a.x * (l2.b.y - l2.a.y) * (l1.b.x - l1.a.x) - l1.a.x * (l1.b.y - l1.a.y) * (l2.b.x - l2.a.x);

        x = tmpRight / tmpLeft;

        tmpLeft = (l1.a.x - l1.b.x) * (l2.b.y - l2.a.y) - (l1.b.y - l1.a.y) * (l2.a.x - l2.b.x);
        tmpRight = l1.b.y * (l1.a.x - l1.b.x) * (l2.b.y - l2.a.y) + (l2.b.x - l1.b.x) * (l2.b.y - l2.a.y) * (l1.a.y - l1.b.y) - l2.b.y * (l2.a.x - l2.b.x) * (l1.b.y - l1.a.y);
        y = tmpRight / tmpLeft;
    }
    return point(x, y);
}

CGL_2_D:Distance

  • 给定两个不相交线段,求两个线段最近距离
  • 很明显,最近距离就是两个端点到另一个线段的距离.
  • 那么两遍点到线段距离就出来了.
    • 点到线段距离有三种
    • 第一种是点在线段正上方,则距离为过点向线段作垂线
    • 第二种是点在左侧,就是左端点和该点连线
    • 第三种同第二种,不过在右侧
代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
double dian_dao_xianduan(Segment l, point c) //点到线段的距离
{
    double L = changdu(l.b - l.a);
    double r = dianji(l.b - l.a, c - l.a) / (L * L);

    if (zhengfu(r) == -1) {
        return changdu(c - l.a);
    } else if (dayu_dengyu(r, 1)) {
        return changdu(c - l.b);
    } else {
        double L = r * changdu(l.b - l.a), l2 = changdu(c - l.a);
        return sqrt(l2 * l2 - L * L);
    }
}
double xianduanjuli(Segment l1,Segment l2) //两线段距离
{
    if (xianduan_xiangjiao(l1,l2))
        return 0.0;
    double minn = inf;
    double l = dian_dao_xianduan(l1, l2.a);
    minn = dayu_dengyu(minn, l) ? l : minn;
    l = dian_dao_xianduan(l1, l2.b);
    minn = dayu_dengyu(minn, l) ? l : minn;
    l = dian_dao_xianduan(l2, l1.a);
    minn = dayu_dengyu(minn, l) ? l : minn;
    l = dian_dao_xianduan(l2, l1.b);
    minn = dayu_dengyu(minn, l) ? l : minn;
    return minn;
}

CGL_3_A:Area

计算多边形面积的方法蛮多的.
最暴力的当属以原点和多边形临近两点构成三角形,然后计算三角形的有向面积.
多边形内外符号不同,最后留下的就是多边形面积,然后fabs一下就完事了.这个地方建议”脑洞大开”或者拿纸画画.
不过要是会三角剖分的话,把多边形按顶点分割成一堆三角形,然后求面积也阔以.

代码内容
1
2
3
4
5
6
7
8
double duobianxingmianji(int n) //多边形面积
{
    double ans = 0;
    for (int i = 1; i <= n; i++)
        ans += chaji(p[i], p[(i + 1) % n]);
    ans = fabs(ans) * 0.5;
    return ans;
}

CGL_3_B:Is-Convex

问所给的多边形是不是凸的.
题目给的方法是计算内角和.
嘛,感觉好麻烦的亚子.
还不如暴力求个凸包,看看所给的多边形的点数是不是和凸包点数相同来的快.

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
void Andrew(int& tail)  //求凸包
{
    sort(p + 1, p + 1 + n);
    tail = 1;
    q[1] = p[1];
    for (int i = 2; i <= n; i++) {
        while (tail > 1 && xuanzhuan(q[tail - 1], q[tail], p[i]) < 0)
            tail--;
        q[++tail] = p[i];
    }
    int basic = tail;
    for (int i = n - 1; i >= 1; i--) {
        while (tail > basic && xuanzhuan(q[tail - 1], q[tail], p[i]) < 0)
            tail--;
        q[++tail] = p[i];
    }
}

CGL_3_C:Polygon-Point Containment

就,判断点和多边形的位置关系.
看网上都是角度和或者射线法.
结果就让我找到一个看起来很$nb$的象限角度法?
不用考虑角度的精度问题,还不用像射线法考虑多??

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
bool zaibianshang(point& t,int n) //点在多边形边上否
{
    for (int i = 1; i <=n; i++)
        if (dian_zai_xianshang(Line(p[i], p[(i + 1)%n]), t))
            return 1;
    return 0;
}
bool duobianxingnei(point& t,int n) //点在多边形内
{
    int t1, t2, sum, i;
    double f;
    p[0] = p[n];
    for (i = 0; i <= n; i++)
        p[i] = p[i] - t; // 坐标平移
    t1 = p[0].x >= 0 ? (p[0].y >= 0 ? 0 : 3) : (p[0].y >= 0 ? 1 : 2); // 计算象限

    for (sum = 0, i = 1; i <= n; i++) {
        if (fabs(p[i].x)<eps && fabs(p[i].y)<eps)
            break;
        // 被测点为多边形顶点
        f = chaji(p[i - 1], p[i - 1]);

        // 计算叉积
        if (fabs(f)<eps && p[i - 1].x * p[i].x <= 0 && p[i - 1].y * p[i].y <= 0)
            break; // 点在边上
        t2 = p[i].x >= 0 ? (p[i].y >= 0 ? 0 : 3) : (p[i].y >= 0 ? 1 : 2); // 计算象限
        if (t2 == (t1 + 1) % 4)
            sum += 1;
        // 情况1
        else if (t2 == (t1 + 3) % 4)
            sum -= 1;
        // 情况2
        else if (t2 == (t1 + 2) % 4)
        // 情况3
            if (f > 0)
                sum += 2;
            else
                sum -= 2;
        t1 = t2;
    }
    bool tf = 0;
    if (i <= n || sum)
        tf = 1;
    for (i = 0; i <= n; i++)
        p[i] = p[i] + t; // 恢复坐标
    return tf;
}

CGL_4_A:Convex Hull

就是…求凸包
这里有个蛋疼的地方,要求是先排y再排x

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
void Andrew(int& tail)
{
    sort(p + 1, p + 1 + n);
    tail = 1;
    q[1] = p[1];
    for (int i = 2; i <= n; i++) {
        while (tail > 1 && xuanzhuan(q[tail - 1], q[tail], p[i]) < 0)
            tail--;
        q[++tail] = p[i];
    }
    int basic = tail;
    for (int i = n - 1; i >= 1; i--) {
        while (tail > basic && xuanzhuan(q[tail - 1], q[tail], p[i]) < 0)
            tail--;
        q[++tail] = p[i];
    }
}

CGL_4_B:Diameter of a Convex Polygon

找到凸包距离最远的一对点.
就是旋转卡壳嘛.$O(n)$复杂度嘛.

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
double tubao_zhijing(int tail) //求出凸包直径
{
    double re = 0;
    if (tail == 2) //仅有两个点
        return changdu(q[2] - q[1]);
    q[0] = q[tail]; //把最后的点放到最前面
    for (int i = 0, j = 2; i < tail; i++) //枚举边
    {
        while (xuanzhuan(q[i], q[i + 1], q[j]) < xuanzhuan(q[i], q[i + 1], q[j + 1]))
            j = (j + 1) % tail;
        re = max(re, max(changdu(q[j] - q[i]), changdu(q[j] - q[i + 1])));
    }
    return re;
}

CGL_4_C:Convex Cut

用一条直线切割凸包,输出得到图形的坐标.
就是逆时针找交点,按照直线的方向,$\overrightarrow{p_1p_2}$,先放入靠近$p_2$的点,然后按照叉积,向左旋转的放入点.

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
int zhixian_xianduan_xiangjiao(Line l1, Segment l2) //直线与线段是否有交点
{
    double c1 = chaji(l1.b - l1.a, l2.a - l1.a), d1 = chaji(l1.b - l1.a, l2.b - l1.a);
    if (zhengfu(c1) == 0 && zhengfu(d1) == 0) { //重合
        return -1;
    } else if (zhengfu(c1 * d1) <= 0) //有交点
        return 1;
    return 0; //平行
}
void qiegetubao(Line l)
{
    int tail = 0;
    q[0] = q[n];
    int tf = -1, x;
    for (int i = 0; i < n; i++) {
        if (zhengfu(xuanzhuan(l.a,l.b,q[i])) == 1)
            p[++tail] = q[i];
        int f = zhixian_xianduan_xiangjiao(l,Segment( q[i], q[i + 1]));

        if (f == 1) {
            tf = 1;
            p[++tail] = zhixian_xianduan_jiaodian(l, Segment(q[i], q[i + 1]));
        } else if (f == -1) {
            tf = 0;
            x = i;
            break;
        }
    }
    if (tf == 0) {
        //cout<<q[x+1]<<' '<<q[x]<<' '<<endl;
        if (zhengfu(dianji(l.b - l.a, q[x + 1] - q[x])) == -1)
            printf("%.8f\n", 0.0);
        else {
            for (int i = 1; i <= n; i++)
                p[i] = q[i];
            printf("%.8f\n", duobianxingmianji(n));
        }
    } else {
        /*for (int i = 1; i <= tail; i++)
            cout << p[i] << endl;*/
        printf("%.8f\n", duobianxingmianji(tail));
    }
    q[0] = point(0, 0);
}

CGL_5_A:Closest Pair

在空间内找到最近的点对
分治,建议到洛谷上搜索”平面最近点对(加强版)”

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
double solve(int l, int r)
{
    if (l == r)
        return inf;
    int mid = (l + r) >> 1;

    double ans = solve(l, mid);
    ans = min(ans, solve(mid + 1, r));

    int tot = 0;
    for (int i = l; i <= r; i++)
        if (fabs(p[mid].x - p[i].x) <= ans)
            temp[tot++] = p[i];

    sort(temp,temp+tot, cmp);

    for (int i = 0; i < tot; i++)
        for (int j = i + 1; j <tot; j++) {
            if (temp[j].y - temp[i].y > ans)
                break;
            ans = min(ans, changdu(temp[j] - temp[i]));
        }
    return ans;
}

CGL_6_A:Segment Intersections: Manhattan Geometry

扫描线算法,例题可在洛谷上搜索”扫描线”
求平面$n$条线段的交点个数.
有空单开一章来整理这个算法.

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
struct SegTree {
    int l, r;
    ll len;
    //  sum: 被完全覆盖的次数；
    //  len: 区间内被截的长度。
} tree[MAXN << 2];
 
void build_tree(int x, int l, int r)
{
    tree[x].l = l, tree[x].r = r;
    tree[x].len = 0;
    if (l == r)
        return;
    int mid = (l + r) >> 1;
    build_tree(lson(x), l, mid);
    build_tree(rson(x), mid + 1, r);
    return;
}
 
void pushup(int x)
{
    tree[x].len = tree[x << 1].len + tree[x << 1 | 1].len;
    //      合并儿子信息
}
 
void edit_tree(int x, ll id, int c)
{
    int l = tree[x].l, r = tree[x].r;
    if (l == id && id == r) {
        tree[x].len += c;
        return;
    }
    int mid = (l + r) >> 1;
    if (id <= mid)
        edit_tree(lson(x), id, c);
    else
        edit_tree(rson(x), id, c);
    pushup(x);
}
 
int query(int x, int L, int R)
{
    int l = tree[x].l, r = tree[x].r;
    if (L <= l && r <= R)
        return tree[x].len;
    int mid = (l + r) >> 1;
    int res = 0;
    if (L <= mid)
        res += query(lson(x), L, R);
    if (R > mid)
        res += query(rson(x), L, R);
    return res;
}
vector<Line> l;
vector<int> X;
int main()
{
    scanf("%d", &n);
    for (int i = 1; i <= n; i++) {
        point a, b;
        cin >> a >> b;
        if (a.x != b.x) {
            if (a.x > b.x)
                swap(a, b);
            X.push_back(a.x), X.push_back(b.x);
            l.push_back(Line(a, b, 2));
        } else {
            if (a.y > b.y)
                swap(a, b);
            X.push_back(a.x);
            l.push_back(Line(a, a, 1));
            l.push_back(Line(b, b, 3));
        }
    }
     
    sort(l.begin(), l.end());
    sort(X.begin(), X.end());
 
    n = unique(X.begin(), X.end()) - X.begin(); //去重
    X.erase(X.begin() + n, X.end());
    if(n==1){
        cout<<0<<endl;
        return 0;
    }
    build_tree(1, 1, n); //根据X[]的坐标建立线段树
 
    int ans = 0;
    for (int i = 0; i < l.size(); i++) {
        if (l[i].mark!=2) {
            int id = lower_bound(X.begin(), X.end(), l[i].a.x) - X.begin();
            int c = l[i].mark == 3 ? -1 : 1;
            edit_tree(1, id+1, c);
        } else {
            int id1 = lower_bound(X.begin(), X.end(), l[i].a.x) - X.begin();
            int id2 = lower_bound(X.begin(), X.end(), l[i].b.x) - X.begin();
            ans += query(1, id1+1, id2+1);
        }
    }
    printf("%d\n", ans);
    return 0;
}

CGL_7_A:Intersection

求圆的切线个数

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
int yuan_yuan_xiangjiao(cirles a, cirles b) //询问圆和圆的切线个数
{
    double l = changdu(b.o - a.o);
    if (bijiao(l, a.r + b.r) == 1)
        return 4;
    else if (bijiao(l, a.r + b.r) == 0)
        return 3;
    else {
        if (bijiao(l + min(a.r, b.r), max(a.r, b.r)) == 1)
            return 2;
        else if (bijiao(l + min(a.r, b.r), max(a.r, b.r)) == 0)
            return 1;
        else
            return 0;
    }
}

CGL_7_D:Cross Points of a Circle and a Line

求圆和直线的交点.
建议阅读”挑战程序设计竞赛2”

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
pair<point, point> zhixian_yuan_jiaodian(Line l, cirles a) //求直线与圆交点
{
    if (dayu_dengyu(a.r, zhixian_yuanxin_juli(l, a))) {
        vec a_i = touying(l, a.o);
        double L = changdu(l.b - l.a);
        double lt = changdu(a.o - a_i);
        double lr = sqrt(a.r * a.r - lt * lt);
        vec p = (l.b - l.a) / L * lr + a_i;
        return make_pair(a_i - (l.b - l.a) / L * lr, a_i + (l.b - l.a) / L * lr);
    } else
        return make_pair(0, 0);
}

CGL_7_E:Cross Points of Circles

求两个圆交点
建议阅读”挑战程序设计竞赛2”

代码内容
1
2
3
4
5
6
7
pair<vec, vec> yuan_yuan_jiaodian(cirles a, cirles b) //求圆和圆的交点
{
    double l = changdu(a.o - b.o);
    double x = acos((a.r * a.r + l * l - b.r * b.r) / (2.0 * a.r * l));
    double t = atan2((b.o - a.o).y, (b.o - a.o).x);
    return make_pair(a.o + vec(cos(t - x) * a.r, sin(t - x) * a.r), a.o + vec(cos(x + t) * a.r, sin(t + x) * a.r));
}

CGL_7_F:Tangent to a Circle

过一个点做圆的切线
根据半径和圆心到点的距离求出夹角,旋转角度,得到交点

代码内容
1
2
3
4
5
6
7
8
9
pair<point, point> guodian_yuan_qiedian(cirles a, point p) //过一点做圆的切线求切点
{
    double l = changdu(a.o - p);
    double t = asin(a.r / l);
    double lb = l * cos(t);
    vec x = a.o - p;
    x = x / l * lb;
    return make_pair(p + xiangliang_xuanzhuan(x, t), p + xiangliang_xuanzhuan(x, -t));
}

CGL_7_G:Common Tangent

求两个圆的公切线
根据圆和圆的位置进行判断

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
int yuanyuan_gongqiexian(cirles a, cirles b, point* u, point* v) //求圆和圆公切线以及切线个数
{
    int cnt = 0;
    if (a.r < b.r) {
        swap(a, b);
        swap(u, v);
    }
    double l = changdu(a.o - b.o);
    double rdiff = a.r - b.r;
    double rsum = a.r + b.r;
    if (zhengfu(l - rdiff) < 0)
        return 0;
    double base = atan2((b.o - a.o).y, (b.o - a.o).x);
    if (zhengfu(l) == 0)
        return -1;
    if (zhengfu(l - rdiff) == 0) {
        u[cnt] = v[cnt] = a.Point(base);
        cnt++;
        return 1;
    }
    double ang = acos((a.r - b.r) / l);
    u[cnt] = a.Point(base + ang);
    v[cnt] = b.Point(base + ang);
    cnt++;
    u[cnt] = a.Point(base - ang);
    v[cnt] = b.Point(base - ang);
    cnt++;
    if (zhengfu(l - rsum) == 0) {
        u[cnt] = v[cnt] = a.Point(base);
        cnt++;
    } else if (zhengfu(l - rsum) > 0) {
        double ang = acos((a.r + b.r) / l);
        u[cnt] = a.Point(base + ang);
        v[cnt] = b.Point(pi + base + ang);
        cnt++;
        u[cnt] = a.Point(base - ang);
        v[cnt] = b.Point(pi + base - ang);
        cnt++;
    }
    return cnt;
}

CGL_7_H:Intersection of a Circle and a Polygon

求多边形和圆相交的面积
将多边形的边的顶点与圆心连接行成三角形.
那么面积便是三角形在圆内的有向面积.
对每个三角形在圆内进行判断来计算.

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
point zhixian_zhixian_jiaodian(Line l1, Line l2) //两直线交点
{
    double t = ((l1.a.x - l2.a.x) * (l2.a.y - l2.b.y) - (l1.a.y - l2.a.y) * (l2.a.x - l2.b.x)) / ((l1.a.x - l1.b.x) * (l2.a.y - l2.b.y) - (l1.a.y - l1.b.y) * (l2.a.x - l2.b.x));
    return l1.a + (l1.b - l1.a) * t;
}

point xianduan_duandian_dian(point p, Segment l) //线段距离点p最近的端点
{
    point t = p;
    t.x += l.a.y - l.b.y;
    t.y += l.b.x - l.a.x;
    if (chaji(l.a - p, t - p) * chaji(l.b - p, t - p) > eps)
        return changdu(p - l.a) < changdu(p - l.b) ? l.a : l.b;
    return zhixian_zhixian_jiaodian(Line(p, t), l);
}

double distp(Line l) //长度的平方
{
    return (l.a.x - l.b.x) * (l.a.x - l.b.x) + (l.a.y - l.b.y) * (l.a.y - l.b.y);
}
double yuanxin_dian_sanjiao(Line l, cirles c) //求圆心与两点所成三角形有向面积
{
    double sign = 1.0;
    l.a = l.a - c.o;
    l.b = l.b - c.o;
    c.o = point(0.0, 0.0);
    if (fabs(chaji(l.a - c.o, l.b - c.o)) < eps)
        return 0.0;
    if (distp(Line(l.a, c.o)) > distp(Line(l.b, c.o))) {
        swap(l.a, l.b);
        sign = -1.0;
    }
    if (distp(Line(l.a, c.o)) < c.r * c.r + eps) { //a在圆内
        if (distp(Line(l.b, c.o)) < c.r * c.r + eps) //b也在圆内,返回叉积/2
            return chaji(l.a - c.o, l.b - c.o) / 2.0 * sign;
        point p1, p2;
        pair<point, point> q = zhixian_yuan_jiaodian(l, c); //oa和ob与圆的交点
        p1 = q.first;
        p2 = q.second;
        if (changdu(p1 - l.b) > changdu(p2 - l.b))
            swap(p1, p2);
        double ret1 = fabs(chaji(l.a - c.o, p1 - c.o));
        double ret2 = acos((p1.x * l.b.x + p1.y * l.b.y) / changdu(p1) / changdu(l.b)) * c.r * c.r;
        double ret = (ret1 + ret2) / 2.0;
        if (chaji(l.a - c.o, l.b - c.o) < eps && sign > 0.0 || chaji(l.a - c.o, l.b - c.o) > eps && sign < 0.0)
            ret = -ret;
        return ret;
    }
    point ins = xianduan_duandian_dian(c.o, l);
    if (distp(Line(c.o, ins)) > c.r * c.r - eps) {
        double ret = acos((l.a.x * l.b.x + l.a.y * l.b.y) / changdu(l.a) / changdu(l.b)) * c.r * c.r / 2.0;
        if (chaji(l.a - c.o, l.b - c.o) < eps && sign > 0.0 || chaji(l.a - c.o, l.b - c.o) > eps && sign < 0.0)
            ret = -ret;
        return ret;
    }
    point p1, p2;
    pair<point, point> q = zhixian_yuan_jiaodian(l, c); //oa和ob与圆的交点
    p1 = q.first;
    p2 = q.second;
    double cm = c.r / (changdu(c.o - l.a) - c.r);
    point m = point((c.o.x + cm * l.a.x) / (1 + cm), (c.o.y + cm * l.a.y) / (1 + cm));
    double cn = c.r / (changdu(c.o - l.b) - c.r);
    point n = point((c.o.x + cn * l.b.x) / (1 + cn), (c.o.y + cn * l.b.y) / (1 + cn));
    double ret1 = acos((m.x * n.x + m.y * n.y) / changdu(m) / changdu(n)) * c.r * c.r;
    double ret2 = acos((p1.x * p2.x + p1.y * p2.y) / changdu(p1) / changdu(p2)) * c.r * c.r - fabs(chaji(p1 - c.o, p2 - c.o));
    double ret = (ret1 - ret2) / 2.0;
    if (chaji(l.a - c.o, l.b - c.o) < eps && sign > 0.0 || chaji(l.a - c.o, l.b - c.o) > eps && sign < 0.0)
        ret = -ret;
    return ret;
}
double duobianxing_yuan_xiangjiao(cirles c, point p[], int n) //多边形与圆相交面积
{
    double sum = 0;
    for (int i = 0; i < n; i++)
        sum += yuanxin_dian_sanjiao(Line(p[i], p[i + 1]), c);
    return sum;
}

其他相关

辛普森法则

证明待补

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
double fun(double x)//原积分函数
{
    return x;
}
double simpson(double l,double r)   //辛普森法则
{
    double mid=(l+r)/2.0;
    return (fun(l)+fun(r)+4*fun(mid))*(r-l)/6.0;
}
double solve(double l,double r,double ans)  //调整精度求答案
{ 
    double mid=(l+r)/2.0;
    double ls=simpson(l,mid),rs=simpson(mid,r);
    if(fabs(ls+rs-ans)<=15.0*eps)
        return ls+rs+(ls+rs-ans)/15.0;
    return solve(l,mid,ls)+solve(mid,r,rs);
}

最小圆覆盖

代码内容
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
point zhixian_zhixian_jiaodian(Line l1, Line l2) //两直线交点
{
    double t = ((l1.a.x - l2.a.x) * (l2.a.y - l2.b.y) - (l1.a.y - l2.a.y) * (l2.a.x - l2.b.x)) / ((l1.a.x - l1.b.x) * (l2.a.y - l2.b.y) - (l1.a.y - l1.b.y) * (l2.a.x - l2.b.x));
    return l1.a + (l1.b - l1.a) * t;
}
point sanjiaoxing_waixin(point a, point b, point c) //三角形外心
{
    Line u, v;
    u.a.x = (a.x + b.x) / 2;
    u.a.y = (a.y + b.y) / 2;
    u.b.x = u.a.x - a.y + b.y;
    u.b.y = u.a.y + a.x - b.x;
    v.a.x = (a.x + c.x) / 2;
    v.a.y = (a.y + c.y) / 2;
    v.b.x = v.a.x - a.y + c.y;
    v.b.y = v.a.y + a.x - c.x;
    return zhixian_zhixian_jiaodian(u, v);
}
cirles zuixiaoyuan_fugai(point p[], int n)  //最小圆覆盖
{
    random_shuffle(p + 1, p + 1 + n);
    cirles c(p[1], 0.0);
    for (int i = 2; i <= n; i++)
        if (zhengfu(changdu(c.o - p[i]) - c.r) > 0) {
            c = cirles(p[i], 0.0);
            for (int j = 1; j < i; j++)
                if (zhengfu(changdu(c.o - p[j]) - c.r) > 0) {
                    c.o = (p[i] + p[j]) / 2;
                    c.r = changdu(c.o - p[i]);
                    for (int k = 1; k < j; k++)
                        if (zhengfu(changdu(c.o - p[k]) - c.r) > 0) {
                            c.o = sanjiaoxing_waixin(p[i], p[j], p[k]);
                            c.r = changdu(c.o - p[i]);
                        }
                }
        }
    return c;
}







 




 


 


 

 //删除Valine核心代码库外链调用 //调用刚下载的本地文件以加速加载速度 //这里改为从本地加载