Code

#include <iostream>
#include <iomanip>
#include <sstream>
#include <fstream>
#include <string>
#include <cstdlib>
#include <cstdio>
#include <climits>
#include <cstring>
#include <cassert>
#include <ctime>
#include <cmath>
#include <vector>
#include <queue>
#include <stack>
#include <list>
#include <set>
#include <map>
#include <bitset>
#include <algorithm>
#include <utility>
#include <numeric>
#include <functional>

#define forn(i, n) for (int i = 0; i < int(n); i++)
#define forl(i, n) for (int i = 1; i <= int(n); i++)
#define ford(i, n) for (int i = int(n) - 1; i >= 0; i--)
#define fore(i, l, r) for (int i = int(l); i <= int(r); i++)
#define correct(x, y, n, m) (0 <= (x) && (x) < (n) && 0 <= (y) && (y) < (m))
#define all(a) (a).begin(), (a).end()
#define sz(a) int((a).size())
#define pb(a) push_back(a)
#define mp(x, y) make_pair((x), (y))
#define ft first
#define sc second
#define x first
#define y second
#define eprintf(...) fprintf(stderr, __VA_ARGS__)

using namespace std;

typedef long long li;
typedef long double ld;
typedef pair<int, int> pti;

template<typename X> inline X abs(const X& a) { return a < 0? -a: a; }
template<typename X> inline X sqr(const X& a) { return a * a; }

const int INF = int(1e9);
const li INF64 = li(1e18);
const ld EPS = 1e-9, PI = 3.1415926535897932384626433832795;
const ld inf = 1e100;

struct pt
{
    ld x, y;
    pt() {}
    pt(ld x, ld y): x(x), y(y) {}
};

inline bool operator< (const pt& a, const pt& b)
{
    if (abs(a.x - b.x) > EPS) return a.x < b.x;
    if (abs(a.y - b.y) > EPS) return a.y < b.y;
    return false;
}

inline bool operator== (const pt& a, const pt& b)
{
    if (abs(a.x - b.x) > EPS) return false;
    if (abs(a.y - b.y) > EPS) return false;
    return true;
}

inline pt operator- (const pt& a, const pt& b) { return pt(a.x - b.x, a.y - b.y); }
inline pt operator+ (const pt& a, const pt& b) { return pt(a.x + b.x, a.y + b.y); }
inline pt operator* (const pt& a, const ld& b) { return pt(a.x * b, a.y * b); }
inline pt operator/ (const pt& a, const ld& b) { return pt(a.x / b, a.y / b); }

inline ld len(const pt& v) { return sqrt(sqr(v.x) + sqr(v.y)); }
inline ld dist(const pt& a, const pt& b) { return len(a - b); }

inline pt normal(const pt& v) { return len(v) < EPS ? v : v / len(v); }
inline pt normal(const pt& v, const ld& l) { return normal(v) * l; }

inline ld dot(const pt& a, const pt&b) { return a.x * b.x + a.y * b.y; }
inline ld cross(const pt& a, const pt&b) { return a.x * b.y - a.y * b.x; }

const int N = 100 + 3;

int n;
pt a[N];

inline bool read()
{
    if (scanf("%d", &n) != 1)
        return false;

    forn(i, n)
    {
        double x, y;
        assert(scanf("%lf%lf", &x, &y) == 2);

        a[i] = pt(x, y);
    }

    return true;
}

int szr;
ld r[2];

inline void solve(const ld& a, const ld& b, const ld& c)
{
    szr = 0;

    if (abs(a) > EPS)
    {
        ld D = sqr(b) - 4 * a * c;
        if (D + EPS < 0) return;

        D = sqrt(max(ld(0), D));
        r[szr++] = (-b + D) / (2 * a);
        r[szr++] = (-b - D) / (2 * a);
        return;
    }

    if (abs(b) > EPS)
    {
        r[szr++] = -c / b;
        return;
    }

    if (abs(c) < EPS)
    {
        r[szr++] = 0;
        return;
    }
}

