2D Point and Predicates
(geometry/point.hpp)

• View this file on GitHub
• Last update: 2026-10-05 22:23:07+09:00
• Include: #include "geometry/point.hpp"

Overview

Point<T> is the base type for the 2D geometry library. It provides vector arithmetic, lexicographic comparison, dot and cross products, distances, orientation, rotation, normalization, and its trivial centroid.

For integral coordinates, dot products, cross products, squared norms, and orientation calculations use signed 128-bit arithmetic. For floating-point coordinates, predicates use long double and accept a scale-aware epsilon.

Coordinate<T> also accepts exact numeric classes such as math::Rational, including Rational<utilities::BigInt>. A custom coordinate type must be copyable, totally ordered, constructible from 0 and 1, support unary signs and closed +, -, *, /, +=, and -=, and provide an explicit conversion to long double. Its arithmetic and ordering must be exact. bool is excluded.

ExactCoordinate<T> means Coordinate<T> && !std::floating_point<T>. wide_type<T> is __int128_t for built-in integers, long double for built-in floating-point types, and T for custom exact types. Rational dot products, cross products, squared distances, and predicates therefore stay rational. Exact predicates ignore eps.

Inputs and intermediate results must fit the chosen arithmetic type. For Rational<long long>, every normalized intermediate fraction must fit its underlying integer type; use Rational<BigInt> when that is insufficient. Complexities below count scalar operations as constant time. Rational operations add the arithmetic and gcd costs described on the rational documentation page.

Point

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

Point supports equality, lexicographic <, unary signs, addition, subtraction, scalar multiplication, and scalar division. Scalar operations return Point<std::common_type_t<T, Scalar>>. Integer or rational scalars keep rational coordinates exact; a floating-point scalar mixed with rational coordinates produces Point<long double>.

Conversions to other point types are explicit. Functions whose signatures return long double or Point<long double> still return approximations with rational input, including distances, projections, intersection coordinates, rotation, normalization, centroids, and division points.

Functions

Function	Description	Complexity
int sign(wide_type<T> value, long double eps = 1e-12L)	Returns the scalar sign; exact types ignore eps.	$O(1)$
Point<long double> centroid(const Point<T>& point)	Returns the point itself as Point<long double>.	$O(1)$
wide_type<T> dot(const Point<T>& a, const Point<T>& b)	Dot product.	$O(1)$
wide_type<T> cross(const Point<T>& a, const Point<T>& b)	Cross product of two vectors.	$O(1)$
wide_type<T> cross(const Point<T>& origin, const Point<T>& a, const Point<T>& b)	Cross product of vectors a - origin and b - origin.	$O(1)$
wide_type<T> norm2(const Point<T>& point)	Squared Euclidean norm.	$O(1)$
wide_type<T> distance2(const Point<T>& a, const Point<T>& b)	Squared Euclidean distance.	$O(1)$
long double norm(const Point<T>& point)	Euclidean norm as long double.	$O(1)$
long double distance(const Point<T>& a, const Point<T>& b)	Euclidean distance as long double.	$O(1)$
Point<long double> internal_division_point(const Point<T>& a, const Point<T>& b, M m, N n)	Returns the point internally dividing AB in ratio AP:PB = m:n.	$O(1)$
Point<long double> external_division_point(const Point<T>& a, const Point<T>& b, M m, N n)	Returns the point externally dividing AB in ratio AP:PB = m:n.	$O(1)$
int orientation(const Point<T>& a, const Point<T>& b, const Point<T>& c, long double eps = 1e-12L)	Returns 1 for counterclockwise, -1 for clockwise, and 0 for collinear.	$O(1)$
bool collinear(const Point<T>& a, const Point<T>& b, const Point<T>& c, long double eps = 1e-12L)	Returns whether three points are collinear.	$O(1)$
Point<long double> rotate(const Point<T>& point, long double angle)	Rotates p counterclockwise by radians.	$O(1)$
Point<long double> normalized(const Point<T>& point)	Returns a unit vector in p’s direction.	$O(1)$

