m1une's library

This documentation is automatically generated by online-judge-tools/verification-helper

View on GitHub

:heavy_check_mark: verify/geometry/minkowski_sum.test.cpp

Depends on

Code

#define PROBLEM "https://judge.yosupo.jp/problem/aplusb"

#include "../../geometry/convex_hull.hpp"
#include "../../geometry/minkowski_sum.hpp"

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <utility>
#include <vector>

namespace {

using m1une::geometry::Point;
using PointType = Point<long long>;

std::vector<PointType> brute_sum(
    const std::vector<PointType>& first,
    const std::vector<PointType>& second
) {
    std::vector<PointType> sums;
    sums.reserve(first.size() * second.size());
    for (const PointType& left : first) {
        for (const PointType& right : second) {
            sums.push_back(left + right);
        }
    }
    return m1une::geometry::convex_hull(std::move(sums));
}

void test_fixed() {
    std::vector<PointType> square;
    square.emplace_back(0, 0);
    square.emplace_back(2, 0);
    square.emplace_back(2, 2);
    square.emplace_back(0, 2);

    std::vector<PointType> triangle;
    triangle.emplace_back(0, 0);
    triangle.emplace_back(3, 0);
    triangle.emplace_back(0, 1);
    assert(
        m1une::geometry::minkowski_sum(square, triangle) ==
        brute_sum(square, triangle)
    );

    std::reverse(square.begin(), square.end());
    triangle.push_back(triangle.front());
    assert(
        m1une::geometry::minkowski_sum(square, triangle) ==
        brute_sum(square, triangle)
    );

    std::vector<PointType> segment;
    segment.emplace_back(3, 0);
    segment.emplace_back(0, 0);
    std::vector<PointType> point;
    point.emplace_back(2, 4);
    std::vector<PointType> expected;
    expected.emplace_back(2, 4);
    expected.emplace_back(5, 4);
    assert(m1une::geometry::minkowski_sum(segment, point) == expected);

    std::vector<PointType> redundant;
    redundant.emplace_back(0, 0);
    redundant.emplace_back(1, 0);
    redundant.emplace_back(2, 0);
    redundant.emplace_back(2, 2);
    redundant.emplace_back(0, 2);
    redundant.emplace_back(0, 0);
    assert(
        m1une::geometry::minkowski_sum(redundant, point) ==
        brute_sum(redundant, point)
    );
}

void test_randomized() {
    std::uint64_t state = 0x243f6a8885a308d3ULL;
    auto random = [&state]() {
        state ^= state << 7;
        state ^= state >> 9;
        return state;
    };

    for (int trial = 0; trial < 5000; ++trial) {
        std::vector<PointType> first_points;
        std::vector<PointType> second_points;
        const int first_size = 1 + int(random() % 12);
        const int second_size = 1 + int(random() % 12);
        for (int index = 0; index < first_size; ++index) {
            first_points.emplace_back(
                static_cast<long long>(random() % 31) - 15,
                static_cast<long long>(random() % 31) - 15
            );
        }
        for (int index = 0; index < second_size; ++index) {
            second_points.emplace_back(
                static_cast<long long>(random() % 31) - 15,
                static_cast<long long>(random() % 31) - 15
            );
        }
        std::vector<PointType> first =
            m1une::geometry::convex_hull(first_points);
        std::vector<PointType> second =
            m1une::geometry::convex_hull(second_points);
        if (random() & 1) std::reverse(first.begin(), first.end());
        if (random() & 1) std::reverse(second.begin(), second.end());

        const std::vector<PointType> expected = brute_sum(first, second);
        const std::vector<PointType> actual =
            m1une::geometry::minkowski_sum(first, second);
        assert(m1une::geometry::convex_hull(actual) == expected);
    }
}

}  // namespace

int main() {
    test_fixed();
    test_randomized();

    long long first;
    long long second;
    std::cin >> first >> second;
    std::cout << first + second << '\n';
}
#line 1 "verify/geometry/minkowski_sum.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/aplusb"

#line 1 "geometry/convex_hull.hpp"



#include <algorithm>
#include <cstddef>
#include <utility>
#include <vector>

