Split points into circles

CMakeLists.txt

cmake_minimum_required (VERSION 3.8)

project(split_into_circles)
set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Werror -pedantic -DTEST")

add_executable(run_test split_into_circles.cpp)

split_into_circles.cpp

#include <iostream>
#include <random>
#include <vector>
#include <cassert>
#include <optional>
#include <algorithm>
#include <random>
#include <map>
#include <chrono>


using Point = std::pair<double, double>;

namespace details {
    using std::vector;

    using Circle = std::pair<Point, double>;
    constexpr double EPS = 1e-7;

    Point minus(const Point &p1, const Point &p2) {
        return {p1.first - p2.first, p1.second - p2.second};
    }

    double norm2(const Point &v1) {
        return v1.first * v1.first + v1.second * v1.second;
    }

    std::optional<Circle> get_circle(
            const Point &p1,
            const Point &p2,
            const Point &p3
    ) {
        double d = 2 * (p1.first * (p2.second - p3.second)
                        - p1.second * (p2.first - p3.first)
                        + p2.first * p3.second - p3.first * p2.second);

        if (std::abs(d) < EPS) {
            return std::nullopt;
        }

        double x = norm2(p1) * (p2.second - p3.second)
                   + norm2(p2) * (p3.second - p1.second)
                   + norm2(p3) * (p1.second - p2.second);

        double y = norm2(p1) * (p3.first - p2.first)
                   + norm2(p2) * (p1.first - p3.first)
                   + norm2(p3) * (p2.first - p1.first);

        Point center{x / d, y / d};
        double radius2 = norm2(minus(center, p1));

        return Circle{center, radius2};
    }

    bool is_point_on_circle(const Point &p, const Circle &circle) {
        return std::abs(norm2(minus(p, circle.first)) - circle.second) < EPS;
    }

    bool split_into_circles_(
            const vector<Point> &points,
            const vector<size_t> &idx,
            vector<size_t> &result,
            size_t n
    ) {
        if (idx.empty()) {
            return true;
        }
        if (n == 0) {
            return false;
        }

        size_t N = 2 * n + 1;

        vector<Point> sel;

        for (auto point_idx : idx) {
            sel.push_back(points[point_idx]);
            if (sel.size() >= N) break;
        }

        for (size_t i = 0; i < N; ++i) {
            for (size_t j = i + 1; j < N; ++j) {
                for (size_t k = j + 1; k < N; ++k) {
                    auto circle_opt = get_circle(sel[i], sel[j], sel[k]);

                    if (circle_opt.has_value()) {
                        const auto &circle = circle_opt.value();
                        vector<size_t> idx2;

                        for (size_t point_idx : idx) {
                            if (is_point_on_circle(points[point_idx], circle)) {
                                result[point_idx] = n;
                            } else {
                                idx2.push_back(point_idx);
                            }
                        }

                        if (split_into_circles_(points, idx2, result, n - 1)) {
                            return true;
                        }
                    }
                }
            }
        }

        return false;
    }

}

std::vector<size_t> split_into_circles(const std::vector<Point> &points, size_t n_circles) {
    std::vector<size_t> idx(points.size()), result(points.size());

    for (size_t i = 0; i < idx.size(); ++i) {
        idx[i] = i;
    }

    if (!(details::split_into_circles_(points, idx, result, n_circles))) {
        result.clear();
    }

    return result;
}


void test() {
    // RANDOM TEST
    details::Circle circles[3];

    std::uniform_real_distribution<double> unif(-1000, 1000);
    std::default_random_engine re;
    auto rand_double = [&unif, &re] { return unif(re); };

    for (auto &circle : circles) {
        double x = rand_double()
        , y = rand_double()
        , r = rand_double();
        circle = {{x, y}, r};
    }

    std::vector<Point> points;
    std::vector<size_t> expected, indexes;

    for (size_t i = 0; i < 3; ++i) {
        size_t n = 20000;
        for (size_t _ = 0; _ < n; ++_) {
            expected.push_back(i);
            indexes.push_back(indexes.size());
            double angle = rand_double() * M_PI / 1000;
            double x = circles[i].first.first + circles[i].second * cos(angle);
            double y = circles[i].first.second + circles[i].second * sin(angle);
            points.emplace_back(x, y);
        }
    }

    std::shuffle(indexes.begin(), indexes.end(), std::mt19937(std::random_device()()));
    std::vector<Point> points_shuffled(points.size());
    std::vector<size_t> expected_shuffled(expected.size());

    for (size_t i = 0; i < indexes.size(); ++i) {
        points_shuffled[indexes[i]] = points[i];
        expected_shuffled[indexes[i]] = expected[i];
    }

    auto start = std::chrono::high_resolution_clock::now();
    auto result = split_into_circles(points_shuffled, 3);
    std::chrono::duration<double> dtime = std::chrono::high_resolution_clock::now() - start;

    std::map<size_t, size_t> map_index;
    for (size_t i = 0; i < result.size(); ++i) {
        if (map_index.find(expected_shuffled[i]) == map_index.end()) {
            map_index[expected_shuffled[i]] = result[i];
        } else {
            assert(result[i] == map_index[expected_shuffled[i]]);
        }
    }

    std::cout << "success at " << points.size() << " points!\n";
    std::cout << "time: " << dtime.count() << "s\n";
}

int main() {
#ifdef TEST
    test();
    return 0;
#else
    size_t n, m;
    std::cin >> n >> m;
    std::vector<Point> points(n);
    for (auto & point : points) {
        double x, y;
        std::cin >> x >> y;
        point = {x, y};
    }

    for (auto idx : split_into_circles(points, m)) {
        std::cout << idx << '\n';
    }

    return 0;
#endif
}

Линейная сложность при фиксированном M