NicheLibrary

This documentation is automatically generated by NotLeonian/competitive-verifier (forked from competitive-verifier/competitive-verifier)

View the Project on GitHub NotLeonian/NicheLibrary

:heavy_check_mark: verify/aoj-cgl-4-c.test.cpp

Depends on

Code

// competitive-verifier: PROBLEM https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_4_C
// competitive-verifier: ERROR 0.00001

#include <cassert>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <vector>

#include "../geometry/line-convex-polygon-intersection.hpp"

struct Point {
    long long x;
    long long y;

    Point() : x(0), y(0) {}
    Point(long long x_, long long y_) : x(x_), y(y_) {}

    Point operator+(const Point &other) const {
        return Point(x + other.x, y + other.y);
    }
    Point operator-(const Point &other) const {
        return Point(x - other.x, y - other.y);
    }
};

struct RealPoint {
    long double x;
    long double y;

    RealPoint() : x(0), y(0) {}
    RealPoint(long double x_, long double y_) : x(x_), y(y_) {}

    RealPoint operator+(const RealPoint &other) const {
        return RealPoint(x + other.x, y + other.y);
    }
    RealPoint operator-(const RealPoint &other) const {
        return RealPoint(x - other.x, y - other.y);
    }
    RealPoint operator*(long double k) const { return RealPoint(x * k, y * k); }
};

__int128_t cross_value(const Point &a, const Point &b) {
    return static_cast<__int128_t>(a.x) * static_cast<__int128_t>(b.y) -
           static_cast<__int128_t>(a.y) * static_cast<__int128_t>(b.x);
}

int sign_value(__int128_t x) {
    if (x < 0) {
        return -1;
    }
    if (x > 0) {
        return 1;
    }
    return 0;
}

Point read_point() {
    long long x, y;
    std::cin >> x >> y;
    return Point(x, y);
}

RealPoint to_real_point(const Point &p) {
    return RealPoint(static_cast<long double>(p.x),
                     static_cast<long double>(p.y));
}

long double cross_ld(const RealPoint &a, const RealPoint &b) {
    return a.x * b.y - a.y * b.x;
}

long double dot_ld(const RealPoint &a, const RealPoint &b) {
    return a.x * b.x + a.y * b.y;
}

bool same_point(const RealPoint &a, const RealPoint &b) {
    return std::abs(a.x - b.x) <= 1e-12L && std::abs(a.y - b.y) <= 1e-12L;
}

void add_unique(std::vector<RealPoint> &points, const RealPoint &point) {
    for (const auto &q : points) {
        if (same_point(q, point)) {
            return;
        }
    }
    points.push_back(point);
}

RealPoint line_intersection(const Point &line_a, const Point &line_b,
                            const Point &segment_a, const Point &segment_b) {
    const RealPoint u = to_real_point(segment_a);
    const RealPoint v = to_real_point(segment_b);
    const RealPoint dir = to_real_point(line_b - line_a);
    const RealPoint edge = v - u;
    const long double denominator = cross_ld(edge, dir);
    const long double numerator =
        cross_ld(to_real_point(line_a - segment_a), dir);
    const long double t = numerator / denominator;
    return u + edge * t;
}

std::vector<RealPoint> brute_intersections(const std::vector<Point> &polygon,
                                           const Point &line_a,
                                           const Point &line_b) {
    const int n = static_cast<int>(polygon.size());
    const Point direction = line_b - line_a;
    std::vector<RealPoint> points;
    for (int i = 0; i < n; ++i) {
        const Point &u = polygon[i];
        const Point &v = polygon[(i + 1) % n];
        const __int128_t hu = cross_value(u - line_a, direction);
        const __int128_t hv = cross_value(v - line_a, direction);
        const int su = sign_value(hu);
        const int sv = sign_value(hv);
        if (su == 0 && sv == 0) {
            add_unique(points, to_real_point(u));
            add_unique(points, to_real_point(v));
        } else if (su == 0) {
            add_unique(points, to_real_point(u));
        } else if (sv == 0) {
            add_unique(points, to_real_point(v));
        } else if (su != sv) {
            add_unique(points, line_intersection(line_a, line_b, u, v));
        }
    }
    if (points.size() <= 2) {
        return points;
    }
    const RealPoint dir = to_real_point(direction);
    int min_index = 0;
    int max_index = 0;
    for (int i = 1; i < static_cast<int>(points.size()); ++i) {
        if (dot_ld(points[i], dir) < dot_ld(points[min_index], dir)) {
            min_index = i;
        }
        if (dot_ld(points[i], dir) > dot_ld(points[max_index], dir)) {
            max_index = i;
        }
    }
    return {points[min_index], points[max_index]};
}