#line 1 "geometry/point.hpp"



#include <cmath>
#include <concepts>
#include <cassert>
#include <type_traits>

#line 1 "geometry/detail/floating_predicate.hpp"



namespace m1une {
namespace geometry {
namespace predicate_detail {

template <typename T>
constexpr T absolute(T value) {
    return value < T(0) ? -value : value;
}

template <typename T>
constexpr T max_value(T first, T second) {
    return first < second ? second : first;
}

template <typename T>
constexpr T vector_scale(T x, T y) {
    return max_value(absolute(x), absolute(y));
}

template <bool Exact, typename T>
constexpr int scaled_sign(T value, T scale, long double eps) {
    if constexpr (Exact) {
        return (value > T(0)) - (value < T(0));
    } else {
        const T tolerance = T(eps) * scale;
        return (value > tolerance) - (value < -tolerance);
    }
}

template <bool Exact, typename T>
constexpr T determinant_scale(T ax, T ay, T bx, T by) {
    if constexpr (Exact) {
        return T(0);
    } else {
        return vector_scale(ax, ay) * vector_scale(bx, by);
    }
}

template <bool Exact, typename T>
constexpr int determinant_sign(
    T ax,
    T ay,
    T bx,
    T by,
    long double eps
) {
    const T determinant = ax * by - ay * bx;
    return scaled_sign<Exact>(
        determinant,
        determinant_scale<Exact>(ax, ay, bx, by),
        eps
    );
}

template <bool Exact, typename T>
constexpr int orientation_sign(
    T direction_x,
    T direction_y,
    T offset_x,
    T offset_y,
    long double eps
) {
    const T determinant =
        direction_x * offset_y - direction_y * offset_x;
    T scale = T(0);
    if constexpr (!Exact) {
        const T direction_scale =
            vector_scale(direction_x, direction_y);
        scale = direction_scale * max_value(
            direction_scale,
            vector_scale(offset_x, offset_y)
        );
    }
    return scaled_sign<Exact>(determinant, scale, eps);
}

template <bool Exact, typename T>
constexpr int dot_sign(
    T ax,
    T ay,
    T bx,
    T by,
    long double eps
) {
    const T value = ax * bx + ay * by;
    T scale = T(0);
    if constexpr (!Exact) {
        scale = vector_scale(ax, ay) * vector_scale(bx, by);
    }
    return scaled_sign<Exact>(value, scale, eps);
}

}  // namespace predicate_detail
}  // namespace geometry
}  // namespace m1une


#line 10 "geometry/point.hpp"

namespace m1une {
namespace geometry {

template <typename T>
concept Coordinate = !std::same_as<std::remove_cv_t<T>, bool> &&
    (std::is_arithmetic_v<T> ||
     (std::copyable<T> && std::totally_ordered<T> && requires(T a, T b) {
         T(0);
         T(1);
         static_cast<long double>(a);
         { +a } -> std::same_as<T>;
         { -a } -> std::same_as<T>;
         { a + b } -> std::same_as<T>;
         { a - b } -> std::same_as<T>;
         { a * b } -> std::same_as<T>;
         { a / b } -> std::same_as<T>;
         { a += b } -> std::same_as<T&>;
         { a -= b } -> std::same_as<T&>;
     }));

// Custom coordinate types keep their own exact arithmetic.
template <typename T>
concept ExactCoordinate = Coordinate<T> && !std::floating_point<T>;

template <Coordinate T>
using wide_type = std::conditional_t<std::integral<T>, __int128_t,
    std::conditional_t<std::floating_point<T>, long double, T>>;

template <Coordinate T>
struct Point {
    T x;
    T y;

    constexpr Point() : x(0), y(0) {}
    constexpr Point(T x_value, T y_value) : x(x_value), y(y_value) {}

    template <Coordinate U>
    explicit constexpr Point(const Point<U>& other)
        : x(static_cast<T>(other.x)), y(static_cast<T>(other.y)) {}

    constexpr Point& operator+=(const Point& other) {
        x += other.x;
        y += other.y;
        return *this;
    }

    constexpr Point& operator-=(const Point& other) {
        x -= other.x;
        y -= other.y;
        return *this;
    }

