Count Points in Triangle
(geometry/count_points_in_triangle.hpp)

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

Overview

T may be a built-in integer or an exact coordinate class such as math::Rational<long long> or math::Rational<utilities::BigInt>. Rational intermediate arithmetic uses T without floating-point conversion. All intermediate fractions must be representable. Listed complexities count scalar operations; rational arithmetic adds its gcd and integer arithmetic costs.

CountPointsInTriangle preprocesses two exact-coordinate point sets. Each query chooses three points from the first set and returns how many points from the second set lie strictly inside their triangle.

Queries take constant time after preprocessing, making the structure useful when the same candidate vertices are used for many triangles. Counted points on an edge or vertex are excluded.

Interface

template <ExactCoordinate T>
class CountPointsInTriangle {
public:
    CountPointsInTriangle(
        const std::vector<Point<T>>& triangle_vertices,
        const std::vector<Point<T>>& counted_points
    );

    int query(int first, int second, int third) const;
};

Member	Description	Complexity
`CountPointsInTriangle(triangle_vertices, counted_points)`	Precomputes the counts needed by every triangle. Neither input is mutated.	$O(N(N+M)\log(N+M))$ time and $O(N^2+N+M)$ memory
`query(first, second, third)`	Counts points strictly inside the triangle formed by the three indexed vertices.	$O(1)$

Here, $N$ is the number of triangle_vertices and $M$ is the number of counted_points.

Behavior and Requirements

T must satisfy ExactCoordinate<T>.
Query indices are zero-based and must refer to triangle_vertices.
The order of the three query indices does not matter.
A query with repeated or collinear triangle vertices returns 0.
Duplicate counted points are counted with multiplicity. Duplicate triangle vertices are also allowed.
The triangle boundary is excluded exactly; no floating-point tolerance is used.
Coordinate differences and products used by orientation tests must fit wide_type<T> (__int128_t for built-in integers, T for rationals).

Algorithm

For every candidate vertex, the constructor performs an angular sweep of the other candidate vertices and counted points above it. A Fenwick tree records how many counted points lie strictly to one side of each upward edge and how many lie on the edge. Counts on the same horizontal line are stored separately.

A query orders its three vertices by y-coordinate and expresses the open triangle as a signed combination of these precomputed edge regions. Separate strict and collinear counts ensure that every boundary point is excluded.

Example

#include "geometry/count_points_in_triangle.hpp"

#include <iostream>
#include <vector>

int main() {
    using Point = m1une::geometry::Point<long long>;

    std::vector<Point> vertices;
    vertices.emplace_back(0, 0);
    vertices.emplace_back(4, 0);
    vertices.emplace_back(0, 4);

    std::vector<Point> points;
    points.emplace_back(1, 1);
    points.emplace_back(2, 1);
    points.emplace_back(0, 1);  // Boundary: excluded.
    points.emplace_back(4, 4);
    points.emplace_back(1, 1);  // Duplicate: counted again.

    m1une::geometry::CountPointsInTriangle<long long> counter(
        vertices,
        points
    );
    std::cout << counter.query(0, 1, 2) << "\n";  // 3
}

Depends on

Required by

Geometry Bundle (geometry/all.hpp)

Verified with

Code

#ifndef M1UNE_GEOMETRY_COUNT_POINTS_IN_TRIANGLE_HPP
#define M1UNE_GEOMETRY_COUNT_POINTS_IN_TRIANGLE_HPP 1

#include <algorithm>
#include <cassert>
#include <concepts>
#include <cstddef>
#include <limits>
#include <numeric>
#include <utility>
#include <vector>

#include "point.hpp"

namespace m1une {
namespace geometry {

namespace count_points_in_triangle_detail {

class FenwickTree {
   private:
    std::vector<int> values;

   public:
    explicit FenwickTree(std::size_t size) : values(size, 0) {}

    void add(std::size_t position) {
        for (
            ++position;
            position <= values.size();
            position += position & -position
        ) {
            values[position - 1]++;
        }
    }

    int prefix_sum(std::size_t right) const {
        int result = 0;
        for (; right > 0; right -= right & -right) {
            result += values[right - 1];
        }
        return result;
    }
};

template <ExactCoordinate T>
struct SweepItem {
    using Wide = wide_type<T>;

    Wide x;
    Wide y;
    int index;
    int kind;
};

}  // namespace count_points_in_triangle_detail

template <ExactCoordinate T>
class CountPointsInTriangle {
   private:
    using Wide = wide_type<T>;
    using SweepItem = count_points_in_triangle_detail::SweepItem<T>;

