Geometry

#include <bits/stdc++.h>
typedef long long ll;
typedef long double ld;
using namespace std;

namespace geoma {
#define Doub ld
#define Long ll
#define Fix cout << fixed << setprecision(10)
    const long double Pi = acos(-1);
    const long double Eps = 1e-8;
    const Long magic = -2423423324;
    bool is_equal(Doub a, Doub b) {
        return fabs(a - b) < Eps;
    }
    struct line {
        Doub A, B, C;
        line() {
            A = B = C = 0;
        }
        line(Doub _A, Doub _B, Doub _C) {
            A = _A, B = _B, C = _C;
        }
        Doub get_y(Doub &x) {
            if (B == 0) {
                return 0;
            }
            return (-C - A * x) / B;
        }
        Doub get_x(Doub &y) {
            if (A == 0) {
                return 0;
            }
            return (-C - B * y) / A;
        }
        void normalize() {
            auto opr = hypot(A, B);
            A /= opr, B /= opr, C /= opr;
        }
        bool operator == (const line &other) {
            return A == other.A && B == other.B && C == other.C;
        }
    };
    istream& operator >> (istream& in, line &t) { return in >> t.A >> t.B >> t.C; }
    ostream& operator << (ostream& out, line t) { Fix; return out << t.A << ' ' << t.B << ' ' << t.C; }
    struct point {
        Doub x, y;
        point() {
            x = y = 0;
        }
        point(Doub _x, Doub _y) {
            x = _x, y = _y;
        }
        Long operator * (point other) {
            return x * other.y - other.x * y;
        }
        Long operator ^ (point other) {
            return x * other.x + y * other.y;
        }
        point operator- (point other) {
            return {x - other.x, y - other.y};
        }
        point operator+ (point other) {
            return {x + other.x, y + other.y};
        }
        point operator*(Long k) {
            return {x * k, y * k};
        }
        point operator*(Doub k) {
            return {x * k, y * k};
        }
        void operator+= (point b) {
            x += b.x;
            y += b.y;
        }
        void operator-= (point b) {
            x -= b.x;
            y -= b.y;
        }
        Doub len() {
            return sqrt(1.0*(x * x + y * y));
        }
        Long klen() {
            return x * x + y * y;
        }
        bool operator == (const point &other) {
            return x == other.x && y == other.y;
        }
        Doub dist_to_line(line t) {
            t.normalize();
            return t.A * x + t.B * y + t.C;
        }
        Doub ang1(const point &b) { // [0..pi]
            point a = (*this);
            Doub num = a * b;
            Doub den = a ^ b;
            if (num == 0 && den > 0)
                return 0;
            ld ret = fabs(atan2(num, den));
            if (ret == 0) {
                ret = Pi;
            }
            return ret;
        }
        Doub ang2(const point &b) { // [-pi..pi]
            point a = (*this);
            Doub num = a * b;
            Doub den = a ^ b;
            return atan2(num, den);
        }
        bool is_lie(line t) {
            return is_equal(t.A * x + t.B * y + t.C, 0);
        }
        point projection(line t) {
            t.normalize();
            auto len = dist_to_line(t);
            point ret = (*this);
            point ve(t.A, t.B);
            ret -= ve * len;
            bool inv = is_lie(t);
            //assert(inv);
            return ret;
        }
        void normalize() {
            auto op = hypot(x, y);
            x /= op, y /= op;
        }
    };
    istream& operator >> (istream& in, point &t) { return in >> t.x >> t.y; }
    ostream& operator << (ostream& out, point t) { Fix; return out << t.x << ' ' << t.y; }
    line toline(const point &a, const point &b) {
        line ret;
        ret.A = a.y - b.y; ret.B = b.x - a.x; ret.C = a.x * b.y - b.x * a.y;
        return ret;
    }
    struct circle {
        point C;
        Doub R;
        auto cross_with_line(line x) {
            auto t = (*this);
            auto d = fabs(t.C.dist_to_line(x));
            vector <point> ans;
            if (d > t.R)
                return ans;
            x.normalize();
            point ve(x.A, x.B);
            Doub len = ve.len();
            point var_1 = t.C; var_1 += ve * d;
            point var_2 = t.C; var_2 -= ve * d;
            auto len_var_1 = fabs(var_1.dist_to_line(x));
            auto len_var_2 = fabs(var_2.dist_to_line(x));
            if (len_var_1 > len_var_2)
                swap(var_1, var_2);
            if (is_equal(d, t.R)) {
                x.normalize();
                ans.push_back(var_1);
                return ans;
            }
            if (is_equal(d, 0)) {
                x.normalize();
                point ve1 = {-x.B, x.A};
                auto ans1 = t.C;
                ans1 += ve1 * t.R;
                ans.push_back(ans1);
                ans1 = t.C;
                ans1 -= ve1 * t.R;
                ans.push_back(ans1);
                return ans;
            }
            Doub k1 = hypot(t.R, len);
            auto ve1 = var_1 - t.C;
            ve1 = {-ve1.y, ve1.x};
            len = ve1.len() * 2;
            point ans1 = {var_1.x + ve1.x * k1 / len , var_1.y + ve1.y * k1 / len};
            point ans2 = {var_1.x - ve1.x * k1 / len , var_1.y - ve1.y * k1 / len};
            ans.push_back(ans1), ans.push_back(ans2);
            return ans;
        }
        auto cross_with_circle(circle b) {
            circle a = (*this);
            vector <point> ret;
            if (a.C == b.C) {
                if (a.R == b.R) {
                    point clu(-2423423324, -42342423);
                    ret.push_back(clu);
                }
                return ret;
            }
            auto l = (a.C - b.C).len();
            Doub v[3];
            v[0] = a.R;
            v[1] = b.R;
            v[2] = l;
            sort(v, v + 3);
            if (v[0] + v[1] < v[2]) {
                return ret;
            }
            if (v[0] + v[1] == v[2]) {
                point ans;
                ans = a.C;
                auto ve = b.C - a.C;
                ve.normalize();
                ans = a.C;
                ans.x += ve.x * (a.R + b.R);
                ans.y += ve.y * (a.R + b.R);
                if (ans == b.C) {
                    ans.x -= ve.x * b.R;
                    ans.y -= ve.y * b.R;
                    ret.push_back(ans);
                    return ret;
                }
                else {
                    ans = a.C;
                    ans.x += ve.x * (a.R - b.R);
                    ans.y += ve.y  * (a.R - b.R);
                    if (ans == b.C) {
                        ans = a.C;
                        ans.x += ve.x * a.R;
                        ans.y += ve.y * a.R;
                        ret.push_back(ans);
                        return ret;
                    }
                    else {
                        swap(a, b);
                        ans = a.C;
                        ve = b.C - a.C;
                        ve.normalize();
                        ans.x += ve.x * a.R;
                        ans.y += ve.y * a.R;
                        ret.push_back(ans);
                        return ret;
                    }
                }
            }
            else {
                Doub r2 = b.R;
                Doub r1 = a.R;
                Doub cos = (l*l + r2*r2 - r1*r1) / (2 * l * r2);
                Doub O2H = r2 * cos;
                Doub PH = sqrt(r2 * r2 - O2H * O2H);
                point H = b.C;
                auto ve = a.C - b.C;
                ve.normalize();
                H.x = H.x + ve.x * O2H;
                H.y = H.y + ve.y * O2H;
                auto ln = toline(b.C, a.C);
                ln.normalize();
                point ans1 = H;
                ans1.x = ans1.x + ln.A * PH;
                ans1.y = ans1.y + ln.B * PH;
                point ans2 = H;
                ans2.x = ans2.x - ln.A * PH;
                ans2.y = ans2.y - ln.B * PH;
                ret.push_back(ans1), ret.push_back(ans2);
                return ret;
            }
        }
    };
    istream &operator >> (istream &in, circle &t) { return in >> t.C >> t.R; }
    point inter(line fi, line se) {
        Doub y = (se.A * fi.C - fi.A * se.C) / (fi.A * se.B - se.A * fi.B);
        Doub x;
        if (fi.A != 0)
            x = fi.get_x(y);
        else
            x = se.get_x(y);
        return {x, y};
    }
    struct segment {
        point begin, end;
        bool is_beetwen(point a) {
            if (a == begin || a == end)
                return 1;
            return (((end - begin) * (a - begin) == 0) && (((begin - a) ^ (end - a)) <= 0));
        }
        bool is_cross(segment b) {
            auto a1 = begin, a2 = end, b1 = b.begin, b2 = b.end;
            if ((a2 - a1) * (b2 - b1) == 0) {
                ll sc = (a2 - a1) * (b2 - a1);
                if (sc != 0) return 0;
                return (b.is_beetwen(a1) || b.is_beetwen(a2) || segment {a1, a2}.is_beetwen(b1) || segment {a1, a2}.is_beetwen(b2));
            }
            else {
                return ((((a2 - a1) * (b2 - a1)) * ((a2 - a1) * (b1 - a1)) <= 0) && (((b1 - b2) * (a1 - b2)) * ((b1 - b2) * (a2 - b2)) <= 0));
            }
        }
        Doub dist_to_point(point s) {
            point get = s.projection(toline(begin, end));
            segment c;
            c.begin = begin, c.end = end;
            if (c.is_beetwen(get)) {
                return s.dist_to_line(toline(begin, end));
            }
            else {
                return min((s - begin).len(), (s - end).len());
            }
        }
        Doub len() {
            return (end - begin).len();
        }
        Doub dist_to_segment(segment b) {
            segment a = (*this);
            Doub ans;
            if (a.len() == 1 && b.len() == 1) {
                ans = (a.begin - b.begin).len();
            }
            else if (a.len() == 1) {
                ans = fabs(b.dist_to_point(a.begin));
            }
            else if (b.len() == 1) {
                ans = fabs(a.dist_to_point(b.begin));
            }
            else {
                if (a.is_cross(b)) {
                    ans = 0;
                }
                else {
                    ans = fabs(a.dist_to_point(b.begin));
                    ans = min(ans, fabs(a.dist_to_point(b.end)));
                    ans = min(ans, fabs(b.dist_to_point(a.begin)));
                    ans = min(ans, fabs(b.dist_to_point(a.end)));
                }
            }
            return ans;
        }
        bool is_point_on_ray(point a) {
            auto b = begin;
            auto c = end;
            return ((((a - b) ^ (c - b)) >= 0 && ((a - b)*(c - b)) == 0));
        }
        Doub dist_to_ray(segment b) {
            auto a = (*this);
            if (b.is_point_on_ray(a.begin) || b.is_point_on_ray(a.end)|| a.is_point_on_ray(b.begin) || a.is_point_on_ray(b.end)) {
                return 0;
            }

            ld ans = a.dist_to_segment(b);

            line fi = toline(a.begin, a.end);
            line se = toline(b.begin, b.end);
            auto teta = inter(fi, se);

            bool parallel = (((a.end - a.begin) ^ (b.end - b.begin)) == 0);

            if (!parallel && b.is_point_on_ray(teta)&& a.is_point_on_ray(teta)) {
                ans = 0;
            }

            // работаем с a
            auto T = a.begin.projection(se);
            if ( b.is_point_on_ray(T))
                ans = min(ans, fabs(a.begin.dist_to_line(se)));
            T = a.end.projection(se);
            if (b.is_point_on_ray(T))
                ans = min(ans, fabs(a.end.dist_to_line(se)));

            // работаем с b
            T = b.begin.projection(fi);
            if (a.is_point_on_ray(T))
                ans = min(ans, fabs(b.begin.dist_to_line(fi)));
            T = b.end.projection(fi);
            if (a.is_point_on_ray(T))
                ans = min(ans, fabs(b.end.dist_to_line(fi)));
            return ans;
        }
    };
    istream& operator >> (istream& in, segment &t) { return in >> t.begin >> t.end; }
    int sign(int a) {
        if (a == 0)
            return 0;
        if (a < 0)
            return -1;
        return 1;
    }
    bool in_ang(point a, point o, point b, point c) { // aob - угол
        if (segment {c, o}.is_beetwen(b)) return 1;
        if (segment {c, o}.is_beetwen(a)) return 1;
        ll t1 = (b - o) * (c - o);
        ll t2 = (c - o) * (a - o);
        ll t3 = (b - o) * (a - o);
        return sign(t1) == sign(t2) && sign(t2) == sign(t3);
    }
    struct polygon {
        vector <point> a;
        Doub square() {
            Doub res = 0;
            for (unsigned i = 0; i < a.size(); i++) {
                point
                    p1 = i ? a[i - 1] : a.back(),
                    p2 = a[i];
                res += (p1.x - p2.x) * (p1.y + p2.y);
            }
            return fabs(res) / 2;
        }
        bool in_polygon(point c) {
            int n = a.size();
            auto z = a[0];
            int l = 1, r = n - 2;
            if (!in_ang(a[l], z, a[n - 1], c)) {
                return 0;
            }
            while (r - l > 1) {
                int mid = (l + r) / 2;
                if (in_ang(a[mid], z, a[n - 1], c))
                    l = mid;
                else
                    r = mid;
            }
            r = l + 1;
            return (in_ang(a[l], z, a[r], c) && in_ang(a[l], a[r], z, c));
        }
    };
    bool cmph(point a, point b) {
        if (a.y != b.y)
            return a.y < b.y;
        return a.x < b.x;
    }
    bool ch_cmp(point a, point b) {
        if ((a * b) == 0) {
            return a.klen() < b.klen();
        }
        return (a * b) > 0;
    }
    auto convex_hull(polygon a) {
        point po(1e18, 1e18);
        point fi = *min_element(a.a.begin(), a.a.end(), cmph);
        for (auto &i : a.a)
            i -= fi;
        sort(a.a.begin(), a.a.end(), ch_cmp);
        vector <point> hull;
        hull.push_back(a.a[0]);
        for (int i = 1; i < a.a.size(); ++i) {
            while (hull.size() > 1 && ((a.a[i] - hull.back()) * (hull[hull.size() - 2] - hull.back())) <= 0) {
                hull.pop_back();
            }
            hull.push_back(a.a[i]);
        }
        for (auto &i : hull)
            i += po;
        return hull;
    }
    bool in_circle(point p, circle c) {
        return (c.C - p).len() < c.R;
    }
    auto tangent(circle a, point b) {
        ld AO = (a.C - b).len();
        vector <point> ret;
        if (AO == a.R) {
            ret.push_back({b.x, b.y});
            return ret;
        }
        if (AO < a.R) {
            return ret;
        }
        ld OD = a.R;
        ld AD = sqrt(AO * AO - OD * OD);
        ld CD = AD * OD / AO;
        ld OC = sqrt(OD * OD - CD * CD);
        point vec = b - a.C;
        auto cur = a.C;
        cur.x += (vec.x * OC / AO);
        cur.y += (vec.y * OC / AO);
        point vec2 = {-vec.y, vec.x};
        auto p1 = cur;
        p1.x += vec2.x * CD / AO;
        p1.y += vec2.y * CD / AO;
        auto p2 = cur;
        p2.x -= vec2.x * CD / AO;
        p2.y -= vec2.y * CD / AO;
        ret.push_back(p1);
        ret.push_back(p2);
        return ret;
    }
    line bisector(point x, point y, point z) {
        if (x.x == y.x && x.x == z.x) {
            auto ans = toline(x, {x.x + 1, x.y});
            return ans;
        }
        if (x.y == y.y && x.y == z.y) {
            auto ans = toline(x, {x.x, x.y + 1});
            cout << fixed << setprecision(10) << ans.A << ' ' << ans.B << ' ' << ans.C << '\n';
            exit(0);
        }
        auto v1 = y - x;
        v1.normalize();
        auto v2 = z - x;
        v2.normalize();
        auto t = x;
        t += v1;
        auto g = x;
        g += v2;
        auto ans = toline(x, {(g.x + t.x) / 2, (g.y + t.y) / 2});
        return ans;
    }
};

using namespace geoma;