    constexpr Point operator+() const {
        return *this;
    }

    constexpr Point operator-() const {
        return Point(-x, -y);
    }

    friend constexpr Point operator+(Point left, const Point& right) {
        return left += right;
    }

    friend constexpr Point operator-(Point left, const Point& right) {
        return left -= right;
    }

    friend constexpr bool operator==(const Point&, const Point&) = default;

    friend constexpr bool operator<(const Point& left, const Point& right) {
        if (left.x != right.x) return left.x < right.x;
        return left.y < right.y;
    }
};

template <Coordinate T>
constexpr Point<long double> centroid(const Point<T>& point) {
    return Point<long double>(point);
}

template <Coordinate T, typename Scalar>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
constexpr auto operator*(const Point<T>& point, Scalar scalar) {
    using Result = std::common_type_t<T, Scalar>;
    return Point<Result>(
        Result(point.x) * Result(scalar),
        Result(point.y) * Result(scalar)
    );
}

template <typename Scalar, Coordinate T>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
constexpr auto operator*(Scalar scalar, const Point<T>& point) {
    return point * scalar;
}

template <Coordinate T, typename Scalar>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
constexpr auto operator/(const Point<T>& point, Scalar scalar) {
    using Result = std::common_type_t<T, Scalar>;
    return Point<Result>(
        Result(point.x) / Result(scalar),
        Result(point.y) / Result(scalar)
    );
}

template <Coordinate T>
constexpr wide_type<T> dot(const Point<T>& a, const Point<T>& b) {
    using W = wide_type<T>;
    return W(a.x) * W(b.x) + W(a.y) * W(b.y);
}

template <Coordinate T>
constexpr wide_type<T> cross(const Point<T>& a, const Point<T>& b) {
    using W = wide_type<T>;
    return W(a.x) * W(b.y) - W(a.y) * W(b.x);
}

template <Coordinate T>
constexpr wide_type<T> cross(
    const Point<T>& origin,
    const Point<T>& a,
    const Point<T>& b
) {
    using W = wide_type<T>;
    W ax = W(a.x) - W(origin.x);
    W ay = W(a.y) - W(origin.y);
    W bx = W(b.x) - W(origin.x);
    W by = W(b.y) - W(origin.y);
    return ax * by - ay * bx;
}

template <Coordinate T>
constexpr wide_type<T> norm2(const Point<T>& point) {
    return dot(point, point);
}

template <Coordinate T>
constexpr wide_type<T> distance2(const Point<T>& a, const Point<T>& b) {
    using W = wide_type<T>;
    W dx = W(a.x) - W(b.x);
    W dy = W(a.y) - W(b.y);
    return dx * dx + dy * dy;
}

template <Coordinate T>
long double norm(const Point<T>& point) {
    return std::hypot(
        static_cast<long double>(point.x),
        static_cast<long double>(point.y)
    );
}

template <Coordinate T>
long double distance(const Point<T>& a, const Point<T>& b) {
    return std::hypot(
        static_cast<long double>(a.x) - static_cast<long double>(b.x),
        static_cast<long double>(a.y) - static_cast<long double>(b.y)
    );
}

template <Coordinate T, typename M, typename N>
requires (std::is_arithmetic_v<M> || Coordinate<M>) &&
         (std::is_arithmetic_v<N> || Coordinate<N>)
constexpr Point<long double> internal_division_point(
    const Point<T>& a,
    const Point<T>& b,
    M m,
    N n
) {
    long double first_ratio = static_cast<long double>(m);
    long double second_ratio = static_cast<long double>(n);
    long double denominator = first_ratio + second_ratio;
    assert(denominator != 0);
    Point<long double> first(a);
    Point<long double> direction = Point<long double>(b) - first;
    return first + direction * (first_ratio / denominator);
}

template <Coordinate T, typename M, typename N>
requires (std::is_arithmetic_v<M> || Coordinate<M>) &&
         (std::is_arithmetic_v<N> || Coordinate<N>)
constexpr Point<long double> external_division_point(
    const Point<T>& a,
    const Point<T>& b,
    M m,
    N n
) {
    long double first_ratio = static_cast<long double>(m);
    long double second_ratio = static_cast<long double>(n);
    long double denominator = first_ratio - second_ratio;
    assert(denominator != 0);
    Point<long double> first(a);
    Point<long double> direction = Point<long double>(b) - first;
    return first + direction * (first_ratio / denominator);
}

template <Coordinate T>
constexpr int sign(wide_type<T> value, long double eps = 1e-12L) {
    return predicate_detail::scaled_sign<ExactCoordinate<T>>(
        value,
        wide_type<T>(1),
        eps
    );
}

template <Coordinate T>
constexpr int orientation(
    const Point<T>& a,
    const Point<T>& b,
    const Point<T>& c,
    long double eps = 1e-12L
) {
    using W = wide_type<T>;
    const W first_x = W(b.x) - W(a.x);
    const W first_y = W(b.y) - W(a.y);
    const W second_x = W(c.x) - W(a.x);
    const W second_y = W(c.y) - W(a.y);
    return predicate_detail::orientation_sign<ExactCoordinate<T>>(
        first_x,
        first_y,
        second_x,
        second_y,
        eps
    );
}

template <Coordinate T>
constexpr bool collinear(
    const Point<T>& a,
    const Point<T>& b,
    const Point<T>& c,
    long double eps = 1e-12L
) {
    return orientation(a, b, c, eps) == 0;
}

template <Coordinate T>
Point<long double> rotate(const Point<T>& point, long double angle) {
    long double cosine = std::cos(angle);
    long double sine = std::sin(angle);
    return Point<long double>(
        static_cast<long double>(point.x) * cosine -
            static_cast<long double>(point.y) * sine,
        static_cast<long double>(point.x) * sine +
            static_cast<long double>(point.y) * cosine
    );
}

template <Coordinate T>
Point<long double> normalized(const Point<T>& point) {
    long double length = norm(point);
    assert(length != 0);
    return Point<long double>(
        static_cast<long double>(point.x) / length,
        static_cast<long double>(point.y) / length
    );
}

}  // namespace geometry
}  // namespace m1une