inline ld solve(ld sx, ld sy, ld sxy, ld sx2, ld sy2, int n)
{
    if (n == 0) return 0;

    pt c(sx / n, sy / n);

    sxy = sxy + c.x * c.y * n - c.x * sy - c.y * sx;
    sx2 = sx2 + sqr(c.x) * n - 2 * c.x * sx;
    sy2 = sy2 + sqr(c.y) * n - 2 * c.y * sy;
    sx = sy = 0;

    ld res = inf;

    forn(iter, 2)
    {
        solve(sxy, (sy2 - sx2) * (iter == 0 ? 1 : -1), -sxy);

        forn(i, szr)
        {
            const ld A = (iter == 0 ? r[i] : 1), B = (iter == 0 ? 1 : r[i]);
            res = min(res, (sqr(A) * sx2 + 2 * A * B * sxy + sqr(B) * sy2) / (sqr(A) + sqr(B)));
        }
    }

    assert(res < inf / 2);
    return res;
}

pt pole, dir;
ld A, B, C;

inline ld proj(const pt& a) { return dot(a - pole, dir); }

inline int sign(const pt& p)
{
    ld val = A * p.x + B * p.y + C;
    return val + EPS < 0 ? -1 : val > EPS ? 1 : 0;
}

inline bool cmp(const pt& a, const pt& b)
{
    ld v1 = proj(a), v2 = proj(b);
    if (abs(v1 - v2) > EPS) return v1 < v2;
    return false;
}

pt p[N];

ld x[2][N];
ld y[2][N];
ld xy[2][N];
ld x2[2][N];
ld y2[2][N];
int m[2][N];

inline ld sumX(int up, int l, int r)
{
    if (l > r) return 0;
    return x[up][r] - (l == 0 ? 0 : x[up][l - 1]);
}

inline ld sumY(int up, int l, int r)
{
    if (l > r) return 0;
    return y[up][r] - (l == 0 ? 0 : y[up][l - 1]);
}

inline ld sumXY(int up, int l, int r)
{
    if (l > r) return 0;
    return xy[up][r] - (l == 0 ? 0 : xy[up][l - 1]);
}

inline ld sumX2(int up, int l, int r)
{
    if (l > r) return 0;
    return x2[up][r] - (l == 0 ? 0 : x2[up][l - 1]);
}

inline ld sumY2(int up, int l, int r)
{
    if (l > r) return 0;
    return y2[up][r] - (l == 0 ? 0 : y2[up][l - 1]);
}

inline ld sumM(int up, int l, int r)
{
    if (l > r) return 0;
    return m[up][r] - (l == 0 ? 0 : m[up][l - 1]);
}

inline void clearSum()
{
    forn(i, 2)
        forn(j, n)
        {
            x[i][j] = 0;
            y[i][j] = 0;
            xy[i][j] = 0;
            x2[i][j] = 0;
            y2[i][j] = 0;
            m[i][j] = 0;
        }
}

inline void makeSum()
{
    forn(i, 2)
        forl(j, n - 1)
        {
            x[i][j] += x[i][j - 1];
            y[i][j] += y[i][j - 1];
            xy[i][j] += xy[i][j - 1];
            x2[i][j] += x2[i][j - 1];
            y2[i][j] += y2[i][j - 1];
            m[i][j] += m[i][j - 1];
        }
}

int szv;
pt v[N];

const int MASK = (1 << 4) + 3;

ld mx[MASK];
ld my[MASK];
ld mxy[MASK];
ld mx2[MASK];
ld my2[MASK];
int mm[MASK];

int least[MASK];

