cp-documentation

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

View the Project on GitHub zawa-tin/cp-documentation

:heavy_check_mark: Test/My/GeometryZ2/Contain/CountingPointsInTrianglesStressTest.test.cpp

Depends on

Code

#define PROBLEM "https://onlinejudge.u-aizu.ac.jp/courses/lesson/2/ITP1/1/ITP1_1_A"

#include "../../../../Src/GeometryZ2/Contain/CountPointsInTriangles.hpp"
#include "../../../../Src/GeometryZ2/Contain/NaiveCountPointsInTriangles.hpp"

#include <iostream>
#include <vector>
#include <random>
using namespace zawa;
using namespace geometryZ2;

std::mt19937 mt{std::random_device{}()};

bool test(const int N, const int M, const int X, const int Y) {
    const int n = mt() % (N - 2) + 3;
    const int m = mt() % (M - 2) + 3;
    std::vector<Point> a(n), b(m);
    for (int i = 0 ; i < n ; i++) {
        a[i].x() = mt() % (X + 1);
        a[i].y() = mt() % (Y + 1);
    }
    for (int i = 0 ; i < m ; i++) {
        b[i].x() = mt() % (X + 1);
        b[i].y() = mt() % (Y + 1);
    }
    NaiveCountPointsInTriangles ans(a, b);
    CountPointsInTriangles my(a, b);
    auto query = [&](int i, int j, int k) -> bool {
        int p = ans(i, j, k), q = my(i, j, k);
        if (p != q) {
            std::cout << n << '\n';
            for (const Point& v : a) std::cout << v.x() << ' ' << v.y() << '\n';
            std::cout << m << '\n';
            for (const Point& v : b) std::cout << v.x() << ' ' << v.y() << '\n';
            std::cout << 1 << '\n';
            std::cout << i << ' ' << j << ' ' << k << std::endl;
            std::cerr << p << " but you output " << q << std::endl;
            return false;
        }
        return true;
    };
    for (int i = 0 ; i < std::ssize(a) ; i++) {
        for (int j = 0 ; j < std::ssize(a) ; j++) {
            for (int k = 0 ; k < std::ssize(a) ; k++) {
                if (!query(i, j, k)) return false;
            }
        }
    }
    return true;
}

int main() {
    for (int t = 0 ; t < 1000 ; t++) { // tiny
        assert(test(5, 5, 5, 5));
    }
    for (int t = 0 ; t < 100 ; t++) { // small dense
        assert(test(20, 20, 10, 10));
    }
    for (int t = 0 ; t < 100 ; t++) { // small
        assert(test(20, 20, 100, 100));
    }
    assert(test(50, 100, 100, 100));
    std::cout << "Hello World\n";
}
#line 1 "Test/My/GeometryZ2/Contain/CountingPointsInTrianglesStressTest.test.cpp"
#define PROBLEM "https://onlinejudge.u-aizu.ac.jp/courses/lesson/2/ITP1/1/ITP1_1_A"

#line 2 "Src/GeometryZ2/Contain/CountPointsInTriangles.hpp"

#line 2 "Src/GeometryZ2/PointCloud.hpp"

#line 2 "Src/GeometryZ2/Point.hpp"

#line 2 "Src/Template/TypeAlias.hpp"

#include <cstdint>
#include <cstddef>

namespace zawa {

using i16 = std::int16_t;
using i32 = std::int32_t;
using i64 = std::int64_t;
using i128 = __int128_t;

using u8 = std::uint8_t;
using u16 = std::uint16_t;
using u32 = std::uint32_t;
using u64 = std::uint64_t;

using usize = std::size_t;

} // namespace zawa
#line 2 "Src/GeometryZ2/Zahlen.hpp"

#line 4 "Src/GeometryZ2/Zahlen.hpp"

#include <cassert>