#line 10 "geometry/convex_hull.hpp"

namespace m1une {
namespace geometry {

// Returns the convex hull counterclockwise from its lexicographically smallest
// point. The first point is not repeated at the end.
template <Coordinate T>
std::vector<Point<T>> convex_hull(
    std::vector<Point<T>> points,
    bool include_collinear = false
) {
    std::sort(points.begin(), points.end());
    points.erase(std::unique(points.begin(), points.end()), points.end());
    std::size_t size = points.size();
    if (size <= 1) return points;

    std::vector<Point<T>> hull;
    hull.reserve(2 * size);
    auto should_pop = [include_collinear](
        const Point<T>& first,
        const Point<T>& second,
        const Point<T>& third
    ) {
        int turn = orientation(first, second, third);
        return include_collinear ? turn < 0 : turn <= 0;
    };

    for (const Point<T>& point : points) {
        while (
            hull.size() >= 2 &&
            should_pop(hull[hull.size() - 2], hull.back(), point)
        ) {
            hull.pop_back();
        }
        hull.push_back(point);
    }

    std::size_t lower_size = hull.size();
    for (std::size_t index = size - 1; index-- > 0;) {
        const Point<T>& point = points[index];
        while (
            hull.size() > lower_size &&
            should_pop(hull[hull.size() - 2], hull.back(), point)
        ) {
            hull.pop_back();
        }
        hull.push_back(point);
    }
    hull.pop_back();

    if (include_collinear && hull.size() == 2 * points.size() - 2) {
        hull = std::move(points);
    }
    return hull;
}

}  // namespace geometry
}  // namespace m1une


#line 1 "geometry/minkowski_sum.hpp"



#line 8 "geometry/minkowski_sum.hpp"

#line 1 "geometry/detail/convex_polygon_normalize.hpp"



#line 8 "geometry/detail/convex_polygon_normalize.hpp"

#line 10 "geometry/detail/convex_polygon_normalize.hpp"

