Generalization to arbitrary dimensions of https://www.reddit.com/r/dailyprogrammer/comments/6ksmh5/20170630_challenge_321_hard_circle_splitter/ Uses Eigen library for linear algebra http://eigen.tuxfamily.org/

smallest_nsphere.cpp

#include <iostream>
#include <ios>
#include <vector>
#include <algorithm>
#include <cmath>
#include <utility>
#include <tuple>

#define EIGEN_DONT_VECTORIZE
#define EIGEN_FAST_MATH 0

#include "Eigen/Dense"

using std::cout;
using std::cin;
using std::endl;
using std::fixed;

using std::vector;
using std::pair;
using std::tuple;

typedef unsigned int uint;

const double epsilon = 1e-10;

void print_vector(Eigen::VectorXd vector)
{
    for(uint i = 0; i < vector.size(); i++)
    {
        cout << vector(i) << " ";
    }
    
    cout << endl;
}

/*
 * return n-dimensional unit vector with first value one, and the rest zero
 */
Eigen::VectorXd base_vector(uint n_dimensions)
{
    Eigen::VectorXd base(n_dimensions);
    base.setZero();
    base(0) = 1;

    return base;
}

/*
 * derive a rotation matrix that maps one vector to another
 */
Eigen::MatrixXd rotation_matrix(const Eigen::VectorXd& a, const Eigen::VectorXd& b)
{
    uint size = a.rows();

    Eigen::VectorXd u = a.normalized();
    Eigen::VectorXd v = (b - u.dot(b) * u).normalized();

    double x = a.dot(b) / (a.norm() * b.norm());
    double y = std::sqrt(1 - std::pow(x, 2));

    Eigen::MatrixXd uu = u * u.transpose();
    Eigen::MatrixXd vv = v * v.transpose();
    Eigen::MatrixXd uv = (Eigen::MatrixXd(size, 2) << u, v).finished();
    Eigen::MatrixXd xy = (Eigen::MatrixXd(2, 2) << x, -y, y, x).finished();

    return Eigen::MatrixXd::Identity(size, size) - (uu + vv) + uv * xy * uv.transpose();
}

/*
 * derive a unit vector which is orthogonal to all provided vectors
 */
Eigen::VectorXd orthogonal_vector(const vector<Eigen::VectorXd>& vectors)
{
    uint n_dimensions = vectors[0].rows();

    Eigen::MatrixXd m(n_dimensions, n_dimensions);
    m.setZero();

    for(uint i = 0; i < vectors.size(); i++)
    {
        m.row(i + 1) = vectors[i];
    }

    /*
     * add in some base vectors if there aren't enough already
     * make sure they're not equivalent to any other vectors
     */
    if(vectors.size() + 1 < n_dimensions)
    {
        vector<Eigen::VectorXd> base_vectors(n_dimensions, Eigen::VectorXd(n_dimensions));

        for(uint i = n_dimensions - 1; i < n_dimensions; i--)
        {
            base_vectors[i].setZero();
            base_vectors[i](i) = 1;

            for(const Eigen::VectorXd& given : vectors)
            {
                if(base_vectors[i].isApprox(given.cwiseAbs().normalized(), epsilon))
                {
                    base_vectors.erase(base_vectors.begin() + i);
                }
            }
        }

        for(uint i = vectors.size() + 1; i < n_dimensions; i++)
        {
            m.row(i) = base_vectors.back();
            base_vectors.pop_back();
        }
    }

    Eigen::VectorXd orthogonal(n_dimensions);

    for(uint i = 0; i < n_dimensions; i++)
    {
        m(0, i) = 1;
        orthogonal(i) = m.determinant();
        m(0, i) = 0;
    }

    return orthogonal.normalized();
}

/*
 * create rotation matrix to map cohyperplanar vectors to a single dimension lower
 * apply rotation matrix to all vectors
 * return rotated vectors with the now orthogonal axis dropped along with the inverse rotation matrix
 */
pair<vector<Eigen::VectorXd>, Eigen::MatrixXd> reduce_dimension(const vector<Eigen::VectorXd>& vectors)
{
    uint n_dimensions = vectors[0].rows();

    Eigen::VectorXd orthogonal = orthogonal_vector(vectors);
    Eigen::VectorXd base = base_vector(n_dimensions);
    Eigen::MatrixXd rotation = rotation_matrix(orthogonal, base);

    vector<Eigen::VectorXd> reduced_vectors;
    reduced_vectors.reserve(vectors.size());

    for(const Eigen::VectorXd& v : vectors)
    {
        reduced_vectors.push_back((rotation * v).tail(n_dimensions - 1));
    }

    return std::make_pair(reduced_vectors, rotation.transpose());
}

/*
 * inverse of reduce_dimension
 */