namespace zawa {

namespace geometryZ2 {

using Zahlen = i64;

namespace internal {

constexpr i32 positive{1};
constexpr i32 zero{0};
constexpr i32 negative{-1};

} // namespace internal

constexpr i32 Sign(Zahlen value) {
    if (value < 0) return internal::negative;
    if (value > 0) return internal::positive;
    return internal::zero;
}

constexpr bool Positive(Zahlen value) {
    return Sign(value) == internal::positive;
}

constexpr bool Zero(Zahlen value) {
    return Sign(value) == internal::zero;
}

constexpr bool Negative(Zahlen value) {
    return Sign(value) == internal::negative;
}

constexpr Zahlen Abs(Zahlen value) {
    return (value > 0 ? value : -value);
}

constexpr Zahlen Square(Zahlen value) {
    return value * value;
}

} // namespace geometryZ2

} // namespace zawa
#line 5 "Src/GeometryZ2/Point.hpp"

#include <algorithm>
#include <iostream>
#line 9 "Src/GeometryZ2/Point.hpp"
#include <limits>

namespace zawa {

namespace geometryZ2 {

class Point {
private:
    Zahlen x_{}, y_{};
    static constexpr i32 origin{0};
    static constexpr i32 firstQuadrant{1};
    static constexpr i32 secondQuadrant{2};
    static constexpr i32 thirdQuadrant{-2};
    static constexpr i32 forthQuadrant{-1};
public:
    /* constructor */
    Point() = default;
    Point(const Point& p) : x_{p.x()}, y_{p.y()} {}
    Point(Zahlen x, Zahlen y) : x_{x}, y_{y} {}

    /* getter setter */
    Zahlen& x() {
        return x_;
    }
    const Zahlen& x() const {
        return x_;
    }
    Zahlen& y() {
        return y_;
    }
    const Zahlen& y() const {
        return y_;
    }

    /* operator */
    Point& operator=(const Point& p) {
        x() = p.x();
        y() = p.y();
        return *this;
    }
    Point& operator+=(const Point& p) {
        x() += p.x();
        y() += p.y();
        return *this;
    }
    friend Point operator+(const Point& p0, const Point& p1) {
        return Point{p0} += p1;
    }
    Point& operator-=(const Point& p) {
        x() -= p.x();
        y() -= p.y();
        return *this;
    }
    friend Point operator-(const Point& p0, const Point& p1) {
        return Point{p0} -= p1;
    }
    Point& operator*=(Zahlen k) {
        x() *= k;
        y() *= k;
        return *this;
    }
    friend Point operator*(const Point& p, Zahlen k) {
        return Point{p} *= k;
    }
    friend Point operator*(Zahlen k, const Point& p) {
        return Point{p} *= k;
    }
    Point& operator/=(Zahlen k) {
        assert(k);
        assert(x() % k == 0);
        assert(y() % k == 0);
        x() /= k;
        y() /= k;
        return *this;
    }
    friend Point operator/(const Point& p, Zahlen k) {
        return Point{p} /= k;
    }
    friend bool operator==(const Point& p0, const Point& p1) {
        return p0.x() == p1.x() and p0.y() == p1.y();
    }
    friend bool operator!=(const Point& p0, const Point& p1) {
        return p0.x() != p1.x() or p0.y() != p1.y();
    }
    friend bool operator<(const Point& p0, const Point& p1) {
        if (p0.x() != p1.x()) return p0.x() < p1.x();
        else return p0.y() < p1.y();
    }
    friend bool operator<=(const Point& p0, const Point& p1) {
        return (p0 < p1) or (p0 == p1);
    }
    friend bool operator>(const Point& p0, const Point& p1) {
        if (p0.x() != p1.x()) return p0.x() > p1.x();
        else return p0.y() > p1.y();
    }
    friend bool operator>=(const Point& p0, const Point& p1) {
        return (p0 > p1) or (p0 == p1);
    }
    friend std::istream& operator>>(std::istream& is, Point& p) {
        is >> p.x() >> p.y();
        return is;
    }
    friend std::ostream& operator<<(std::ostream& os, const Point& p) {
        os << '(' << p.x() << ',' << p.y() << ')';
        return os;
    }

    /* member function */
    Zahlen normSquare() const {
        return Square(x()) + Square(y());
    }
    bool isNormSquareOver(Zahlen d) const {
        assert(!Negative(d));
        auto [mn, mx]{std::minmax({ Abs(x()), Abs(y()) })};
        if (mx and mx > d / mx) {
            return true;
        }
        long long s1{Square(mn)}, s2{Square(mx)};
        if (s1 > d - s2) {
            return true;
        }
        return false;
    }
    bool isNormSquareOverflow() const {
        return isNormSquareOver(std::numeric_limits<Zahlen>::max());
    }