namespace m1une {
namespace geometry {
namespace convex_polygon_detail {

template <Coordinate T>
wide_type<T> boundary_area2(const std::vector<Point<T>>& polygon) {
    wide_type<T> result = 0;
    for (std::size_t index = 0; index < polygon.size(); ++index) {
        result += cross(
            polygon[index],
            polygon[(index + 1) % polygon.size()]
        );
    }
    return result;
}

template <Coordinate T>
std::vector<Point<T>> normalize_convex_boundary(
    std::vector<Point<T>> polygon,
    long double eps
) {
    if (polygon.size() >= 2 && polygon.front() == polygon.back()) {
        polygon.pop_back();
    }
    polygon.erase(
        std::unique(polygon.begin(), polygon.end()),
        polygon.end()
    );
    if (polygon.size() >= 2 && polygon.front() == polygon.back()) {
        polygon.pop_back();
    }
    if (polygon.size() <= 1) return polygon;
    if (
        polygon.size() >= 3 &&
        sign<T>(boundary_area2(polygon), eps) < 0
    ) {
        std::reverse(polygon.begin(), polygon.end());
    }

    const auto start = std::min_element(
        polygon.begin(),
        polygon.end(),
        [](const Point<T>& first, const Point<T>& second) {
            if (first.y != second.y) return first.y < second.y;
            return first.x < second.x;
        }
    );
    std::rotate(polygon.begin(), start, polygon.end());

    if (polygon.size() >= 3) {
        std::vector<Point<T>> cleaned;
        const std::size_t size = polygon.size();
        cleaned.reserve(size);
        for (std::size_t index = 0; index < size; ++index) {
            const Point<T>& previous = polygon[(index + size - 1) % size];
            const Point<T>& current = polygon[index];
            const Point<T>& next = polygon[(index + 1) % size];
            if (
                orientation(previous, current, next, eps) != 0 ||
                sign<T>(dot(current - previous, next - current), eps) < 0
            ) {
                cleaned.push_back(current);
            }
        }
        polygon = std::move(cleaned);
    }
    return polygon;
}

}  // namespace convex_polygon_detail
}  // namespace geometry
}  // namespace m1une


#line 10 "geometry/minkowski_sum.hpp"

namespace m1une {
namespace geometry {

// Returns the normalized boundary of the Minkowski sum of two nonempty
// ordered convex polygons.
template <Coordinate T>
std::vector<Point<T>> minkowski_sum(
    std::vector<Point<T>> first,
    std::vector<Point<T>> second,
    long double eps = 1e-12L
) {
    assert(!first.empty());
    assert(!second.empty());
    first = convex_polygon_detail::normalize_convex_boundary(
        std::move(first),
        eps
    );
    second = convex_polygon_detail::normalize_convex_boundary(
        std::move(second),
        eps
    );

    if (first.size() == 1 || second.size() == 1) {
        if (second.size() == 1) std::swap(first, second);
        for (Point<T>& point : second) point += first[0];
        return convex_polygon_detail::normalize_convex_boundary(
            std::move(second),
            eps
        );
    }

    std::vector<Point<T>> first_edges;
    std::vector<Point<T>> second_edges;
    first_edges.reserve(first.size());
    second_edges.reserve(second.size());
    for (std::size_t index = 0; index < first.size(); ++index) {
        first_edges.push_back(
            first[(index + 1) % first.size()] - first[index]
        );
    }
    for (std::size_t index = 0; index < second.size(); ++index) {
        second_edges.push_back(
            second[(index + 1) % second.size()] - second[index]
        );
    }

    Point<T> current = first.front() + second.front();
    std::vector<Point<T>> result;
    result.reserve(first.size() + second.size());
    result.push_back(current);
    std::size_t first_index = 0;
    std::size_t second_index = 0;
    while (
        first_index < first_edges.size() ||
        second_index < second_edges.size()
    ) {
        Point<T> step;
        if (first_index == first_edges.size()) {
            step = second_edges[second_index++];
        } else if (second_index == second_edges.size()) {
            step = first_edges[first_index++];
        } else {
            const auto turn = cross(
                first_edges[first_index],
                second_edges[second_index]
            );
            if (turn > 0) {
                step = first_edges[first_index++];
            } else if (turn < 0) {
                step = second_edges[second_index++];
            } else {
                step = first_edges[first_index++] +
                       second_edges[second_index++];
            }
        }
        current += step;
        if (
            first_index < first_edges.size() ||
            second_index < second_edges.size()
        ) {
            result.push_back(current);
        }
    }
    return convex_polygon_detail::normalize_convex_boundary(
        std::move(result),
        eps
    );
}

}  // namespace geometry
}  // namespace m1une