std::vector<RealPoint> convex_cut(const std::vector<Point> &polygon,
                                  const Point &line_a, const Point &line_b) {
    const int n = static_cast<int>(polygon.size());
    const Point direction = line_b - line_a;
    std::vector<RealPoint> result;
    for (int i = 0; i < n; ++i) {
        const Point &now = polygon[i];
        const Point &next = polygon[(i + 1) % n];

        // AOJ CGL_4_C は有向直線 line_a -> line_b の左側を残す。
        const __int128_t current = cross_value(direction, now - line_a);
        const __int128_t next_value = cross_value(direction, next - line_a);

        if (sign_value(current) >= 0) {
            result.push_back(to_real_point(now));
        }
        if ((current < 0 && next_value > 0) ||
            (current > 0 && next_value < 0)) {
            result.push_back(line_intersection(line_a, line_b, now, next));
        }
    }
    return result;
}

long double polygon_area(const std::vector<RealPoint> &polygon) {
    const int n = static_cast<int>(polygon.size());
    long double area = 0;
    for (int i = 0; i < n; ++i) {
        area += cross_ld(polygon[i], polygon[(i + 1) % n]);
    }
    return std::abs(area) / 2;
}

std::vector<Point>
remove_collinear_vertices(const std::vector<Point> &polygon) {
    const int n = static_cast<int>(polygon.size());
    std::vector<Point> result;
    result.reserve(n);
    for (int i = 0; i < n; ++i) {
        const Point &prev = polygon[(i + n - 1) % n];
        const Point &cur = polygon[i];
        const Point &next = polygon[(i + 1) % n];
        if (sign_value(cross_value(cur - prev, next - cur)) == 0) {
            continue;
        }
        result.push_back(cur);
    }
    assert(result.size() >= 3);
    return result;
}

void verify_intersection(const std::vector<Point> &polygon, const Point &line_a,
                         const Point &line_b) {
    const auto fast = line_polygon_intersection(polygon, line_a, line_b);
    const auto slow = brute_intersections(polygon, line_a, line_b);
    assert(fast.size() == slow.size());
    if (fast.empty()) {
        return;
    }

    auto to_real = [](const auto &point) {
        return RealPoint(point.template x_as<long double>(),
                         point.template y_as<long double>());
    };

    if (fast.size() == 1) {
        assert(same_point(to_real(fast[0]), slow[0]));
        return;
    }
    const RealPoint p0 = to_real(fast[0]);
    const RealPoint p1 = to_real(fast[1]);
    const bool ok = (same_point(p0, slow[0]) && same_point(p1, slow[1])) ||
                    (same_point(p0, slow[1]) && same_point(p1, slow[0]));
    assert(ok);
}

int main() {
    int n;
    std::cin >> n;
    std::vector<Point> polygon(n);
    for (auto &p : polygon) {
        p = read_point();
    }
    polygon = remove_collinear_vertices(polygon);

    int q;
    std::cin >> q;
    std::cout << std::fixed << std::setprecision(20);
    while (q--) {
        const Point line_a = read_point();
        const Point line_b = read_point();
        verify_intersection(polygon, line_a, line_b);
        std::cout << polygon_area(convex_cut(polygon, line_a, line_b)) << '\n';
    }
}

#line 1 "verify/aoj-cgl-4-c.test.cpp"
// competitive-verifier: PROBLEM https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_4_C
// competitive-verifier: ERROR 0.00001

#include <cassert>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <vector>

#line 1 "geometry/line-convex-polygon-intersection.hpp"



// 反時計回りの面積正の狭義凸多角形と直線の共通部分を求める。
// 共通部分が空集合なら 0 個、1 点なら 1 個、線分ならその端点 2 個を返す。
// 前提: polygon は面積正の狭義凸多角形で、反時計回りで、
//       3 頂点連続で一直線にならず、line_a != line_b。
// 注意: convex_hull の結果が 1 点または 2 点になる退化ケースは前提外。
//       呼び出し側で点または線分として処理すること。
// 座標型が整数なら内部は整数演算で処理し、交点を有理表現で返す。
// 主な対象は 10^9 級制約であり、非常に大きい整数座標、特に 10^18 級では
// 返り値生成時にオーバーフローする可能性がある。
// 計算量 O(log N)。