    i32 area() const {
        if (x_ == 0 and y_ == 0) return origin;
        if (x_ <= 0 and y_ < 0) return thirdQuadrant;
        if (x_ > 0 and y_ <= 0) return forthQuadrant;
        if (x_ >= 0 and y_ > 0) return firstQuadrant;
        return secondQuadrant;
    }

    /* static member */
    static bool ArgComp(const Point& p0, const Point& p1) {
        if (p0.area() != p1.area()) return p0.area() < p1.area();
        Zahlen cross{Cross(p0, p1)};
        return (!Zero(cross) ? Positive(cross) : p0.normSquare() < p1.normSquare());
    }

    /* friend function */
    friend Zahlen Dot(const Point& p0, const Point& p1) {
        return p0.x() * p1.x() + p0.y() * p1.y();
    }
    friend Zahlen Cross(const Point& p0, const Point& p1) {
        return p0.x() * p1.y() - p0.y() * p1.x();
    }
};
using Vector = Point;

} // namespace geometryZ2

} // namespace zawa
#line 4 "Src/GeometryZ2/PointCloud.hpp"

#line 6 "Src/GeometryZ2/PointCloud.hpp"
#include <vector>

namespace zawa {

namespace geometryZ2 {

using PointCloud = std::vector<Point>;

void ArgSort(PointCloud& p) {
    std::sort(p.begin(), p.end(), Point::ArgComp);
}

} // namespace geometryZ2 

} // namespace zawa
#line 4 "Src/GeometryZ2/Contain/CountPointsInTriangles.hpp"

#line 6 "Src/GeometryZ2/Contain/CountPointsInTriangles.hpp"
#include <numeric>
#include <utility>
#line 9 "Src/GeometryZ2/Contain/CountPointsInTriangles.hpp"

namespace zawa {

namespace geometryZ2 {

class CountPointsInTriangles {
public:

    CountPointsInTriangles(PointCloud a, PointCloud b) 
        : m_a(a.size()), m_under(a.size()), m_eq(a.size()), m_on(a.size()), m_cover(a.size()), m_inv(a.size()) {
        {
            std::vector<std::pair<Point, usize>> idxs(m_a.size());
            for (usize i = 0 ; i < m_a.size() ; i++) {
                m_on[i].resize(m_a.size() - i);
                m_cover[i].resize(m_a.size() - i);
                idxs[i].first = std::move(a[i]);
                idxs[i].second = i;
            }
            std::sort(idxs.begin(), idxs.end());
            for (usize i = 0 ; i < m_a.size() ; i++) {
                m_a[i] = std::move(idxs[i].first);
                m_inv[idxs[i].second] = i;
            }
        }
        std::sort(b.begin(), b.end());
        // calc m_under and m_eq
        for (usize i = 0 ; i < m_a.size() ; i++) {
            for (usize j = 0 ; j < b.size() and b[j].x() <= m_a[i].x() ; j++) {
                if (b[j].x() == m_a[i].x()) {
                    if (b[j].y() < m_a[i].y()) m_under[i]++;
                    else if (b[j].y() == m_a[i].y()) m_eq[i]++;
                }
            } 
        }
        // calc m_cover, m_on
        for (usize i = 0 ; i < m_a.size() ; i++) {
            std::vector<std::pair<Point, usize>> dj;
            std::vector<Vector> dirs;
            {
                std::vector<std::pair<Point, usize>> da, db;
                for (usize j = i + 1 ; j < m_a.size() ; j++) if (m_a[i] != m_a[j]) {
                    da.push_back({m_a[j]-m_a[i], j});
                }
                for (const auto& q : b) if (m_a[i].x() <= q.x() and m_a[i] != q) {
                    db.emplace_back(q-m_a[i],std::numeric_limits<usize>::max());
                    dirs.push_back(q-m_a[i]);
                }
                std::ranges::sort(db);
                dj.resize(db.size() + da.size());
                for (usize j = 0, k = 0, t = 0 ; t < dj.size() ; t++) {
                    if (k == db.size()) dj[t] = std::move(da[j++]);
                    else if (j == da.size()) dj[t] = std::move(db[k++]);
                    else if (da[j] < db[k]) dj[t] = std::move(da[j++]);
                    else dj[t] = std::move(db[k++]);
                }
            }
            std::ranges::sort(dirs, Point::ArgComp);
            dirs.erase(
                    std::unique(dirs.begin(), dirs.end(), [&](const auto& l, const auto& r) { return Zero(Cross(l, r)) and !Negative(Dot(l, r)); }), 
                    dirs.end());
            std::vector<u32> fen(dirs.size() + 1), dat(dirs.size());
            for (const auto& [d, j] : dj) {
                if (j == std::numeric_limits<usize>::max()) {
                    auto it = std::distance(
                            dirs.begin(), 
                            std::ranges::lower_bound(dirs, d, [](const Point& l, const Point& r) {
                                    if (l.area() != r.area()) return l.area() < r.area();
                                    return Positive(Cross(l, r));
                                })
                            );
                    dat[it]++;
                    for (it++ ; it < std::ssize(fen) ; it += it & -it) fen[it]++;
                }
                else {
                    auto it = std::distance(dirs.begin(), std::ranges::upper_bound(dirs, d, Point::ArgComp));
                    if (it and Zero(Cross(d, dirs[it - 1])) and !Negative(Dot(d, dirs[it - 1]))) {
                        m_on[i][j - i] = dat[it - 1];
                    }
                    for ( ; it ; it -= it & -it) m_cover[i][j - i] += fen[it];
                }
            }
        }
    }

