m1une's library

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

View on GitHub

:heavy_check_mark: Angle Sort
(geometry/angle_sort.hpp)

Overview

This header sorts points counterclockwise by their angle around a chosen origin, without calling atan2.

sort_by_angle(points, origin, start);
auto sorted = angle_sorted(points, origin, start);

sort_by_angle modifies its argument. angle_sorted returns a sorted copy. Both take $O(N\log N)$ time.

Complexity Notation

sort_by_angle and angle_sorted take $O(N\log N)$ time. The comparator uses $O(1)$ additional memory per comparison.

Ordering Convention

The default, AngleSortStart::NegativeXAxis, sorts by the usual atan2(y, x) value from $-\pi$ to $\pi$. Imagine starting just below the negative x-axis and rotating counterclockwise:

Order Direction Angle
1 down-left close to $-\pi$
2 down $-\pi/2$
3 right $0$
4 up $\pi/2$
5 left $\pi$

In short:

down-left -> down -> right -> up -> up-left -> left

The origin is treated as having angle zero, so it appears with points pointing right. This convention matches Library Checker’s Sort Points by Argument problem.

AngleSortStart::PositiveXAxis uses the more common circular ordering that starts on the positive x-axis:

right -> up -> left -> down -> back to right

Points on the same ray are ordered by increasing distance from origin. Equivalent duplicate points may appear in either relative order.

Comparator

AngleLess<T> is the comparator used by the helpers:

m1une::geometry::AngleLess<long long> less(origin);
std::sort(points.begin(), points.end(), less);

For integral coordinates, comparisons use signed 128-bit arithmetic. As with the other exact geometry predicates, coordinate differences and products must fit that type.

Example

#include "geometry/angle_sort.hpp"

#include <iostream>
#include <vector>

int main() {
    using Point = m1une::geometry::Point<long long>;
    std::vector<Point> points;
    points.emplace_back(1, 0);
    points.emplace_back(0, 1);
    points.emplace_back(-1, 0);
    points.emplace_back(0, -1);

    m1une::geometry::sort_by_angle(
        points,
        Point(0, 0),
        m1une::geometry::AngleSortStart::PositiveXAxis
    );

    for (const Point& point : points) {
        std::cout << point.x << ' ' << point.y << "\n";
    }
}

Depends on

Required by

Verified with

Code

#ifndef M1UNE_GEOMETRY_ANGLE_SORT_HPP
#define M1UNE_GEOMETRY_ANGLE_SORT_HPP 1

#include <algorithm>
#include <vector>

#include "point.hpp"

namespace m1une {
namespace geometry {

enum class AngleSortStart {
    NegativeXAxis,
    PositiveXAxis,
};

template <Coordinate T>
struct AngleLess {
    Point<T> origin;
    AngleSortStart start;

    constexpr explicit AngleLess(
        Point<T> origin_value = Point<T>(),
        AngleSortStart start_value = AngleSortStart::NegativeXAxis
    ) : origin(origin_value), start(start_value) {}

    constexpr bool operator()(
        const Point<T>& first,
        const Point<T>& second
    ) const {
        using W = wide_type<T>;
        W first_x = W(first.x) - W(origin.x);
        W first_y = W(first.y) - W(origin.y);
        W second_x = W(second.x) - W(origin.x);
        W second_y = W(second.y) - W(origin.y);
        W first_distance = first_x * first_x + first_y * first_y;
        W second_distance = second_x * second_x + second_y * second_y;

        // atan2(0, 0) is treated as angle zero.
        if (first_distance == 0) first_x = 1;
        if (second_distance == 0) second_x = 1;

        auto half = [this](W x, W y) {
            if (start == AngleSortStart::PositiveXAxis) {
                return y < 0 || (y == 0 && x < 0);
            }
            return y > 0 || (y == 0 && x < 0);
        };

        bool first_half = half(first_x, first_y);
        bool second_half = half(second_x, second_y);
        if (first_half != second_half) return first_half < second_half;

        W product = first_x * second_y - first_y * second_x;
        if (product != 0) return product > 0;

        return first_distance < second_distance;
    }
};

// Sorts points counterclockwise by angle around `origin`.
template <Coordinate T>
void sort_by_angle(
    std::vector<Point<T>>& points,
    Point<T> origin = Point<T>(),
    AngleSortStart start = AngleSortStart::NegativeXAxis
) {
    std::sort(points.begin(), points.end(), AngleLess<T>(origin, start));
}

// Returns a counterclockwise angle-sorted copy.
template <Coordinate T>
std::vector<Point<T>> angle_sorted(
    std::vector<Point<T>> points,
    Point<T> origin = Point<T>(),
    AngleSortStart start = AngleSortStart::NegativeXAxis
) {
    sort_by_angle(points, origin, start);
    return points;
}

}  // namespace geometry
}  // namespace m1une

#endif  // M1UNE_GEOMETRY_ANGLE_SORT_HPP
#line 1 "geometry/angle_sort.hpp"



#include <algorithm>
#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 8 "geometry/angle_sort.hpp"

namespace m1une {
namespace geometry {

enum class AngleSortStart {
    NegativeXAxis,
    PositiveXAxis,
};

template <Coordinate T>
struct AngleLess {
    Point<T> origin;
    AngleSortStart start;

    constexpr explicit AngleLess(
        Point<T> origin_value = Point<T>(),
        AngleSortStart start_value = AngleSortStart::NegativeXAxis
    ) : origin(origin_value), start(start_value) {}

    constexpr bool operator()(
        const Point<T>& first,
        const Point<T>& second
    ) const {
        using W = wide_type<T>;
        W first_x = W(first.x) - W(origin.x);
        W first_y = W(first.y) - W(origin.y);
        W second_x = W(second.x) - W(origin.x);
        W second_y = W(second.y) - W(origin.y);
        W first_distance = first_x * first_x + first_y * first_y;
        W second_distance = second_x * second_x + second_y * second_y;

        // atan2(0, 0) is treated as angle zero.
        if (first_distance == 0) first_x = 1;
        if (second_distance == 0) second_x = 1;

        auto half = [this](W x, W y) {
            if (start == AngleSortStart::PositiveXAxis) {
                return y < 0 || (y == 0 && x < 0);
            }
            return y > 0 || (y == 0 && x < 0);
        };

        bool first_half = half(first_x, first_y);
        bool second_half = half(second_x, second_y);
        if (first_half != second_half) return first_half < second_half;

        W product = first_x * second_y - first_y * second_x;
        if (product != 0) return product > 0;

        return first_distance < second_distance;
    }
};

// Sorts points counterclockwise by angle around `origin`.
template <Coordinate T>
void sort_by_angle(
    std::vector<Point<T>>& points,
    Point<T> origin = Point<T>(),
    AngleSortStart start = AngleSortStart::NegativeXAxis
) {
    std::sort(points.begin(), points.end(), AngleLess<T>(origin, start));
}

// Returns a counterclockwise angle-sorted copy.
template <Coordinate T>
std::vector<Point<T>> angle_sorted(
    std::vector<Point<T>> points,
    Point<T> origin = Point<T>(),
    AngleSortStart start = AngleSortStart::NegativeXAxis
) {
    sort_by_angle(points, origin, start);
    return points;
}

}  // namespace geometry
}  // namespace m1une
Back to top page