计算几何模板

庆祝该板子大概 20k 了/se！！！

另外新增了已验证的动态半平面交（肝了 4 天啊（悲。

#include <bits/stdc++.h>
template <typename T> inline void ckmax(T &x, T y) { x = x > y ? x : y; }
template <typename T> inline void ckmin(T &x, T y) { x = x < y ? x : y; }
typedef long long ll;
typedef unsigned long long ull;

namespace Geometry {
std::mt19937 rnd((unsigned)time(NULL));
const double eps = 1e-6, pi = acos(-1.0), inf = 1e100;
inline double sign(double x) { return x > eps ? 1 : (x < -eps ? -1 : 0); }
struct Vector {
    double x, y;
    Vector() { x = y = 0; }
    Vector(const double &x, const double &y) : x(x), y(y) {}
    Vector(const Vector &t) : x(t.x), y(t.y) {}
    friend Vector operator-(const Vector &x, const Vector &y) { return Vector(x.x - y.x, x.y - y.y); }
    Vector &operator-=(const Vector &t) { x -= t.x, y -= t.y; return *this; }
    friend Vector operator+(const Vector &x, const Vector &y) { return Vector(x.x + y.x, x.y + y.y); }
    Vector &operator+=(const Vector &t) { x += t.x; y += t.y; return *this; }
    friend Vector operator/(const Vector &x, const double &y) { return Vector(x.x / y, x.y / y); }
    Vector &operator/=(const double t) { x /= t; y /= t; return *this; }
    friend Vector operator*(const Vector &x, const double &y) { return Vector(x.x * y, x.y * y); }
    friend Vector operator*(const double &y, const Vector &x) { return Vector(x.x * y, x.y * y); }
    Vector &operator*=(const double t) { x *= t; y *= t; return *this; }
    double mod() const { return hypot(x, y); }
    double mod2() const { return x * x + y * y; }
    friend double operator*(Vector a, Vector b) { return a.x * b.x + a.y * b.y; } //Dot Product
    friend double operator^(Vector a, Vector b) { return a.x * b.y - a.y * b.x; } //cross Product
    friend bool operator==(Vector a, Vector b) { return sign(a.x - b.x) == 0 && sign(a.y - b.y) == 0; }
    friend bool operator!=(Vector a, Vector b) { return sign(a.x - b.x) != 0 || sign(a.y - b.y) != 0; }
    friend bool operator<(Vector a, Vector b) { return sign(a.x - b.x) ? a.x < b.x : sign(a.y - b.y) < 0; }
    friend bool operator>(Vector a, Vector b) { return sign(a.x - b.x) ? a.x > b.x : sign(a.y - b.y) > 0; }
    friend bool operator||(Vector a, Vector b) { return sign(a ^ b) == 0; }
    Vector rotate(const double &theta) const {
        double c = cos(theta), s = sin(theta);
        return Vector(x * c + y * s, -x * s + y * c);
    }
    void read() { std::cin >> x >> y; }
    Vector normal() const { return Vector(-y, x); }
    double slope() const { return sign(y) ? (sign(x) ? y / x : (y > 0 ? inf : -inf)) : 0; } 
    double angle() const { return atan2(y, x); }
};
typedef Vector Point;
double dot(Vector a, Vector b) { return a.x * b.x + a.y * b.y; }
double cross(Vector a, Vector b) { return a.x * b.y - a.y * b.x; }
double angle(Vector a, Vector b) {
    return acos((a * b) / a.mod() / b.mod());
}
double dis(Point x, Point y) { return (x - y).mod(); }
//it's a line segment when seg = true
struct Line {
    Point A, B; bool seg;
    Line() { seg = false; }
    Line(const Point &A, const Point &B, bool seg = false) : A(A), B(B), seg(seg) {}
    Line(const Line &t) : A(t.A), B(t.B), seg(t.seg) {}
    Line change() const { return Line(A, B, seg ^ 1); }
    double mod() const{ return (B - A).mod(); }
    double mod2() const { return (B - A).mod2(); }
    Vector direct() const { return B - A; }
    double slope() const { return direct().slope(); }
    void read() { A.read(); B.read(); }
    bool left(const Point &x) const { return sign(cross(B - A, x - A)) > 0; }
    bool right(const Point &x) const { return sign(cross(B - A, x - A)) < 0; }
    friend bool operator||(const Line &x, const Line &y) {
        return (x.B - x.A) || (y.B - y.A);
    }
    friend bool operator<(const Line &x, const Line &y) {
        Vector a = x.B - x.A, b = y.B - y.A;
        double A = atan2(a.y, a.x), B = atan2(b.y, b.x);
        if (sign(A - B) == 0) return y.left(x.A);
        return A > B; 
    }
    double angle() const { return direct().angle(); }
};
bool cap(const Point &p, const Line &l) { return sign(cross(p - l.A, l.B - l.A)) == 0 && (!l.seg || sign(dot(p - l.A, p - l.B)) <= 0); }
bool cap(const Line &l, const Point &p) { return cap(p, l); }
Point footPoint(const Point &p, const Line &l) {
    Vector Ap = p - l.A, Bp = p - l.B, AB = l.B - l.A;
    double l1 = (Ap * AB) / AB.mod(), l2 = -(Bp * AB) / AB.mod();
    return l.A + AB * (l1 / (l1 + l2));
}
Point closestPoint(const Point &p, const Line &l) {
    Vector pA = l.A - p, pB = l.B - p, AB = l.B - l.A;
    if (sign(pA * AB) > 0) return l.A;
    if (sign(pB * AB) < 0) return l.B;
    return footPoint(p, l);
}
double dis(const Point &p, const Line &l) {
    Vector pA = l.A - p, pB = l.B - p, AB = l.B - l.A;
    if (l.seg) {
        if (sign(pA * AB) > 0) return pA.mod();
        if (sign(pB * AB) < 0) return pB.mod();
    }
    return fabs(pA ^ pB) / AB.mod();
}
Point symmetryPoint(const Point &p, const Line &l) {
    return (footPoint(p, l) * 2) - p;
}
Point intersection(const Line &x, const Line &y) {//ensure they aren't parallel
    Vector AB = y.A - x.A, v1 = x.B - x.A, v2 = y.B - y.A;
    return x.A + (((AB ^ v2) / (v1 ^ v2)) * v1);
}
int cap(const Line &x, const Line &y) {//updated
    if (x || y) {
        if (sign(cross(x.A - y.A, x.B - y.A)) == 0) return 0;
        if (x.seg && y.seg) {
            return cap(x.A, y) || cap(x.B, y) || cap(y.A, x) || cap(y.B, x) ? -1 : 0;
        }
        if (x.seg || y.seg) {
            return (x.seg && (cap(x.A, y) || cap(x.B, y))) || (y.seg && (cap(y.A, x) || cap(y.B, x))) ? -1 : 0;
        }
        return -1;
    }
    if (x.seg && y.seg) {
        if (!(std::min(x.A.x, x.B.x) <= std::max(y.A.x, y.B.x) && std::min(y.A.x, y.B.x) <= std::max(x.A.x, x.B.x) && std::min(x.A.y, x.B.y) <= std::max(y.A.y, y.B.y) && std::min(y.A.y, y.B.y) <= std::max(x.A.y, x.B.y))) {
            return false;
        }
        return sign(cross(x.A - y.A, y.B - y.A)) * sign(cross(x.B - y.A, y.B - y.A)) <= 0 && sign(cross(y.A - x.A, x.B - x.A)) * sign(cross(y.B - x.A, x.B - x.A)) <= 0;
    }
    if (x.seg || y.seg) {
        Point p = intersection(x, y);
        return x.seg ? cap(p, x) : cap(p, y);
    }
    return 1;
}

typedef std::vector <Point> Polygon;
// Polygon, clockwise 
double area(const Polygon &P) {
    double res = 0;
    int n = (int)P.size();
    for (int i = 0; i < n; ++i) {
        res += cross(P[(i + 1) % n], P[i]);
    }
    return res * 0.5;
}
int in(const Point &p, const Polygon &P) {
    int n = (int)P.size(), cnt = 0;
    for (int i = 0; i < n; ++i) { 
        int nx = (i + 1) % n;
        if (cap(p, Line(P[i], P[nx], true))) return 2;
        if (p.y >= std::min(P[i].y, P[nx].y) && p.y < std::max(P[i].y, P[nx].y)) {
            double tx = P[i].x + (P[nx].x - P[i].x) / (P[nx].y - P[i].y) * (p.y - P[i].y);
            cnt += sign(tx - p.x) > 0; 
        }
    }
    return cnt & 1;
}

// keep only the right part, clockwise
Polygon halfCut(std::vector <Line> L) {
    sort(L.begin(), L.end());
    std::vector <Line> q(L.size());
    auto judge = [&](const Line &l, const Point &p) {
        return l.left(p) || cap(p, l);
    };
    int head = 0, tail = -1;
    for (Line l : L) {
        while (head <= tail && (l || q[tail])) --tail;
        while (head < tail && judge(l, intersection(q[tail - 1], q[tail]))) --tail;
        while (head < tail && judge(l, intersection(q[head], q[head + 1]))) ++head;
        q[++tail] = l;
    }
    while (head < tail - 1 && judge(q[head], intersection(q[tail - 1], q[tail]))) --tail;
    Polygon res(tail - head + 1);
    for (int i = 0; i < tail - head + 1; ++i) {
        res[i] = intersection(q[i + head], q[i + head == tail ? head : i + head + 1]);
    }
    return res;
}

struct ConvexHull {
    Polygon U, D;//��͹�Ǻ���͹�ǡ� 
    ConvexHull() {}
    ConvexHull(const Polygon &U, const Polygon &D) : U(U), D(D) {}
    ConvexHull(const ConvexHull &H) : U(H.U), D(H. D) {}
    ConvexHull(const Polygon &con) { construct(con); }
    void construct(const Polygon &con) {
        U.clear(); D.clear();
        if (con.empty()) return;
        if ((int)con.size() == 1) {
            U = D = con;
            return;
        }
        int mn = 0, mx = 0, n = (int)con.size();
        for (int i = 1; i < n; ++i) {
            if (con[i] < con[mn]) mn = i;
            if (con[mx] < con[i]) mx = i;
        }
        U.emplace_back(con[mn]);
        for (int i = (mn + 1) % n; i != mx; i = (i + 1) % n) {
            U.emplace_back(con[i]);
        }
        U.emplace_back(con[mx]);
        D.emplace_back(con[mx]);
        for (int i = (mx + 1) % n; i != mn; i = (i + 1) % n) {
            D.emplace_back(con[i]);
        }
        D.emplace_back(con[mn]);
    }
    void clear() { U.clear(); D.clear(); }
    int size() const { return std::max((int)(U.size() + D.size()) - 2, 0); }
    Polygon merge() const {
        Polygon res(U);
        for (int i = 1; i < (int)D.size() - 1; ++i) {
            res.emplace_back(D[i]);
        }
        return res;
    }
    Point operator[](const int &i) const {
        if (i < (int)U.size()) return U[i];
        return D[i - U.size() + 1];
    }
};
// clockwise right
int inConvex(const Point &p, const ConvexHull &P) {
    int n = (int)P.size();
    if (n < 3) return 0;
    if (cap(p, Line(P[0], P[1], true)) || cap(p, Line(P[n - 1], P[0], true))) return 2;
    if (Line(P[0], P[1], true).left(p) || Line(P[n - 1], P[0], true).left(p)) return 0;
    int A = 2, B = n - 2, mid, pos = 1;
    while (A <= B) {
        mid = (A + B) >> 1;
        if (Line(P[0], P[mid], true).right(p)) pos = mid, A = mid + 1;
        else B = mid - 1;
    }
    if (cap(p, Line(P[pos], P[pos + 1], true))) return 2;
    if (Line(P[pos], P[pos + 1], true).left(p)) return 0;
    return 1;
}
Polygon getConvexHull(const Polygon &P) {
    Polygon res(0);
    int n = (int)P.size(), top = 2;
    res.emplace_back(P[0]);
    res.emplace_back(P[1]);
    for (int i = 2; i < n; ++i) {
        while (top > 1 && sign(cross(res[top - 1] - res[top - 2], P[i] - res[top - 1])) >= 0) {
            --top; res.pop_back();
        }
        res.emplace_back(P[i]); ++top;
    }
    return res;
};
ConvexHull andrew(Polygon P) {
    if ((int)P.size() < 2) return ConvexHull(P, P);
    ConvexHull H;
    sort(P.begin(), P.end());
    H.U = getConvexHull(P);
    reverse(P.begin(), P.end());
    H.D = getConvexHull(P);
    return H;
}
double radius(const ConvexHull &C) {
    Polygon T = C.merge();
    int n = (int)T.size();
    double ans = (T[0] - T[1]).mod();
    if (n == 2) return ans;
    for (int i = 0, j = 2; i < n; ++i) {
        while (sign(cross(T[i] - T[j], T[j] - T[(i + 1) % n]) - cross(T[i] - T[(j + 1) % n], T[(j + 1) % n] - T[(i + 1) % n])) < 0) {
            j = (j + 1) % n;
        }
        ckmax(ans, std::max(dis(T[j], T[i]), dis(T[j], T[(i + 1) % n])));
    }
    return ans;
}
Polygon mergeU(const Polygon &P, const Polygon &Q) {
    Polygon res(P.size() + Q.size());
    int now = 0;
    for (Point x : P) res[now++] = x;
    for (Point x : Q) res[now++] = x;
    inplace_merge(res.begin(), res.begin() + P.size(), res.end());
    return getConvexHull(res);
}
Polygon mergeD(const Polygon &P, const Polygon &Q) {
    Polygon res(P.size() + Q.size());
    int now = 0;
    for (Point x : P) res[now++] = x;
    for (Point x : Q) res[now++] = x;
    inplace_merge(res.begin(), res.begin() + P.size(), res.end(), [&](Point a, Point b) {
        return b < a;
    });
    return getConvexHull(res);
}
ConvexHull merge(const ConvexHull &P, const ConvexHull &Q) {
    return ConvexHull(mergeU(P.U, Q.U), mergeD(P.D, Q.D));
}
ConvexHull mincowski(const ConvexHull &F, const ConvexHull &G) {
    auto mincowski = [&](const Polygon &F, const Polygon &G) {
        int n = (int)F.size(), m = (int)G.size();
        std::vector <Vector> V1(n - 1), V2(m - 1);
        for (int i = 0; i < n - 1; ++i) V1[i] = F[i + 1] - F[i];
        for (int i = 0; i < m - 1; ++i) V2[i] = G[i + 1] - G[i];
        int i, j = 0;
        Polygon res;
        res.emplace_back(F[0] + G[0]);
        for (i = 0; i < n - 1; ++i) {
            while (j < m - 1 && sign(cross(V1[i], V2[j])) > 0) {
                res.emplace_back(res.back() + V2[j++]);
            }
            if (j < m - 1 && (V1[i] || V2[j])) {
                res.emplace_back(res.back() + V1[i] + V2[j++]);
            } else {
                res.emplace_back(res.back() + V1[i]);
            }
        }
        for (; j < m - 1; ++j) {
            res.emplace_back(res.back() + V2[j]);
        }
        return res;
    };//�ϲ�͹�� 
    return ConvexHull(mincowski(F.U, G.U), mincowski(F.D, G.D));
}
struct DynamicConvexHull {
    std::set <Point> U, D;
    int in(const std::set <Point> &s, const Point &p, int op) const {
        auto it = s.lower_bound(p);
        if (it == s.end()) return 0;
        if (p == *it) return 2;
        if (it == s.begin()) return 0;
        auto pv = prev(it);
        double d = cross(p - *pv, *it - *pv) * op;
        return sign(d) == 0 ? 2 : (sign(d) > 0 ? 1 : 0);
    }
    int in(const Point &p) const {
        int x = in(U, p, 1), y = in(D, p, -1);
        if (x == 2 || y == 2) return 2;
        return x == 1 && y == 1 ? 1 : 0;
    }
    void insert(std::set <Point> &s, const Point &p, int op) {
        if (in(s, p, op)) return;
        auto it = s.lower_bound(p);
        if (it != s.end()) {
            auto nx = next(it);
            while (nx != s.end() && sign(cross((*it) - p, (*nx) - (*it))) * op >= 0) {
                s.erase(it++);
                nx = next(it);
            }
        }
        if (it != s.begin()) {
            --it;
            std::set <Point>::iterator pv;
            while (it != s.begin()) {
                pv = prev(it);
                if (sign(cross((*it) - (*pv), p - (*it))) * op < 0) break;
                s.erase(it--);
            }
        }
        s.insert(p);
    }
    void insert(const Point &p) {
        insert(U, p, 1);
        insert(D, p, -1);
    }
};


Polygon changeConvex(std::vector<Line> vec) {
    int n = vec.size();
    Polygon res(n);
    for (int i = 0; i < n; ++i) {
        res[i] = intersection(vec[i], vec[i == n - 1 ? 0 : i + 1]);
    }
    return res;
}

// keep only the right part

struct DynamicHalfCut {
    struct info {
        Line v, mx;
        int rk;
        info *lc, *rc;
        info() { clear(); }
        info(const Line &x) : v(x), mx(x) { lc = rc = nullptr; }
        void clear() {
            v = mx = Line();
            lc = rc = nullptr;
        }
        void pushup() {
            if (rc) mx = rc -> mx;
            else mx = v;
        }
    } *rt;
    DynamicHalfCut() {
        rt = nullptr;
    }
    void split(info *&o, const Line &v, info *&x, info *&y, bool flag = false) {
        if (!o) return x = y = nullptr, void();
        if (flag ? !(v < o -> v) : o -> v < v) {
            x = o;
            split(o -> rc, v, x -> rc, y, flag);
            x -> pushup();
        } else {
            y = o;
            split(o -> lc, v, x, y -> lc, flag);
            y -> pushup();
        }
    }
    info* merge(info *x, info *y) {
        if (!x || !y) return x ? x : y;
        if (x -> rk < y -> rk) {
            x -> rc = merge(x -> rc, y);
            x -> pushup();
            return x;
        } else {
            y -> lc = merge(x, y -> lc);
            y -> pushup();
            return y;
        }
    }
    