#line 16 "geometry/line-convex-polygon-intersection.hpp"
#include <complex>
#include <type_traits>
#include <utility>
#line 20 "geometry/line-convex-polygon-intersection.hpp"

namespace line_convex_polygon_intersection_internal {
template <class T> struct is_integral : std::is_integral<T> {};
template <class T> struct is_signed : std::is_signed<T> {};
#ifdef __SIZEOF_INT128__
template <> struct is_integral<__int128_t> : std::true_type {};
template <> struct is_integral<__uint128_t> : std::true_type {};
template <> struct is_signed<__int128_t> : std::true_type {};
template <> struct is_signed<__uint128_t> : std::false_type {};
#endif

template <class T>
using enable_if_integral_t = std::enable_if_t<is_integral<T>::value, int>;

template <class Point, class = void> struct has_member_xy : std::false_type {};
template <class Point>
struct has_member_xy<Point,
                     std::void_t<decltype(std::declval<const Point &>().x),
                                 decltype(std::declval<const Point &>().y)>>
    : std::true_type {};

template <class Point> decltype(auto) get_x(const Point &point) {
    if constexpr (has_member_xy<Point>::value) {
        return point.x;
    } else {
        using std::real;
        return real(point);
    }
}

template <class Point> decltype(auto) get_y(const Point &point) {
    if constexpr (has_member_xy<Point>::value) {
        return point.y;
    } else {
        using std::imag;
        return imag(point);
    }
}

template <class Point>
using coordinate_t =
    std::remove_cvref_t<decltype(get_x(std::declval<const Point &>()))>;

template <class T, class = void> struct wide_integer {};
template <class T>
struct wide_integer<
    T, std::enable_if_t<is_integral<T>::value && is_signed<T>::value>> {
    using type = __int128_t;
};
template <class T>
struct wide_integer<
    T, std::enable_if_t<is_integral<T>::value && !is_signed<T>::value>> {
    using type = __uint128_t;
};

template <class T> using wide_integer_t = typename wide_integer<T>::type;

template <class Point, bool = is_integral<coordinate_t<Point>>::value>
struct calc_type {
    using type = long double;
};
template <class Point> struct calc_type<Point, true> {
    using type = wide_integer_t<coordinate_t<Point>>;
};

template <class Point> using calc_t = typename calc_type<Point>::type;

template <class T> int sign_value(const T &value) {
    if constexpr (is_integral<T>::value) {
        if (value < 0) {
            return -1;
        }
        if (value > 0) {
            return 1;
        }
        return 0;
    } else {
        constexpr long double eps = 1e-12L;
        if (value < -eps) {
            return -1;
        }
        if (value > eps) {
            return 1;
        }
        return 0;
    }
}

template <class T> bool equals_value(const T &lhs, const T &rhs) {
    return sign_value(lhs - rhs) == 0;
}

template <class Point> calc_t<Point> to_calc_x(const Point &point) {
    return static_cast<calc_t<Point>>(get_x(point));
}

template <class Point> calc_t<Point> to_calc_y(const Point &point) {
    return static_cast<calc_t<Point>>(get_y(point));
}

template <class Point> calc_t<Point> cross(const Point &lhs, const Point &rhs) {
    return to_calc_x(lhs) * to_calc_y(rhs) - to_calc_y(lhs) * to_calc_x(rhs);
}

template <class Point>
Point make_point(const calc_t<Point> &x, const calc_t<Point> &y) {
    using Coord = coordinate_t<Point>;
    return Point(static_cast<Coord>(x), static_cast<Coord>(y));
}

template <class T, enable_if_integral_t<T> = 0> T abs_value(T value) {
    if constexpr (is_signed<T>::value) {
        return value < 0 ? -value : value;
    } else {
        return value;
    }
}

template <class T, enable_if_integral_t<T> = 0> T gcd_value(T a, T b) {
    a = abs_value(a);
    b = abs_value(b);
    while (b != 0) {
        T t = a % b;
        a = b;
        b = t;
    }
    return a;
}

template <class T, enable_if_integral_t<T> = 0> T gcd3_value(T a, T b, T c) {
    return gcd_value(gcd_value(a, b), c);
}

template <class Func> int first_non_negative(int length, Func &&func) {
    int left = 0;
    int right = length;
    while (left + 1 < right) {
        int mid = (left + right) >> 1;
        if (sign_value(func(mid)) < 0) {
            left = mid;
        } else {
            right = mid;
        }
    }
    return right;
}

template <class Func> int last_non_positive(int length, Func &&func) {
    int left = 0;
    int right = length;
    while (left < right) {
        int mid = (left + right + 1) >> 1;
        if (sign_value(func(mid)) <= 0) {
            left = mid;
        } else {
            right = mid - 1;
        }
    }
    return left;
}

template <class Func> int last_non_negative(int length, Func &&func) {
    int left = 0;
    int right = length;
    while (left < right) {
        int mid = (left + right + 1) >> 1;
        if (sign_value(func(mid)) >= 0) {
            left = mid;
        } else {
            right = mid - 1;
        }
    }
    return left;
}

inline int positive_mod(int index, int mod) {
    int ret = index % mod;
    if (ret < 0) {
        ret += mod;
    }
    return ret;
}

inline int distance_forward(int from, int to, int n) {
    if (from <= to) {
        return to - from;
    }
    return to + n - from;
}

template <class Point, class Func>
std::pair<int, int> find_extreme_vertices(int n, Func &&height) {
    assert(n >= 3);

    int minimum_index = -1;
    int maximum_index = -1;
    const auto first = height(0);
    const auto last = height(n - 1);
    if (equals_value(first, last)) {
        const auto second = height(1);
        if (sign_value(first - second) < 0) {
            minimum_index = n - 1;
        } else {
            maximum_index = n - 1;
        }
    } else {
        const auto second = height(1);
        if (sign_value(last - first) > 0 && sign_value(first - second) <= 0) {
            minimum_index = 0;
        } else if (sign_value(last - first) < 0 &&
                   sign_value(first - second) >= 0) {
            maximum_index = 0;
        } else {
            const auto denom = static_cast<calc_t<Point>>(n - 1);
            auto g_scaled = [&](int index) {
                const auto base =
                    first * static_cast<calc_t<Point>>(n - 1 - index) +
                    last * static_cast<calc_t<Point>>(index);
                const auto value = height(index) * denom;
                return sign_value(base - value) < 0 ? base : value;
            };
            minimum_index = first_non_negative(n - 1, [&](int index) {
                return g_scaled(index + 1) - g_scaled(index);
            });
        }
    }

    if (maximum_index == -1) {
        const int base = minimum_index - n;
        const int k = first_non_negative(n, [&](int index) {
            return height(base + index) - height(base + index + 1);
        });
        maximum_index = positive_mod(minimum_index + k, n);
    } else {
        const int base = maximum_index - n;
        const int k = first_non_negative(n, [&](int index) {
            return height(base + index + 1) - height(base + index);
        });
        minimum_index = positive_mod(maximum_index + k, n);
    }
    return {minimum_index, maximum_index};
}
} // namespace line_convex_polygon_intersection_internal