#line 5 "verify/geometry/minkowski_sum.test.cpp"

#line 8 "verify/geometry/minkowski_sum.test.cpp"
#include <cstdint>
#include <iostream>
#line 12 "verify/geometry/minkowski_sum.test.cpp"

namespace {

using m1une::geometry::Point;
using PointType = Point<long long>;

std::vector<PointType> brute_sum(
    const std::vector<PointType>& first,
    const std::vector<PointType>& second
) {
    std::vector<PointType> sums;
    sums.reserve(first.size() * second.size());
    for (const PointType& left : first) {
        for (const PointType& right : second) {
            sums.push_back(left + right);
        }
    }
    return m1une::geometry::convex_hull(std::move(sums));
}

void test_fixed() {
    std::vector<PointType> square;
    square.emplace_back(0, 0);
    square.emplace_back(2, 0);
    square.emplace_back(2, 2);
    square.emplace_back(0, 2);

    std::vector<PointType> triangle;
    triangle.emplace_back(0, 0);
    triangle.emplace_back(3, 0);
    triangle.emplace_back(0, 1);
    assert(
        m1une::geometry::minkowski_sum(square, triangle) ==
        brute_sum(square, triangle)
    );

    std::reverse(square.begin(), square.end());
    triangle.push_back(triangle.front());
    assert(
        m1une::geometry::minkowski_sum(square, triangle) ==
        brute_sum(square, triangle)
    );

    std::vector<PointType> segment;
    segment.emplace_back(3, 0);
    segment.emplace_back(0, 0);
    std::vector<PointType> point;
    point.emplace_back(2, 4);
    std::vector<PointType> expected;
    expected.emplace_back(2, 4);
    expected.emplace_back(5, 4);
    assert(m1une::geometry::minkowski_sum(segment, point) == expected);

    std::vector<PointType> redundant;
    redundant.emplace_back(0, 0);
    redundant.emplace_back(1, 0);
    redundant.emplace_back(2, 0);
    redundant.emplace_back(2, 2);
    redundant.emplace_back(0, 2);
    redundant.emplace_back(0, 0);
    assert(
        m1une::geometry::minkowski_sum(redundant, point) ==
        brute_sum(redundant, point)
    );
}

void test_randomized() {
    std::uint64_t state = 0x243f6a8885a308d3ULL;
    auto random = [&state]() {
        state ^= state << 7;
        state ^= state >> 9;
        return state;
    };

    for (int trial = 0; trial < 5000; ++trial) {
        std::vector<PointType> first_points;
        std::vector<PointType> second_points;
        const int first_size = 1 + int(random() % 12);
        const int second_size = 1 + int(random() % 12);
        for (int index = 0; index < first_size; ++index) {
            first_points.emplace_back(
                static_cast<long long>(random() % 31) - 15,
                static_cast<long long>(random() % 31) - 15
            );
        }
        for (int index = 0; index < second_size; ++index) {
            second_points.emplace_back(
                static_cast<long long>(random() % 31) - 15,
                static_cast<long long>(random() % 31) - 15
            );
        }
        std::vector<PointType> first =
            m1une::geometry::convex_hull(first_points);
        std::vector<PointType> second =
            m1une::geometry::convex_hull(second_points);
        if (random() & 1) std::reverse(first.begin(), first.end());
        if (random() & 1) std::reverse(second.begin(), second.end());

        const std::vector<PointType> expected = brute_sum(first, second);
        const std::vector<PointType> actual =
            m1une::geometry::minkowski_sum(first, second);
        assert(m1une::geometry::convex_hull(actual) == expected);
    }
}

}  // namespace

int main() {
    test_fixed();
    test_randomized();

    long long first;
    long long second;
    std::cin >> first >> second;
    std::cout << first + second << '\n';
}
Back to top page