    info* getFirst(info *&o) {
        if (!o) return nullptr;
        info *x = o, *y;
        while (x -> lc) x = x -> lc;
        split(o, x -> v, x, y, true);
        o = y;
        return x;
    }
    info* getLast(info *&o) {
        if (!o) return nullptr;
        info *x, *y = o;
        while (y -> rc) y = y -> rc;
        split(o, y -> v, x, y);
        o = x;
        return y;
    }
    void insert(const Line &v) {
        info *x, *z;
        split(rt, v, x, z);
        auto pre = [&]() {
            if (x != nullptr) return std::make_pair(getLast(x), 0);
            return std::make_pair(getLast(z), 1);
        };
        auto push_pre = [&](const std::pair<info*, int> &o) {
            if (o.second) {
                z = merge(z, o.first);
            } else {
                x = merge(x, o.first);
            }
        };
        auto nxt = [&]() {
            if (z != nullptr) return std::make_pair(getFirst(z), 1);
            return std::make_pair(getFirst(x), 0);
        };
        auto push_nxt = [&](const std::pair<info*, int> &o) {
            if (o.second) {
                z = merge(o.first, z);
            } else {
                x = merge(o.first, x);
            }
        };
        auto check = [&](const Line &l, const Point &p) {
            return l.left(p) || cap(p, l);
        };
        auto judge = [&](const Line &a, const Line &b, const Line &c) {
            return check(a, intersection(b, c)) && check(c, intersection(b, a));
        };
        std::pair<info*, int> a, b;
        a = nxt();
        if (a.first && ((a.first -> v) || v) && check(a.first -> v, v.A)) {
            push_nxt(a);
            rt = merge(x, z); 
            return;
        }
        b = pre();
        if (a.first && b.first && judge(a.first -> v, v, b.first -> v)) {
            push_pre(b), push_nxt(a);
            rt = merge(x, z);
            return;
        }
        push_pre(b), push_nxt(a);
        a = pre();
        while (a.first && ((a.first -> v) || v)) {
            delete a.first;
            a = pre();
        }
        
        b = pre();
        while (b.first && judge(v, a.first -> v, b.first -> v)) {
            delete a.first;
            a = b;
            b = pre();
        }
        push_pre(b), push_pre(a);
        
        a = nxt(), b = nxt();
        while (b.first && judge(b.first -> v, a.first -> v, v)) {
            delete a.first;
            a = b;
            b = nxt();
        }
        push_nxt(b), push_nxt(a);
        info *y = new info(v);
        y -> rk = rnd();
        rt = merge(merge(x, y), z);
    }
    