    static constexpr int before_counted_points = 0;
    static constexpr int counted_point = 1;
    static constexpr int after_counted_points = 2;

    std::vector<Point<T>> vertices;
    std::vector<int> horizontal_left;
    std::vector<int> horizontal_equal;
    std::vector<std::vector<int>> edge_left;
    std::vector<std::vector<int>> edge_equal;

    static bool yx_less(const Point<T>& first, const Point<T>& second) {
        if (first.y != second.y) return first.y < second.y;
        return first.x < second.x;
    }

    SweepItem make_item(
        const Point<T>& origin,
        const Point<T>& point,
        int index,
        int kind
    ) const {
        return SweepItem{
            Wide(point.x) - Wide(origin.x),
            Wide(point.y) - Wide(origin.y),
            index,
            kind
        };
    }

    void build(const std::vector<Point<T>>& counted_points) {
        const int vertex_count = int(vertices.size());
        const int counted_point_count = int(counted_points.size());
        for (int anchor = 0; anchor < vertex_count; ++anchor) {
            for (const Point<T>& point : counted_points) {
                if (point.y != vertices[anchor].y) continue;
                if (point.x < vertices[anchor].x) {
                    horizontal_left[anchor]++;
                } else if (point.x == vertices[anchor].x) {
                    horizontal_equal[anchor]++;
                }
            }

            std::vector<SweepItem> items;
            items.reserve(2 * vertices.size() + counted_points.size());
            for (int index = 0; index < vertex_count; ++index) {
                if (vertices[anchor].y < vertices[index].y) {
                    items.push_back(make_item(
                        vertices[anchor],
                        vertices[index],
                        index,
                        before_counted_points
                    ));
                }
            }
            for (int index = 0; index < counted_point_count; ++index) {
                if (vertices[anchor].y < counted_points[index].y) {
                    items.push_back(make_item(
                        vertices[anchor],
                        counted_points[index],
                        index,
                        counted_point
                    ));
                }
            }
            for (int index = 0; index < vertex_count; ++index) {
                if (vertices[anchor].y < vertices[index].y) {
                    items.push_back(make_item(
                        vertices[anchor],
                        vertices[index],
                        index,
                        after_counted_points
                    ));
                }
            }

            std::sort(
                items.begin(),
                items.end(),
                [](const SweepItem& first, const SweepItem& second) {
                    const Wide determinant =
                        first.x * second.y - first.y * second.x;
                    if (determinant != 0) return determinant < 0;
                    return first.kind < second.kind;
                }
            );

            std::vector<std::size_t> height_order(items.size());
            std::iota(height_order.begin(), height_order.end(), 0);
            std::sort(
                height_order.begin(),
                height_order.end(),
                [&items](std::size_t first, std::size_t second) {
                    if (items[first].y != items[second].y) {
                        return items[first].y < items[second].y;
                    }
                    return items[first].kind % 2 < items[second].kind % 2;
                }
            );

            count_points_in_triangle_detail::FenwickTree fenwick(items.size());
            for (std::size_t position : height_order) {
                const SweepItem& item = items[position];
                if (item.kind == before_counted_points) {
                    edge_left[anchor][item.index] =
                        fenwick.prefix_sum(position + 1);
                } else if (item.kind == counted_point) {
                    fenwick.add(position);
                } else {
                    edge_equal[anchor][item.index] =
                        fenwick.prefix_sum(position + 1);
                }
            }
            for (int index = 0; index < vertex_count; ++index) {
                edge_equal[anchor][index] -= edge_left[anchor][index];
            }
        }
    }

   public:
    CountPointsInTriangle(
        const std::vector<Point<T>>& triangle_vertices,
        const std::vector<Point<T>>& counted_points
    )
        : vertices(triangle_vertices),
          horizontal_left(vertices.size(), 0),
          horizontal_equal(vertices.size(), 0),
          edge_left(vertices.size(), std::vector<int>(vertices.size(), 0)),
          edge_equal(vertices.size(), std::vector<int>(vertices.size(), 0)) {
        assert(
            vertices.size() <=
            static_cast<std::size_t>(std::numeric_limits<int>::max())
        );
        assert(
            counted_points.size() <=
            static_cast<std::size_t>(std::numeric_limits<int>::max())
        );
        build(counted_points);
    }