normalized requires a nonzero vector.

Floating-point tolerance

For floating-point orientation and collinear, eps is a dimensionless relative tolerance. If u = b - a, v = c - a, and s(w) = max(abs(w.x), abs(w.y)), their determinant is treated as zero when

\[|u_xv_y-u_yv_x| \leq \mathtt{eps}\,s(u)\max(s(u),s(v)).\]

Consequently, uniformly scaling all coordinates does not change the result. The s(u)^2 floor also absorbs roundoff near either endpoint of the directed baseline a–b. Passing eps = 0 requests a strict comparison of the computed long double determinant. Integral and custom exact predicates ignore eps.

The lower-level sign(value, eps) function retains absolute-tolerance semantics because it has no operands from which to derive a scale.

Internal and external division

For internal_division_point(a, b, m, n), the returned point $P$ satisfies $AP:PB=m:n$ and is computed as

\[P = A + \frac{m}{m+n}(B-A).\]

For positive m and n, P lies between a and b. The function requires m + n != 0.

external_division_point uses the same ratio convention:

\[P = A + \frac{m}{m-n}(B-A).\]

For positive unequal ratios, the point lies outside the segment. It is beyond b when m > n and beyond a when m < n. The function requires m != n. Both functions return Point<long double>.

Example

#include "geometry/point.hpp"

#include <iostream>

int main() {
    using Point = m1une::geometry::Point<long long>;
    Point a(0, 0);
    Point b(3, 0);
    Point c(1, 2);

    std::cout << m1une::geometry::orientation(a, b, c) << "\n"; // 1
    std::cout << m1une::geometry::cross(a, b, c) << "\n";       // 6
}

Rational coordinates

#include "geometry/convex_hull.hpp"
#include "math/rational.hpp"

using Fraction = m1une::math::Rational<>;
using Point = m1une::geometry::Point<Fraction>;

Point a(0, 0);
Point b(Fraction(1, 2), 0);
Point c(0, Fraction(1, 3));
auto twice_b = b * 2;                    // Point<Fraction>(1, 0)
auto area2 = m1une::geometry::cross(a, b, c); // Fraction(1, 6)
auto hull = m1une::geometry::convex_hull(std::vector<Point>{a, b, c});

Depends on

• geometry/detail/floating_predicate.hpp

Required by

• Geometry Bundle (geometry/all.hpp)
• Geometry Bundle (geometry/all.hpp)
• Angle Sort (geometry/angle_sort.hpp)
• Circles (geometry/circle.hpp)
• Circle Coverage Areas (geometry/circle_coverage_areas.hpp)
• Area of Union of Circles (geometry/circle_union_area.hpp)
• Closest Pair of Points (geometry/closest_pair.hpp)
• Convex Decomposition (geometry/convex_decomposition.hpp)
• Convex Hull (geometry/convex_hull.hpp)
• Convex Layers (geometry/convex_layers.hpp)
• Convex Polygons (geometry/convex_polygon.hpp)
• Convex Polygons (geometry/convex_polygon.hpp)
• Count Points in Triangle (geometry/count_points_in_triangle.hpp)
• Delaunay Triangulation (geometry/delaunay_triangulation.hpp)
• geometry/detail/convex_polygon_normalize.hpp
• Euclidean Minimum Spanning Tree (geometry/euclidean_mst.hpp)
• Farthest Pair of Points (geometry/farthest_pair.hpp)
• Half-Plane Intersection (geometry/half_plane_intersection.hpp)
• Linear Objects (geometry/linear.hpp)
• Manhattan Minimum Spanning Tree (geometry/manhattan_mst.hpp)
• Manhattan Segment Intersections (geometry/manhattan_segment_intersections.hpp)
• Minimum Enclosing Circle (geometry/minimum_enclosing_circle.hpp)
• Minkowski Sum (geometry/minkowski_sum.hpp)
• Perpendicular Bisector (geometry/perpendicular_bisector.hpp)
• Polygons (geometry/polygon.hpp)
• Area of Union of Rectangles (geometry/rectangle_union_area.hpp)
• Steiner Convex Decomposition (geometry/steiner_convex_decomposition.hpp)
• Voronoi Diagram (geometry/voronoi_diagram.hpp)