    //check if the point P in all halfcuts
    bool doCheck(info *o, const Point &P, bool op) {
        if (!o) return true;
        if ((o -> v).left(P)) return false;
        if (!(o -> lc)) return doCheck(o -> rc, P, op);
        if (!(o -> rc)) return doCheck(o -> lc, P, op);
        Point t = intersection(o -> v, o -> lc -> mx);
        if ((t.y < P.y) ^ op) {
            return doCheck(o -> lc, P, op);
        } else {
            return doCheck(o -> rc, P, op);
        }
    }
    bool check(const Point &P) {
        info *x, *y;
        split(rt, Line(Point(0, 0), Point(1, 0)), x, y);
        bool flag = doCheck(x, P, 1) & doCheck(y, P, 0);
        rt = merge(x, y);
        return flag;
    }
    void dfs(info *o, std::vector<Line> &vec) {
        if (!o) return;
        dfs(o -> lc, vec);
        vec.push_back(o -> v);
        dfs(o -> rc, vec);
    }
    auto getAll() {
        std::vector<Line> vec;
        dfs(rt, vec);
        return vec;
    }
};


#define pw(x) ((x) * (x))
struct Circle {
    Point O; double R;
    Circle() { R = 0; }
    Circle(Point O, double R = 0) : O(O), R(R) {}
    Circle(Point A, Point B) { construct(A, B); }
    Circle(Point A, Point B, Point C) { construct(A, B, C); }
    void construct(Point A, double _ = 0) { O = A; R = _; }
    void construct(Point A, Point B) {
        O = (A + B) * 0.5;
        R = (A - O).mod();
    }
    void construct(Point A, Point B, Point C) {
        assert(sign(cross(A - B, A - C)) != 0);
        Point P1 = (A + B) * 0.5, P2 = (A + C) * 0.5;
        O = intersection(Line(P1, P1 + (B - A).normal()), Line(P2, P2 + (C - A).normal()));
        R = (A - O).mod();
    }
};
double area(const Circle &O) { return pi * O.R * O.R; }
double area(double theta, double R) { return 0.5 * theta * R; }
int in(const Point &p, const Circle &C) {
    return sign(C.R - dis(p, C.O)) > 0 ? 1 : (sign(C.R - dis(p, C.O)) == 0 ? 2 : 0); 
}
int cap(Circle A, Circle B) {
    if (A.R < B.R) std::swap(A, B);
    double d = dis(A.O, B.O);
    if (d < A.R) return sign(d + B.R - A.R) > 0 ? -1 : (sign(d + B.R - A.R) == 0 ? -2 : -3);
    return sign(A.R + B.R - d) > 0 ? 0 : (sign(A.R + B.R - d) == 0 ? 2 : 1);
}
int cap(Point p, Circle C) {
    double d = sign(dis(p, C.O) - C.R);
    return d == 0 ? 2 : (d < 0 ? 1 : 0);
}
int cap(Line l, Circle C) {
    double d = sign(dis(C.O, l) - C.R);
    return d == 0 ? 2 : (d < 0 ? 1 : 0);
}
double caplen(Line l, Circle C) {
    double d = dis(C.O, l);
    if (d >= C.R) return 0;
    return sqrt(C.R * C.R - d * d) * 2;
}
Circle getMinCircleCover(Polygon P) { //Circle((0,0),0) when P is empty
    int n = (int)P.size();
    shuffle(P.begin(), P.end(), rnd);
    Circle C(Point(0, 0));
    for (int i = 0; i < n; ++i) if (!in(P[i], C)) {
        C.construct(P[i]);
        for (int j = 0; j < i; ++j) if (!in(P[j], C)) {
            C.construct(P[i], P[j]);
            for (int k = 0; k < j; ++k) if (!in(P[k], C)) {
                C.construct(P[i], P[j], P[k]);
            }
        }
    }
    return C;
}

// size(res) = 0 means they not 'cap'
// size(res) = 1 means there is exactly one point
// size(res) = 2 means 2, and from res[0] to res[1] clockwise is covered by B
std::vector<Point> intersection(Circle A, Circle B) {
    std::vector<Point> res;
    int type = cap(A, B);
    if (type == -3 || type == 0) return res;
    if (type == 2) {
        res.emplace_back((B.O - A.O) * A.R / (A.R + B.R) + A.O);
        return res;
    }
    if (type == -2) {
        res.emplace_back((B.O - A.O) * B.R / (A.R - B.R) + B.O);
        return res;
    }
    double dis = (B.O - A.O).mod();
    double theta = acos((dis * dis + A.R * A.R - B.R * B.R) / (2.0 * A.R * dis));
    res.emplace_back((B.O - A.O).rotate(theta) + A.O);
    res.emplace_back((B.O - A.O).rotate(-theta) + A.O);
    return res;
}
//δ��֤ 
}

int main() {
    
    return 0;
}