    int query(int first, int second, int third) const {
        const int size = int(vertices.size());
        assert(0 <= first && first < size);
        assert(0 <= second && second < size);
        assert(0 <= third && third < size);

        if (yx_less(vertices[second], vertices[first])) {
            std::swap(first, second);
        }
        if (yx_less(vertices[third], vertices[second])) {
            std::swap(second, third);
        }
        if (yx_less(vertices[second], vertices[first])) {
            std::swap(first, second);
        }

        const Wide determinant = cross(
            vertices[first],
            vertices[second],
            vertices[third]
        );
        if (determinant == 0) return 0;

        long long result;
        if (vertices[first].y == vertices[second].y) {
            result =
                static_cast<long long>(edge_left[second][third]) -
                edge_left[first][third] -
                edge_equal[first][third];
        } else if (vertices[second].y == vertices[third].y) {
            result =
                static_cast<long long>(edge_left[first][third]) -
                edge_left[first][second] -
                edge_equal[first][second];
        } else if (determinant < 0) {
            result =
                static_cast<long long>(edge_left[first][third]) -
                edge_left[second][third] -
                edge_equal[second][third] -
                edge_left[first][second] -
                edge_equal[first][second] -
                horizontal_left[second] -
                horizontal_equal[second];
        } else {
            result =
                static_cast<long long>(edge_left[first][second]) +
                edge_left[second][third] +
                horizontal_left[second] -
                edge_left[first][third] -
                edge_equal[first][third];
        }
        assert(
            0 <= result &&
            result <= std::numeric_limits<int>::max()
        );
        return int(result);
    }
};

}  // namespace geometry
}  // namespace m1une

#endif  // M1UNE_GEOMETRY_COUNT_POINTS_IN_TRIANGLE_HPP

#line 1 "geometry/count_points_in_triangle.hpp"



#include <algorithm>
#include <cassert>
#include <concepts>
#include <cstddef>
#include <limits>
#include <numeric>
#include <utility>
#include <vector>

#line 1 "geometry/point.hpp"



#include <cmath>
#line 7 "geometry/point.hpp"
#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 14 "geometry/count_points_in_triangle.hpp"

namespace m1une {
namespace geometry {

namespace count_points_in_triangle_detail {

class FenwickTree {
   private:
    std::vector<int> values;

   public:
    explicit FenwickTree(std::size_t size) : values(size, 0) {}

    void add(std::size_t position) {
        for (
            ++position;
            position <= values.size();
            position += position & -position
        ) {
            values[position - 1]++;
        }
    }

    int prefix_sum(std::size_t right) const {
        int result = 0;
        for (; right > 0; right -= right & -right) {
            result += values[right - 1];
        }
        return result;
    }
};

template <ExactCoordinate T>
struct SweepItem {
    using Wide = wide_type<T>;

    Wide x;
    Wide y;
    int index;
    int kind;
};

}  // namespace count_points_in_triangle_detail

template <ExactCoordinate T>
class CountPointsInTriangle {
   private:
    using Wide = wide_type<T>;
    using SweepItem = count_points_in_triangle_detail::SweepItem<T>;

    static constexpr int before_counted_points = 0;
    static constexpr int counted_point = 1;
    static constexpr int after_counted_points = 2;

    std::vector<Point<T>> vertices;
    std::vector<int> horizontal_left;
    std::vector<int> horizontal_equal;
    std::vector<std::vector<int>> edge_left;
    std::vector<std::vector<int>> edge_equal;

    static bool yx_less(const Point<T>& first, const Point<T>& second) {
        if (first.y != second.y) return first.y < second.y;
        return first.x < second.x;
    }

    SweepItem make_item(
        const Point<T>& origin,
        const Point<T>& point,
        int index,
        int kind
    ) const {
        return SweepItem{
            Wide(point.x) - Wide(origin.x),
            Wide(point.y) - Wide(origin.y),
            index,
            kind
        };
    }