template <class T> struct LinePolygonIntersectionPoint {
    T x_numerator;
    T y_numerator;
    T denominator;

    template <class Real> Real x_as() const {
        return static_cast<Real>(x_numerator) / static_cast<Real>(denominator);
    }

    template <class Real> Real y_as() const {
        return static_cast<Real>(y_numerator) / static_cast<Real>(denominator);
    }

    template <class Point> Point to_point() const {
        using Coord =
            line_convex_polygon_intersection_internal::coordinate_t<Point>;
        if constexpr (line_convex_polygon_intersection_internal::is_integral<
                          Coord>::value) {
            assert(denominator == 1);
            return Point(static_cast<Coord>(x_numerator),
                         static_cast<Coord>(y_numerator));
        } else {
            return Point(static_cast<Coord>(x_as<long double>()),
                         static_cast<Coord>(y_as<long double>()));
        }
    }
};

namespace line_convex_polygon_intersection_internal {
template <class Point, bool = is_integral<coordinate_t<Point>>::value>
struct result_value_type {
    using type = Point;
};
template <class Point> struct result_value_type<Point, true> {
    using type = LinePolygonIntersectionPoint<calc_t<Point>>;
};

template <class Point>
using result_value_t = typename result_value_type<Point>::type;
} // namespace line_convex_polygon_intersection_internal

template <class Point>
using LinePolygonIntersectionValue =
    line_convex_polygon_intersection_internal::result_value_t<Point>;

template <class Point>
using LinePolygonIntersectionResult =
    std::vector<LinePolygonIntersectionValue<Point>>;

