Given n points within the unit square, calculate the smallest circle that lies completely within the unit square and contains exactly half of the given points.

main.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define PRECISION "10"

typedef struct __point point;
typedef struct __circle circle;

struct __point
{
    double x, y;
};

struct __circle
{
    double x, y, r;
};

unsigned next_permutation(unsigned* array, unsigned length)
{
    unsigned i = length - 1;
    while(i > 0 && array[i - 1] >= array[i])
    {
        i--;
    }

    if(i == 0)
    {
        return 0;
    }

    unsigned j = length - 1;
    while(array[j] <= array[i - 1])
    {
        j--;
    }

    unsigned temp = array[i - 1];
    array[i - 1] = array[j];
    array[j] = temp;

    j = length - 1;
    while(i < j)
    {
        temp = array[i];
        array[i] = array[j];
        array[j] = temp;
        i++;
        j--;
    }

    return 1;
}

void fprintarray(FILE* file, unsigned* array, unsigned length)
{
    fprintf(file, "[ ");

    for(unsigned i = 0; i < length; i++)
    {
        fprintf(file, "%u ", array[i]);
    }

    fprintf(file, "]\n");
}

inline double dist(point a, point b)
{
    return hypot(a.x - b.x, a.y - b.y);
}

inline circle smallest_circle(point a, point b)
{
    return (circle) { (a.x + b.x) / 2, (a.y + b.y) / 2, dist(a, b) / 2 };
}

inline circle circumcircle(point a, point b, point c)
{
    // Thanks Wikipedia
    double d = 2 * (a.x * (b.y - c.y) + b.x * (c.y - a.y) + c.x * (a.y - b.y));

    double x = ((a.x * a.x + a.y * a.y) * (b.y - c.y) + (b.x * b.x + b.y * b.y) * (c.y - a.y) + (c.x * c.x + c.y * c.y) * (a.y - b.y)) / d;
    double y = ((a.x * a.x + a.y * a.y) * (c.x - b.x) + (b.x * b.x + b.y * b.y) * (a.x - c.x) + (c.x * c.x + c.y * c.y) * (b.x - a.x)) / d;
    double r = dist((point) { x, y }, a);

    return (circle) { x, y, r };
}

inline unsigned circle_contains(circle c, point p)
{
    return hypot(c.x - p.x, c.y - p.y) <= c.r;
}

inline unsigned valid_circle(circle c)
{
    return c.x - c.r >= 0.0 && c.x + c.r <= 1.0 && c.y - c.r >= 0.0 && c.y + c.r <= 1.0;
}

/* 
 * point_count < 6
 */
circle small_input(point* points, unsigned point_count)
{
    if(point_count == 2)
    {
        return (circle) { points[0].x, points[0].y, 0.0 };
    }
    else if(point_count == 4)
    {
        circle* circles = malloc(6 * sizeof(circle));

        circles[0] = smallest_circle(points[0], points[1]);
        circles[1] = smallest_circle(points[0], points[2]);
        circles[2] = smallest_circle(points[0], points[3]);
        circles[3] = smallest_circle(points[1], points[2]);
        circles[4] = smallest_circle(points[1], points[3]);
        circles[5] = smallest_circle(points[2], points[3]);

        unsigned* valid = malloc(6 * sizeof(unsigned));

        valid[0] = !circle_contains(circles[0], points[2]) && !circle_contains(circles[0], points[3]) && valid_circle(circles[0]);
        valid[1] = !circle_contains(circles[1], points[1]) && !circle_contains(circles[1], points[3]) && valid_circle(circles[1]);
        valid[2] = !circle_contains(circles[2], points[1]) && !circle_contains(circles[2], points[2]) && valid_circle(circles[2]);
        valid[3] = !circle_contains(circles[3], points[0]) && !circle_contains(circles[3], points[3]) && valid_circle(circles[3]);
        valid[4] = !circle_contains(circles[4], points[0]) && !circle_contains(circles[4], points[2]) && valid_circle(circles[4]);
        valid[5] = !circle_contains(circles[5], points[0]) && !circle_contains(circles[5], points[1]) && valid_circle(circles[5]);

        unsigned smallest_i = -1;
        double smallest_r = 1.0 / 0.0;

        for(unsigned i = 0; i < 6; i++)
        {
            fprintf(stdout, "%.5f %.5f %.5f\n", circles[i].x, circles[i].y, circles[i].r);

            if(valid[i] && circles[i].r < smallest_r)
            {
                smallest_i = i;
                smallest_r = circles[i].r;
            }
        }

        circle answer;

        if(smallest_i == -1)
        {
            answer = (circle) { 0.0, 0.0, -1.0 };
        }
        else
        {
            answer = circles[smallest_i];
        }

        free(circles);
        free(valid);

        return answer;
    }

    return (circle) { 0, 0, -1.0 / 0.0 };
}