    inline usize size() const {
        return m_a.size();
    }

    // strictly inside
    u32 operator()(usize i, usize j, usize k) const {
        assert(i < size() and j < size() and k < size());
        if (i == j or j == k or i == k) return 0;
        i = m_inv[i];
        j = m_inv[j];
        k = m_inv[k];
        if (m_a[i] == m_a[j] or m_a[j] == m_a[k] or m_a[k] == m_a[i]) return 0;
        if (m_a[i] > m_a[j]) std::swap(i, j);
        if (m_a[j] > m_a[k]) std::swap(j, k);
        if (m_a[i] > m_a[j]) std::swap(i, j);
        // std::cout << i << ' ' << j << ' ' << k << std::endl;
        const Zahlen crs = Cross(m_a[j] - m_a[i], m_a[k] - m_a[i]);
        if (crs == 0) {
            return 0;
        }
        else if (crs > 0) {
            return cover(i, k) - cover(i, j) - cover(j, k) - m_eq[j] + m_under[j] - on(i, k);
        }
        else {
            return cover(i, j) + cover(j, k) - cover(i, k) - m_under[j] - on(i, j) - on(j, k);
        }
    }
    
private:

    std::vector<Point> m_a;

    std::vector<u32> m_under, m_eq;

    std::vector<std::vector<u32>> m_on, m_cover;
    
    std::vector<usize> m_inv;

    u32 cover(usize i, usize j) const {
        if (i > j) std::swap(i, j);
        return m_cover[i][j - i];
    }

    u32 on(usize i, usize j) const {
        if (i > j) std::swap(i, j);
        return m_on[i][j - i];
    }
};

} // namespace geometryZ2

} // namespace zawa
#line 2 "Src/GeometryZ2/Contain/NaiveCountPointsInTriangles.hpp"

#line 2 "Src/GeometryZ2/Relation.hpp"

#line 5 "Src/GeometryZ2/Relation.hpp"

namespace zawa {

namespace geometryZ2 {

enum RELATION {
    // p0 -> p1 -> p2の順で直線上に並んでいる
    ONLINE_FRONT        = -2,
    // (p1 - p0) -> (p2 - p0)が時計回りになっている
    CLOCKWISE           = -1,
    // p0 -> p2 -> p1の順で直線上に並んでいる
    ON_SEGMENT          =  0,
    // (p1 - p0) -> (p2 - p0)が反時計回りになっている
    COUNTER_CLOCKWISE   = +1,
    // p2 -> p0 -> p1、またはp1 -> p0 -> p2の順で直線上に並んでいる
    ONLINE_BACK         = +2
};

RELATION Relation(const Point& p0, const Point& p1, const Point& p2) {
    Point a{p1 - p0}, b{p2 - p0};
    if (Positive(Cross(a, b))) return COUNTER_CLOCKWISE;
    if (Negative(Cross(a, b))) return CLOCKWISE;
    if (Negative(Dot(a, b))) return ONLINE_BACK;
    if (a.normSquare() < b.normSquare()) return ONLINE_FRONT;
    return ON_SEGMENT;
};

} // namespace geometryZ2

} // namespace zawa
#line 5 "Src/GeometryZ2/Contain/NaiveCountPointsInTriangles.hpp"