namespace line_convex_polygon_intersection_internal {
template <class Point>
LinePolygonIntersectionPoint<calc_t<Point>>
make_integral_point(const Point &point) {
    using Calc = calc_t<Point>;
    return {static_cast<Calc>(get_x(point)), static_cast<Calc>(get_y(point)),
            Calc(1)};
}

template <class Point>
LinePolygonIntersectionPoint<calc_t<Point>>
line_edge_intersection_integral(const Point &line_a, const Point &line_b,
                                const Point &segment_a,
                                const Point &segment_b) {
    using Calc = calc_t<Point>;
    const Point direction = line_b - line_a;
    const Point edge = segment_b - segment_a;
    Calc denominator = cross(edge, direction);
    assert(denominator != 0);
    Calc numerator = cross(line_a - segment_a, direction);

    Calc x_numerator = static_cast<Calc>(get_x(segment_a)) * denominator +
                       static_cast<Calc>(get_x(edge)) * numerator;
    Calc y_numerator = static_cast<Calc>(get_y(segment_a)) * denominator +
                       static_cast<Calc>(get_y(edge)) * numerator;

    if (denominator < 0) {
        denominator = -denominator;
        x_numerator = -x_numerator;
        y_numerator = -y_numerator;
    }

    const Calc g = gcd3_value(x_numerator, y_numerator, denominator);
    if (g != 0) {
        x_numerator /= g;
        y_numerator /= g;
        denominator /= g;
    }
    return {x_numerator, y_numerator, denominator};
}

template <class Point>
Point line_edge_intersection_floating(const Point &line_a, const Point &line_b,
                                      const Point &segment_a,
                                      const Point &segment_b) {
    using Calc = calc_t<Point>;
    const Point direction = line_b - line_a;
    const Point edge = segment_b - segment_a;
    const Calc denominator = cross(edge, direction);
    assert(sign_value(denominator) != 0);
    const Calc numerator = cross(line_a - segment_a, direction);
    const Calc t = numerator / denominator;
    const Calc x = static_cast<Calc>(get_x(segment_a)) +
                   static_cast<Calc>(get_x(edge)) * t;
    const Calc y = static_cast<Calc>(get_y(segment_a)) +
                   static_cast<Calc>(get_y(edge)) * t;
    return make_point<Point>(x, y);
}

template <class Point>
result_value_t<Point> make_vertex_result(const Point &point) {
    if constexpr (is_integral<coordinate_t<Point>>::value) {
        return make_integral_point(point);
    } else {
        return point;
    }
}

template <class Point>
result_value_t<Point> make_edge_result(const Point &line_a, const Point &line_b,
                                       const Point &segment_a,
                                       const Point &segment_b) {
    if constexpr (is_integral<coordinate_t<Point>>::value) {
        return line_edge_intersection_integral(line_a, line_b, segment_a,
                                               segment_b);
    } else {
        return line_edge_intersection_floating(line_a, line_b, segment_a,
                                               segment_b);
    }
}
} // namespace line_convex_polygon_intersection_internal