    void build(const std::vector<Point<T>>& counted_points) {
        const int vertex_count = int(vertices.size());
        const int counted_point_count = int(counted_points.size());
        for (int anchor = 0; anchor < vertex_count; ++anchor) {
            for (const Point<T>& point : counted_points) {
                if (point.y != vertices[anchor].y) continue;
                if (point.x < vertices[anchor].x) {
                    horizontal_left[anchor]++;
                } else if (point.x == vertices[anchor].x) {
                    horizontal_equal[anchor]++;
                }
            }

            std::vector<SweepItem> items;
            items.reserve(2 * vertices.size() + counted_points.size());
            for (int index = 0; index < vertex_count; ++index) {
                if (vertices[anchor].y < vertices[index].y) {
                    items.push_back(make_item(
                        vertices[anchor],
                        vertices[index],
                        index,
                        before_counted_points
                    ));
                }
            }
            for (int index = 0; index < counted_point_count; ++index) {
                if (vertices[anchor].y < counted_points[index].y) {
                    items.push_back(make_item(
                        vertices[anchor],
                        counted_points[index],
                        index,
                        counted_point
                    ));
                }
            }
            for (int index = 0; index < vertex_count; ++index) {
                if (vertices[anchor].y < vertices[index].y) {
                    items.push_back(make_item(
                        vertices[anchor],
                        vertices[index],
                        index,
                        after_counted_points
                    ));
                }
            }

            std::sort(
                items.begin(),
                items.end(),
                [](const SweepItem& first, const SweepItem& second) {
                    const Wide determinant =
                        first.x * second.y - first.y * second.x;
                    if (determinant != 0) return determinant < 0;
                    return first.kind < second.kind;
                }
            );

            std::vector<std::size_t> height_order(items.size());
            std::iota(height_order.begin(), height_order.end(), 0);
            std::sort(
                height_order.begin(),
                height_order.end(),
                [&items](std::size_t first, std::size_t second) {
                    if (items[first].y != items[second].y) {
                        return items[first].y < items[second].y;
                    }
                    return items[first].kind % 2 < items[second].kind % 2;
                }
            );

            count_points_in_triangle_detail::FenwickTree fenwick(items.size());
            for (std::size_t position : height_order) {
                const SweepItem& item = items[position];
                if (item.kind == before_counted_points) {
                    edge_left[anchor][item.index] =
                        fenwick.prefix_sum(position + 1);
                } else if (item.kind == counted_point) {
                    fenwick.add(position);
                } else {
                    edge_equal[anchor][item.index] =
                        fenwick.prefix_sum(position + 1);
                }
            }
            for (int index = 0; index < vertex_count; ++index) {
                edge_equal[anchor][index] -= edge_left[anchor][index];
            }
        }
    }

   public:
    CountPointsInTriangle(
        const std::vector<Point<T>>& triangle_vertices,
        const std::vector<Point<T>>& counted_points
    )
        : vertices(triangle_vertices),
          horizontal_left(vertices.size(), 0),
          horizontal_equal(vertices.size(), 0),
          edge_left(vertices.size(), std::vector<int>(vertices.size(), 0)),
          edge_equal(vertices.size(), std::vector<int>(vertices.size(), 0)) {
        assert(
            vertices.size() <=
            static_cast<std::size_t>(std::numeric_limits<int>::max())
        );
        assert(
            counted_points.size() <=
            static_cast<std::size_t>(std::numeric_limits<int>::max())
        );
        build(counted_points);
    }

    int query(int first, int second, int third) const {
        const int size = int(vertices.size());
        assert(0 <= first && first < size);
        assert(0 <= second && second < size);
        assert(0 <= third && third < size);

        if (yx_less(vertices[second], vertices[first])) {
            std::swap(first, second);
        }
        if (yx_less(vertices[third], vertices[second])) {
            std::swap(second, third);
        }
        if (yx_less(vertices[second], vertices[first])) {
            std::swap(first, second);
        }

        const Wide determinant = cross(
            vertices[first],
            vertices[second],
            vertices[third]
        );
        if (determinant == 0) return 0;

        long long result;
        if (vertices[first].y == vertices[second].y) {
            result =
                static_cast<long long>(edge_left[second][third]) -
                edge_left[first][third] -
                edge_equal[first][third];
        } else if (vertices[second].y == vertices[third].y) {
            result =
                static_cast<long long>(edge_left[first][third]) -
                edge_left[first][second] -
                edge_equal[first][second];
        } else if (determinant < 0) {
            result =
                static_cast<long long>(edge_left[first][third]) -
                edge_left[second][third] -
                edge_equal[second][third] -
                edge_left[first][second] -
                edge_equal[first][second] -
                horizontal_left[second] -
                horizontal_equal[second];
        } else {
            result =
                static_cast<long long>(edge_left[first][second]) +
                edge_left[second][third] +
                horizontal_left[second] -
                edge_left[first][third] -
                edge_equal[first][third];
        }
        assert(
            0 <= result &&
            result <= std::numeric_limits<int>::max()
        );
        return int(result);
    }
};

}  // namespace geometry
}  // namespace m1une

m1une's library

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

Count Points in Triangle
(geometry/count_points_in_triangle.hpp)

Overview

Interface

Behavior and Requirements

Algorithm

Example

Depends on

Required by

Verified with

Code

Count Points in Triangle (geometry/count_points_in_triangle.hpp)

Overview

Interface

Behavior and Requirements

Algorithm

Example

Depends on

Required by

Verified with

Code

Count Points in Triangle
(geometry/count_points_in_triangle.hpp)