namespace zawa {

namespace geometryZ2 {

class NaiveCountPointsInTriangles {
public:

    NaiveCountPointsInTriangles(PointCloud a, PointCloud b) : m_a{std::move(a)}, m_b{std::move(b)} {}

    u32 operator()(u32 p, u32 q, u32 r) const {
        assert(p < size() and q < size() and r < size());
        if (m_a[p] > m_a[q]) std::swap(p, q);
        if (m_a[q] > m_a[r]) std::swap(q, r);
        if (m_a[p] > m_a[q]) std::swap(p, q);
        RELATION R = Relation(m_a[p], m_a[q], m_a[r]);
        if (R == RELATION::ONLINE_FRONT or R == RELATION::ONLINE_BACK or R == RELATION::ON_SEGMENT) return 0;
        bool ctr = R == RELATION::COUNTER_CLOCKWISE;
        u32 res = 0;
        for (const Point& i : m_b) {
            Zahlen a = Cross(m_a[q] - m_a[p], i - m_a[p]), b = Cross(m_a[r] - m_a[q], i - m_a[q]), c = Cross(m_a[p] - m_a[r], i - m_a[r]);
            if (ctr and a > 0 and b > 0 and c > 0) res++;
            else if (!ctr and a < 0 and b < 0 and c < 0) res++;
        }
        return res;
    }

    inline usize size() const {
        return m_a.size();
    }

private:

    std::vector<Point> m_a, m_b;
};

} // namespace geometryZ2

} // namespace zawa
#line 5 "Test/My/GeometryZ2/Contain/CountingPointsInTrianglesStressTest.test.cpp"

#line 8 "Test/My/GeometryZ2/Contain/CountingPointsInTrianglesStressTest.test.cpp"
#include <random>
using namespace zawa;
using namespace geometryZ2;

std::mt19937 mt{std::random_device{}()};

bool test(const int N, const int M, const int X, const int Y) {
    const int n = mt() % (N - 2) + 3;
    const int m = mt() % (M - 2) + 3;
    std::vector<Point> a(n), b(m);
    for (int i = 0 ; i < n ; i++) {
        a[i].x() = mt() % (X + 1);
        a[i].y() = mt() % (Y + 1);
    }
    for (int i = 0 ; i < m ; i++) {
        b[i].x() = mt() % (X + 1);
        b[i].y() = mt() % (Y + 1);
    }
    NaiveCountPointsInTriangles ans(a, b);
    CountPointsInTriangles my(a, b);
    auto query = [&](int i, int j, int k) -> bool {
        int p = ans(i, j, k), q = my(i, j, k);
        if (p != q) {
            std::cout << n << '\n';
            for (const Point& v : a) std::cout << v.x() << ' ' << v.y() << '\n';
            std::cout << m << '\n';
            for (const Point& v : b) std::cout << v.x() << ' ' << v.y() << '\n';
            std::cout << 1 << '\n';
            std::cout << i << ' ' << j << ' ' << k << std::endl;
            std::cerr << p << " but you output " << q << std::endl;
            return false;
        }
        return true;
    };
    for (int i = 0 ; i < std::ssize(a) ; i++) {
        for (int j = 0 ; j < std::ssize(a) ; j++) {
            for (int k = 0 ; k < std::ssize(a) ; k++) {
                if (!query(i, j, k)) return false;
            }
        }
    }
    return true;
}

int main() {
    for (int t = 0 ; t < 1000 ; t++) { // tiny
        assert(test(5, 5, 5, 5));
    }
    for (int t = 0 ; t < 100 ; t++) { // small dense
        assert(test(20, 20, 10, 10));
    }
    for (int t = 0 ; t < 100 ; t++) { // small
        assert(test(20, 20, 100, 100));
    }
    assert(test(50, 100, 100, 100));
    std::cout << "Hello World\n";
}
Back to top page