inline void solve()
{
    ld SX = 0, SY = 0, SXY = 0, SX2 = 0, SY2 = 0;
    int N = n;

    forn(i, n)
    {
        SX += a[i].x;
        SY += a[i].y;
        SXY += a[i].x * a[i].y;
        SX2 += sqr(a[i].x);
        SY2 += sqr(a[i].y);
    }

    ld ans = inf;

    forn(ii, n)
        forn(jj, ii)
        {
            pole = a[ii], dir = normal(a[jj] - a[ii]);
            A = -dir.y, B = dir.x, C = -(A * pole.x + B * pole.y);

            {
                ld v = sqr(A) + sqr(B);
                A /= v, B /= v, C /= v;
            }

            forn(i, n) p[i] = a[i];
            sort(p, p + n, cmp);

            clearSum();

            forn(i, n)
            {
                int s = sign(p[i]);
                if (s == 0) continue;

                int up = s > 0;

                x[up][i] = p[i].x;
                y[up][i] = p[i].y;
                xy[up][i] = p[i].x * p[i].y;
                x2[up][i] = sqr(p[i].x);
                y2[up][i] = sqr(p[i].y);
                m[up][i] = 1;
            }

            makeSum();

            for (int i = 0, j = 0; i < n; i = j)
            {
                szv = 0;
                v[szv++] = a[ii], v[szv++] = a[jj];

                while (abs(proj(p[i]) - proj(p[j])) < EPS)
                {
                    v[szv++] = p[j];
                    j++;
                }

                sort(v, v + szv);
                szv = int(unique(v, v + szv) - v);

                assert(szv >= 2 && szv <= 4);

                mx[0] = my[0] = mxy[0] = mx2[0] = my2[0] = mm[0] = 0;

                forl(mask, (1 << szv) - 1)
                {
                    int idx = least[mask];
                    int pmask = mask ^ (1 << idx);

                    mx[mask] = mx[pmask] + v[idx].x;
                    my[mask] = my[pmask] + v[idx].y;
                    mxy[mask] = mxy[pmask] + v[idx].x * v[idx].y;
                    mx2[mask] = mx2[pmask] + sqr(v[idx].x);
                    my2[mask] = my2[pmask] + sqr(v[idx].y);
                    mm[mask] = mm[pmask] + 1;
                }

                ld curx = sumX(0, 0, i - 1) + sumX(1, j, n - 1);
                ld cury = sumY(0, 0, i - 1) + sumY(1, j, n - 1);
                ld curxy = sumXY(0, 0, i - 1) + sumXY(1, j, n - 1);
                ld curx2 = sumX2(0, 0, i - 1) + sumX2(1, j, n - 1);
                ld cury2 = sumY2(0, 0, i - 1) + sumY2(1, j, n - 1);
                int curn = sumM(0, 0, i - 1) + sumM(1, j, n - 1);

                forn(mask, 1 << szv)
                {
                    curx += mx[mask];
                    cury += my[mask];
                    curxy += mxy[mask];
                    curx2 += mx2[mask];
                    cury2 += my2[mask];
                    curn += mm[mask];

                    ans = min(ans, solve(curx, cury, curxy, curx2, cury2, curn) + solve(SX - curx, SY - cury, SXY - curxy, SX2 - curx2, SY2 - cury2, N - curn));

                    curx -= mx[mask];
                    cury -= my[mask];
                    curxy -= mxy[mask];
                    curx2 -= mx2[mask];
                    cury2 -= my2[mask];
                    curn -= mm[mask];
                }
            }
        }

    assert(ans < inf / 2);
    cout << ans / n << endl;
}

int main()
{
#ifdef SU2_PROJ
    assert(freopen("input.txt", "rt", stdin));
    //assert(freopen("output.txt", "wt", stdout));
#endif

    forl(i, MASK - 1)
    {
        forn(j, 4)
        {
            if (i & (1 << j))
            {
                least[i] = j;
                break;
            }
        }
    }

    cout << setprecision(10) << fixed;
    cerr << setprecision(10) << fixed;

    assert(read());
    solve();

    return 0;
}