/* 
 * point_count >= 6
 */
circle large_input(point* points, unsigned point_count)
{
    unsigned* triangle_combination = malloc(point_count * sizeof(unsigned));
    circle global_smallest = (circle) { 0.0, 0.0, 1.0 / 0.0 };
    unsigned smallest_found = 0;

    for(unsigned i = 0; i < point_count; i++)
    {
        triangle_combination[i] = point_count - i <= 3;
    }

    point* triangle_selection = malloc(3 * sizeof(point));

    do
    {
        for(unsigned i = 0, j = 0; i < point_count && j < 3; i++)
        {
            if(triangle_combination[i])
            {
                triangle_selection[j++] = points[i];
            }
        }

        double side_min_1, side_min_2, side_max;
        unsigned side_max_label;

        {
            double side_a = dist(triangle_selection[0], triangle_selection[1]);
            double side_b = dist(triangle_selection[0], triangle_selection[2]);
            double side_c = dist(triangle_selection[1], triangle_selection[2]);

            if(side_a > side_b)
            {
                if(side_a > side_c)
                {
                    side_min_1 = side_b;
                    side_min_2 = side_c;
                    side_max = side_a;
                    side_max_label = 0;
                }
                else
                {
                    side_min_1 = side_a;
                    side_min_2 = side_b;
                    side_max = side_c;
                    side_max_label = 2;
                }
            }
            else
            {
                if(side_b > side_c)
                {
                    side_min_1 = side_a;
                    side_min_2 = side_c;
                    side_max = side_b;
                    side_max_label = 1;
                }
                else
                {
                    side_min_1 = side_a;
                    side_min_2 = side_b;
                    side_max = side_c;
                    side_max_label = 2;
                }
            }
        }

        circle local_smallest;

        if(hypot(side_min_1, side_min_2) <= side_max)
        {
            switch(side_max_label)
            {
                case 0:
                {
                    local_smallest = smallest_circle(triangle_selection[0], triangle_selection[1]);
                    break;
                }
                case 1:
                {
                    local_smallest = smallest_circle(triangle_selection[0], triangle_selection[2]);
                    break;
                }
                case 2:
                {
                    local_smallest = smallest_circle(triangle_selection[1], triangle_selection[2]);
                    break;
                }
            }
        }
        else
        {
            local_smallest = circumcircle(triangle_selection[0], triangle_selection[1], triangle_selection[2]);
        }

        if(!valid_circle(local_smallest) || local_smallest.r >= global_smallest.r)
        {
            continue;
        }

        unsigned in_circle = 3;

        for(unsigned i = 0; i < point_count && in_circle <= point_count / 2; i++)
        {
            if(!triangle_combination[i])
            {
                in_circle += circle_contains(local_smallest, points[i]);
            }
        }

        if(in_circle != point_count / 2)
        {
            continue;
        }

        global_smallest = local_smallest;
        smallest_found = 1;
    }
    while(next_permutation(triangle_combination, point_count));

    free(triangle_combination);
    free(triangle_selection);

    if(smallest_found)
    {
        return global_smallest;
    }
    else
    {
        return (circle) { 0.0, 0.0, -1.0 };
    }
}

int main(int argc, char** argv)
{
    point* points;
    unsigned point_count;

    fscanf(stdin, "%u", &point_count);
    points = malloc(point_count * sizeof(point));

    for(unsigned i = 0; i < point_count; i++)
    {
        fscanf(stdin, "%lf %lf", &points[i].x, &points[i].y);
    }

    circle answer;

    /* 
     * large_input requires at least 6 points, as the algorithm only checks circles with 3 or more points within them
     */
    if(point_count < 6)
    {
        answer = small_input(points, point_count);
    }
    else
    {
        answer = large_input(points, point_count);
    }

    free(points);

    if(answer.r >= 0)
    {
        fprintf(stdout, "%." PRECISION "lf %." PRECISION "lf\n%." PRECISION "lf\n", answer.x, answer.y, answer.r);
    }
    else
    {
        fprintf(stdout, "No solution found\n");
    }
}