template <class Point>
LinePolygonIntersectionResult<Point>
line_polygon_intersection(const std::vector<Point> &polygon,
                          const Point &line_a, const Point &line_b) {
    namespace lpi = line_convex_polygon_intersection_internal;
    using Coord = lpi::coordinate_t<Point>;
    using Calc = lpi::calc_t<Point>;

    static_assert(!lpi::is_integral<Coord>::value ||
                      lpi::is_signed<Coord>::value,
                  "integral coordinate type must be signed");

    const int n = static_cast<int>(polygon.size());
    assert(n >= 3);
    assert(!lpi::equals_value(static_cast<Calc>(lpi::get_x(line_a)),
                              static_cast<Calc>(lpi::get_x(line_b))) ||
           !lpi::equals_value(static_cast<Calc>(lpi::get_y(line_a)),
                              static_cast<Calc>(lpi::get_y(line_b))));

    const Point direction = line_b - line_a;
    auto vertex = [&](int index) -> const Point & {
        return polygon[lpi::positive_mod(index, n)];
    };
    auto height = [&](int index) -> Calc {
        return lpi::cross(vertex(index) - line_a, direction);
    };
    auto chain_index = [&](int start, int step, int offset) {
        return lpi::positive_mod(start + step * offset, n);
    };

    const auto [minimum_index, maximum_index] =
        lpi::find_extreme_vertices<Point>(n, height);
    const Calc minimum_value = height(minimum_index);
    const Calc maximum_value = height(maximum_index);

    LinePolygonIntersectionResult<Point> result;
    if (lpi::sign_value(minimum_value) > 0 ||
        lpi::sign_value(maximum_value) < 0) {
        return result;
    }

    const int forward_length =
        lpi::distance_forward(minimum_index, maximum_index, n);
    const int backward_length = n - forward_length;

    if (lpi::sign_value(minimum_value) == 0) {
        const int forward_zero =
            lpi::last_non_positive(forward_length, [&](int offset) {
                return height(minimum_index + offset);
            });
        const int backward_zero =
            lpi::last_non_positive(backward_length, [&](int offset) {
                return height(minimum_index - offset);
            });
        const int forward_index = chain_index(minimum_index, +1, forward_zero);
        const int backward_index =
            chain_index(minimum_index, -1, backward_zero);
        result.push_back(lpi::make_vertex_result(vertex(forward_index)));
        if (forward_index != backward_index) {
            result.push_back(lpi::make_vertex_result(vertex(backward_index)));
        }
        return result;
    }

    if (lpi::sign_value(maximum_value) == 0) {
        const int forward_length_from_max =
            lpi::distance_forward(maximum_index, minimum_index, n);
        const int backward_length_from_max = n - forward_length_from_max;
        const int forward_zero =
            lpi::last_non_negative(forward_length_from_max, [&](int offset) {
                return height(maximum_index + offset);
            });
        const int backward_zero =
            lpi::last_non_negative(backward_length_from_max, [&](int offset) {
                return height(maximum_index - offset);
            });
        const int forward_index = chain_index(maximum_index, +1, forward_zero);
        const int backward_index =
            chain_index(maximum_index, -1, backward_zero);
        result.push_back(lpi::make_vertex_result(vertex(forward_index)));
        if (forward_index != backward_index) {
            result.push_back(lpi::make_vertex_result(vertex(backward_index)));
        }
        return result;
    }

    const int first_cross =
        lpi::first_non_negative(forward_length, [&](int offset) {
            return height(minimum_index + offset);
        });
    const int first_cross_index = chain_index(minimum_index, +1, first_cross);
    if (lpi::sign_value(height(first_cross_index)) == 0) {
        result.push_back(lpi::make_vertex_result(vertex(first_cross_index)));
    } else {
        const int prev_index = chain_index(minimum_index, +1, first_cross - 1);
        result.push_back(lpi::make_edge_result(
            line_a, line_b, vertex(prev_index), vertex(first_cross_index)));
    }

    const int second_cross =
        lpi::first_non_negative(backward_length, [&](int offset) {
            return height(minimum_index - offset);
        });
    const int second_cross_index = chain_index(minimum_index, -1, second_cross);
    if (lpi::sign_value(height(second_cross_index)) == 0) {
        result.push_back(lpi::make_vertex_result(vertex(second_cross_index)));
    } else {
        const int prev_index = chain_index(minimum_index, -1, second_cross - 1);
        result.push_back(lpi::make_edge_result(
            line_a, line_b, vertex(prev_index), vertex(second_cross_index)));
    }

    return result;
}


#line 11 "verify/aoj-cgl-4-c.test.cpp"

struct Point {
    long long x;
    long long y;

    Point() : x(0), y(0) {}
    Point(long long x_, long long y_) : x(x_), y(y_) {}

    Point operator+(const Point &other) const {
        return Point(x + other.x, y + other.y);
    }
    Point operator-(const Point &other) const {
        return Point(x - other.x, y - other.y);
    }
};

struct RealPoint {
    long double x;
    long double y;

    RealPoint() : x(0), y(0) {}
    RealPoint(long double x_, long double y_) : x(x_), y(y_) {}

    RealPoint operator+(const RealPoint &other) const {
        return RealPoint(x + other.x, y + other.y);
    }
    RealPoint operator-(const RealPoint &other) const {
        return RealPoint(x - other.x, y - other.y);
    }
    RealPoint operator*(long double k) const { return RealPoint(x * k, y * k); }
};

__int128_t cross_value(const Point &a, const Point &b) {
    return static_cast<__int128_t>(a.x) * static_cast<__int128_t>(b.y) -
           static_cast<__int128_t>(a.y) * static_cast<__int128_t>(b.x);
}

int sign_value(__int128_t x) {
    if (x < 0) {
        return -1;
    }
    if (x > 0) {
        return 1;
    }
    return 0;
}

Point read_point() {
    long long x, y;
    std::cin >> x >> y;
    return Point(x, y);
}