vector<Eigen::VectorXd> increase_dimension(const vector<Eigen::VectorXd>& vectors, const Eigen::MatrixXd& rotation)
{
    if(vectors.size() == 0)
    {
        return vector<Eigen::VectorXd>();
    }

    uint n_dimensions = vectors[0].size();

    vector<Eigen::VectorXd> increased_vectors;
    increased_vectors.reserve(vectors.size());

    for(const Eigen::VectorXd& vector : vectors)
    {
        Eigen::VectorXd increased_vector(n_dimensions + 1);
        increased_vector.setZero();
        increased_vector.tail(n_dimensions) = vector;
        increased_vector = rotation * increased_vector;

        increased_vectors.push_back(increased_vector);
    }

    return increased_vectors;
}

/*
 * calculate the vector to the circumcenter of this simplex
 * all vertices are given as vectors from a single vertex
 * therefore, one vertex of the simplex already lies at the origin
 */
Eigen::VectorXd circumcenter_simplex(const vector<Eigen::VectorXd>& simplex)
{
    uint n_dimensions = simplex[0].rows();

    if(simplex.size() < n_dimensions)
    {
        vector<Eigen::VectorXd> reduced_simplex;
        Eigen::MatrixXd rotation;

        std::tie(reduced_simplex, rotation) = reduce_dimension(simplex);

        Eigen::VectorXd reduced_circumcenter = circumcenter_simplex(reduced_simplex);

        return increase_dimension({ reduced_circumcenter }, rotation)[0];
    }

    Eigen::VectorXd circumcenter(n_dimensions);
    Eigen::MatrixXd m(n_dimensions + 1, n_dimensions + 1);
    m.setZero();

    for(uint i = 0; i < n_dimensions; i++)
    {
        m(i + 1, 0) = simplex[i].squaredNorm();
        m.row(i + 1).tail(n_dimensions) = simplex[i];
    }

    for(uint i = 0; i < n_dimensions; i++)
    {
        m(0, i + 1) = 1;
        circumcenter(i) = -m.determinant();
        m(0, i + 1) = 0;
    }

    m.resize(n_dimensions, n_dimensions);
    m.setZero();

    for(uint i = 0; i < n_dimensions; i++)
    {
        m.row(i) = simplex[i];
    }

    return circumcenter / (2 * m.determinant());
}

bool nsphere_valid(Eigen::VectorXd center, double radius)
{
    for(uint i = 0; i < center.size(); i++)
    {
        if(center(i) - radius < 0 || center(i) + radius > 1)
        {
            return false;
        }
    }

    return true;
}

int main(int argc, char** argv)
{
    uint n_dimensions, n_points;

    cin >> n_dimensions >> n_points;
    
    vector<Eigen::VectorXd> points;
    points.reserve(n_points);

    for(uint i = 0; i < n_points; i++)
    {
        Eigen::VectorXd point(n_dimensions);

        for(uint j = 0; j < n_dimensions; j++)
        {
            cin >> point[j];
        }

        points.push_back(point);
    }

    vector<bool> choosing_vector(n_points, 0);

    double min_radius = 1.0 / 0.0;
    Eigen::VectorXd min_center(n_dimensions);
    min_center.fill(0);

    vector<uint> chosen_indices;
    chosen_indices.reserve(n_chosen);
    vector<Eigen::VectorXd> chosen_points;
    chosen_points.reserve(n_chosen);

    for(uint n_chosen = 1; n_chosen <= n_dimensions + 1 && n_chosen <= n_points / 2; n_chosen++)
    {
        choosing_vector.rbegin()[n_chosen - 1] = 1;

        do
        {
            chosen_indices.clear();
            chosen_points.clear();

            for(uint i = 0; i < n_points; i++)
            {
                if(choosing_vector[i])
                {
                    chosen_indices.push_back(i);
                    chosen_points.push_back(points[i]);
                }
            }

            Eigen::VectorXd& origin = chosen_points[0];
            vector<Eigen::VectorXd> facet;
            facet.reserve(n_chosen - 1);

            for(uint j = 1; j < chosen_points.size(); j++)
            {
                facet.push_back(chosen_points[j] - origin);
            }

            Eigen::VectorXd relative_circumcenter = circumcenter_simplex(facet);

            Eigen::VectorXd nsphere_center = origin + relative_circumcenter;
            double nsphere_radius = relative_circumcenter.norm();

            if(nsphere_radius < min_radius && nsphere_valid(nsphere_center, nsphere_radius))
            {
                uint points_contained = 0;

                for(uint i = 0; i < n_points && points_contained <= n_points / 2; i++)
                {
                    if((nsphere_center - points[i]).norm() <= nsphere_radius + epsilon)
                    {
                        points_contained++;
                    }
                }

                if(points_contained == n_points / 2)
                {
                    min_radius = nsphere_radius;
                    min_center = nsphere_center;
                }
            }
        }
        while(std::next_permutation(choosing_vector.begin(), choosing_vector.end()));
    }

    if(min_radius != 1.0 / 0.0)
    {
        cout.precision(8);
        cout << fixed;
        
        print_vector(min_center);
        cout << min_radius << endl;
    }
    else
    {
        cout << "No solution" << endl;
    }
}