Verified with

• verify/geometry/angle_sort.test.cpp
• verify/geometry/centroid.test.cpp
• verify/geometry/centroid.test.cpp
• verify/geometry/circle_boundary_intersection.test.cpp
• verify/geometry/circle_circle_intersection.test.cpp
• verify/geometry/circle_circle_intersection_area.test.cpp
• verify/geometry/circle_common_tangents.test.cpp
• verify/geometry/circle_coverage_areas.test.cpp
• verify/geometry/circle_filled.test.cpp
• verify/geometry/circle_line_intersection.test.cpp
• verify/geometry/circle_operations.test.cpp
• verify/geometry/circle_polygon_intersection_area.test.cpp
• verify/geometry/circle_ray.test.cpp
• verify/geometry/circle_relation.test.cpp
• verify/geometry/circle_tangent_points.test.cpp
• verify/geometry/circle_union_area.test.cpp
• verify/geometry/circumcircle.test.cpp
• verify/geometry/closest_pair.test.cpp
• verify/geometry/closest_points.test.cpp
• verify/geometry/convex_decomposition.test.cpp
• verify/geometry/convex_decomposition.test.cpp
• verify/geometry/convex_diameter.test.cpp
• verify/geometry/convex_diameter.test.cpp
• verify/geometry/convex_hull.test.cpp
• verify/geometry/convex_layers.test.cpp
• verify/geometry/convex_layers.test.cpp
• verify/geometry/convex_polygon.test.cpp
• verify/geometry/convex_polygon.test.cpp
• verify/geometry/count_points_in_triangle.test.cpp
• verify/geometry/delaunay_triangulation.test.cpp
• verify/geometry/euclidean_mst.test.cpp
• verify/geometry/farthest_pair.test.cpp
• verify/geometry/floating_predicates.test.cpp
• verify/geometry/geometry_algorithms.test.cpp
• verify/geometry/geometry_algorithms.test.cpp
• verify/geometry/half_plane_intersection.test.cpp
• verify/geometry/half_plane_intersection_random.test.cpp
• verify/geometry/incircle.test.cpp
• verify/geometry/is_convex_polygon.test.cpp
• verify/geometry/is_convex_polygon.test.cpp
• verify/geometry/linear_intersection.test.cpp
• verify/geometry/manhattan_mst.test.cpp
• verify/geometry/manhattan_segment_intersections.test.cpp
• verify/geometry/minimum_enclosing_circle.test.cpp
• verify/geometry/minkowski_sum.test.cpp
• verify/geometry/minkowski_sum.test.cpp
• verify/geometry/perpendicular_bisector.test.cpp
• verify/geometry/point_in_polygon.test.cpp
• verify/geometry/polygon_area.test.cpp
• verify/geometry/polygon_clipping.test.cpp
• verify/geometry/polygon_filled.test.cpp
• verify/geometry/polygon_operations.test.cpp
• verify/geometry/polygon_operations.test.cpp
• verify/geometry/projection.test.cpp
• verify/geometry/rational.test.cpp
• verify/geometry/rational.test.cpp
• verify/geometry/ray.test.cpp
• verify/geometry/rectangle_union_area.test.cpp
• verify/geometry/segment_intersection.test.cpp
• verify/geometry/segment_intersection_point.test.cpp
• verify/geometry/steiner_convex_decomposition.test.cpp
• verify/geometry/steiner_convex_decomposition.test.cpp
• verify/geometry/voronoi_diagram.test.cpp

Code

#ifndef M1UNE_GEOMETRY_POINT_HPP
#define M1UNE_GEOMETRY_POINT_HPP 1

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

#include "detail/floating_predicate.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

#endif  // M1UNE_GEOMETRY_POINT_HPP

#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

Back to top page