RealPoint to_real_point(const Point &p) {
    return RealPoint(static_cast<long double>(p.x),
                     static_cast<long double>(p.y));
}

long double cross_ld(const RealPoint &a, const RealPoint &b) {
    return a.x * b.y - a.y * b.x;
}

long double dot_ld(const RealPoint &a, const RealPoint &b) {
    return a.x * b.x + a.y * b.y;
}

bool same_point(const RealPoint &a, const RealPoint &b) {
    return std::abs(a.x - b.x) <= 1e-12L && std::abs(a.y - b.y) <= 1e-12L;
}

void add_unique(std::vector<RealPoint> &points, const RealPoint &point) {
    for (const auto &q : points) {
        if (same_point(q, point)) {
            return;
        }
    }
    points.push_back(point);
}

RealPoint line_intersection(const Point &line_a, const Point &line_b,
                            const Point &segment_a, const Point &segment_b) {
    const RealPoint u = to_real_point(segment_a);
    const RealPoint v = to_real_point(segment_b);
    const RealPoint dir = to_real_point(line_b - line_a);
    const RealPoint edge = v - u;
    const long double denominator = cross_ld(edge, dir);
    const long double numerator =
        cross_ld(to_real_point(line_a - segment_a), dir);
    const long double t = numerator / denominator;
    return u + edge * t;
}

std::vector<RealPoint> brute_intersections(const std::vector<Point> &polygon,
                                           const Point &line_a,
                                           const Point &line_b) {
    const int n = static_cast<int>(polygon.size());
    const Point direction = line_b - line_a;
    std::vector<RealPoint> points;
    for (int i = 0; i < n; ++i) {
        const Point &u = polygon[i];
        const Point &v = polygon[(i + 1) % n];
        const __int128_t hu = cross_value(u - line_a, direction);
        const __int128_t hv = cross_value(v - line_a, direction);
        const int su = sign_value(hu);
        const int sv = sign_value(hv);
        if (su == 0 && sv == 0) {
            add_unique(points, to_real_point(u));
            add_unique(points, to_real_point(v));
        } else if (su == 0) {
            add_unique(points, to_real_point(u));
        } else if (sv == 0) {
            add_unique(points, to_real_point(v));
        } else if (su != sv) {
            add_unique(points, line_intersection(line_a, line_b, u, v));
        }
    }
    if (points.size() <= 2) {
        return points;
    }
    const RealPoint dir = to_real_point(direction);
    int min_index = 0;
    int max_index = 0;
    for (int i = 1; i < static_cast<int>(points.size()); ++i) {
        if (dot_ld(points[i], dir) < dot_ld(points[min_index], dir)) {
            min_index = i;
        }
        if (dot_ld(points[i], dir) > dot_ld(points[max_index], dir)) {
            max_index = i;
        }
    }
    return {points[min_index], points[max_index]};
}

std::vector<RealPoint> convex_cut(const std::vector<Point> &polygon,
                                  const Point &line_a, const Point &line_b) {
    const int n = static_cast<int>(polygon.size());
    const Point direction = line_b - line_a;
    std::vector<RealPoint> result;
    for (int i = 0; i < n; ++i) {
        const Point &now = polygon[i];
        const Point &next = polygon[(i + 1) % n];

        // AOJ CGL_4_C は有向直線 line_a -> line_b の左側を残す。
        const __int128_t current = cross_value(direction, now - line_a);
        const __int128_t next_value = cross_value(direction, next - line_a);

        if (sign_value(current) >= 0) {
            result.push_back(to_real_point(now));
        }
        if ((current < 0 && next_value > 0) ||
            (current > 0 && next_value < 0)) {
            result.push_back(line_intersection(line_a, line_b, now, next));
        }
    }
    return result;
}

long double polygon_area(const std::vector<RealPoint> &polygon) {
    const int n = static_cast<int>(polygon.size());
    long double area = 0;
    for (int i = 0; i < n; ++i) {
        area += cross_ld(polygon[i], polygon[(i + 1) % n]);
    }
    return std::abs(area) / 2;
}

std::vector<Point>
remove_collinear_vertices(const std::vector<Point> &polygon) {
    const int n = static_cast<int>(polygon.size());
    std::vector<Point> result;
    result.reserve(n);
    for (int i = 0; i < n; ++i) {
        const Point &prev = polygon[(i + n - 1) % n];
        const Point &cur = polygon[i];
        const Point &next = polygon[(i + 1) % n];
        if (sign_value(cross_value(cur - prev, next - cur)) == 0) {
            continue;
        }
        result.push_back(cur);
    }
    assert(result.size() >= 3);
    return result;
}

void verify_intersection(const std::vector<Point> &polygon, const Point &line_a,
                         const Point &line_b) {
    const auto fast = line_polygon_intersection(polygon, line_a, line_b);
    const auto slow = brute_intersections(polygon, line_a, line_b);
    assert(fast.size() == slow.size());
    if (fast.empty()) {
        return;
    }

    auto to_real = [](const auto &point) {
        return RealPoint(point.template x_as<long double>(),
                         point.template y_as<long double>());
    };

    if (fast.size() == 1) {
        assert(same_point(to_real(fast[0]), slow[0]));
        return;
    }
    const RealPoint p0 = to_real(fast[0]);
    const RealPoint p1 = to_real(fast[1]);
    const bool ok = (same_point(p0, slow[0]) && same_point(p1, slow[1])) ||
                    (same_point(p0, slow[1]) && same_point(p1, slow[0]));
    assert(ok);
}

int main() {
    int n;
    std::cin >> n;
    std::vector<Point> polygon(n);
    for (auto &p : polygon) {
        p = read_point();
    }
    polygon = remove_collinear_vertices(polygon);

    int q;
    std::cin >> q;
    std::cout << std::fixed << std::setprecision(20);
    while (q--) {
        const Point line_a = read_point();
        const Point line_b = read_point();
        verify_intersection(polygon, line_a, line_b);
        std::cout << polygon_area(convex_cut(polygon, line_a, line_b)) << '\n';
    }
}

Test cases

Env Name Status Elapsed Memory
g++ 00_sample_00 :heavy_check_mark: AC 2 ms 4 MB
g++ 01_small_00 :heavy_check_mark: AC 2 ms 4 MB
g++ 01_small_01 :heavy_check_mark: AC 2 ms 4 MB
g++ 02_empty_00 :heavy_check_mark: AC 2 ms 4 MB
g++ 03_general_00 :heavy_check_mark: AC 2 ms 4 MB
g++ 04_minimum_00 :heavy_check_mark: AC 2 ms 4 MB
g++ 05_extreme_00 :heavy_check_mark: AC 2 ms 4 MB
g++ 06_rand_00 :heavy_check_mark: AC 3 ms 4 MB
g++ 06_rand_01 :heavy_check_mark: AC 3 ms 4 MB
g++ 06_rand_02 :heavy_check_mark: AC 2 ms 4 MB
g++ 06_rand_03 :heavy_check_mark: AC 2 ms 4 MB
g++ 06_rand_04 :heavy_check_mark: AC 3 ms 4 MB
g++ 06_rand_05 :heavy_check_mark: AC 4 ms 4 MB
g++ 06_rand_06 :heavy_check_mark: AC 3 ms 4 MB
g++ 06_rand_07 :heavy_check_mark: AC 3 ms 4 MB
g++ 07_corner_00 :heavy_check_mark: AC 2 ms 4 MB
clang++ 00_sample_00 :heavy_check_mark: AC 2 ms 4 MB
clang++ 01_small_00 :heavy_check_mark: AC 2 ms 4 MB
clang++ 01_small_01 :heavy_check_mark: AC 2 ms 4 MB
clang++ 02_empty_00 :heavy_check_mark: AC 2 ms 4 MB
clang++ 03_general_00 :heavy_check_mark: AC 2 ms 4 MB
clang++ 04_minimum_00 :heavy_check_mark: AC 2 ms 4 MB
clang++ 05_extreme_00 :heavy_check_mark: AC 2 ms 4 MB
clang++ 06_rand_00 :heavy_check_mark: AC 3 ms 4 MB
clang++ 06_rand_01 :heavy_check_mark: AC 3 ms 4 MB
clang++ 06_rand_02 :heavy_check_mark: AC 2 ms 4 MB
clang++ 06_rand_03 :heavy_check_mark: AC 3 ms 4 MB
clang++ 06_rand_04 :heavy_check_mark: AC 3 ms 4 MB
clang++ 06_rand_05 :heavy_check_mark: AC 3 ms 4 MB
clang++ 06_rand_06 :heavy_check_mark: AC 3 ms 4 MB
clang++ 06_rand_07 :heavy_check_mark: AC 3 ms 4 MB
clang++ 07_corner_00 :heavy_check_mark: AC 2 ms 4 MB
Back to top page