#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=a;i<b;i++)
#define rrep(i,a,b) for(int i=a;i>=b;i--)
#define fore(i,a) for(auto &i:a)
#define all(x) (x).begin(),(x).end()
#pragma GCC optimize ("-O3")
using namespace std; void _main(); int main() { cin.tie(0); ios::sync_with_stdio(false); _main(); }
typedef long long ll; const int inf = INT_MAX / 2; const ll infl = 1LL << 60;
template<class T>bool chmax(T &a, const T &b) { if (a<b) { a = b; return 1; } return 0; }
template<class T>bool chmin(T &a, const T &b) { if (b<a) { a = b; return 1; } return 0; }
//---------------------------------------------------------------------------------------------------
const double EPS = 1e-8, INF = 1e12, PI = 2 * acos(0.0);
typedef complex<double> P;
namespace std {
    bool operator < (const P& a, const P& b) { return real(a) != real(b) ? real(a) < real(b) : imag(a) < imag(b); } 
    bool operator == (const P& a, const P& b) { return abs(real(a) - real(b)) < EPS && abs(imag(a) - imag(b)) < EPS; }
    P operator / (const P& a, const double& b) { return P(real(a) / b, imag(a) / b); }
    P operator * (const P& a, const double& b) { return P(real(a) * b, imag(a) * b); }
}
double cross(const P& a, const P& b) { return imag(conj(a)*b); }
double dot(const P& a, const P& b) { return real(conj(a)*b); }
struct L : public vector<P> { L() {} L(const P &a, const P &b) { push_back(a); push_back(b); } };
typedef vector<P> G;
struct C { P p; double r; C() {} C(const P &p, double r) : p(p), r(r) {} };
int ccw(P a, P b, P c) {
	b -= a; c -= a;
	if (cross(b, c) > 0)   return 0;       // counter clockwise
	if (cross(b, c) < 0)   return 1;       // clockwise
	if (dot(b, c) < 0)     return +2;       // c--a--b on line
	if (norm(b) < norm(c)) return -2;       // a--b--c on line
	return 0;
}
P rotate(P vec, double ang) {
    double x = real(vec), y = imag(vec);
    return P(x * cos(ang) - y * sin(ang), x * sin(ang) + y * cos(ang));
}
bool intersectLL(const L &l, const L &m) {
    return abs(cross(l[1] - l[0], m[1] - m[0])) > EPS || // non-parallel
        abs(cross(l[1] - l[0], m[0] - l[0])) < EPS;   // same line
}
bool intersectLS(const L &l, const L &s) {
    return cross(l[1] - l[0], s[0] - l[0])*       // s[0] is left of l
        cross(l[1] - l[0], s[1] - l[0]) < EPS; // s[1] is right of l
}
bool intersectLP(const L &l, const P &p) {
    return abs(cross(l[1] - p, l[0] - p)) < EPS;
}
bool intersectSS(const L &s, const L &t) {
    return ccw(s[0], s[1], t[0])*ccw(s[0], s[1], t[1]) <= 0 &&
        ccw(t[0], t[1], s[0])*ccw(t[0], t[1], s[1]) <= 0;
}
bool intersectSP(const L &s, const P &p) {
    return abs(s[0] - p) + abs(s[1] - p) - abs(s[1] - s[0]) < EPS; // triangle inequality
}
bool intersectSSwithoutPoint(const L &s, const L &t) { // not verified
    return ccw(s[0], s[1], t[0])*ccw(s[0], s[1], t[1]) < 0 &&
        ccw(t[0], t[1], s[0])*ccw(t[0], t[1], s[1]) < 0;
}
P crosspoint(const L &l, const L &m) {
    double A = cross(l[1] - l[0], m[1] - m[0]);
    double B = cross(l[1] - l[0], l[1] - m[0]);
    if (abs(A) < EPS && abs(B) < EPS) return m[0]; // same line
    if (abs(A) < EPS) assert(false); // !!!PRECONDITION NOT SATISFIED!!!
    return m[0] + B / A * (m[1] - m[0]);
}
double argument(const P &a, const P &b) { // argument for A->B[-PI,PI]
    double ax = real(a), ay = imag(a), bx = real(b), by = imag(b);
    return atan2(by - ay, bx - ax);
}
// A->B->C時計回りの時の外角
// argument for A->B->C
double threePointArgumentTokei(const P &a, const P &b, const P &c) {
    P ba = b - a;
    P cb = c - b;

    double r1 = atan2(real(ba), imag(ba)) - PI / 2.0;
    while (r1<-PI) r1 += 2 * PI;
    while (PI<r1) r1 -= 2 * PI;
    double r2 = atan2(real(cb), imag(cb)) - PI / 2.0;
    while (r2<-PI) r2 += 2 * PI;
    while (PI<r2) r2 -= 2 * PI;

    //cout << id << "\t" << r1 << "\t" << r2 << endl;

    double range = 0;
    if (r1 >= 0 && r2>0) range = r2 - r1;
    else if (r1<0 && r2 >= 0) range = -r1 + r2;
    else if (r1 >= 0 && r2<0) range = (PI - r1) + (PI + r2);
    else range = r2 - r1;

    return range;
}
// A-B-Cの成す各のうち小さい方
double threePointArgument(const P &a, const P &b, const P &c) {
    double x = argument(b, a);
    double y = argument(b, c);

    double res = x - y;
    while (res < 0) res += PI * 2;
    while (PI * 2 <= res) res -= PI * 2;

    return min(res, 2 * PI - res);
}
// 三点より中心導出
P three_point_circle(const P& a, const P& b, const P& c) {
    P x = 1.0 / conj(b - a), y = 1.0 / conj(c - a);
    return (y - x) / (conj(x)*y - x * conj(y)) + a;
}
/*---------------------------------------------------------------------------------------------------
　　　　　　　　　　　 ∧＿∧  
　　　　　 ∧＿∧ 　（´<_｀ ）　 Welcome to My Coding Space!
　　　　 （ ´_ゝ`）　/　 ⌒i     
　　　　／　　　＼　 　  |　|     
　　　 /　　 /￣￣￣￣/　　|  
　 ＿_(__ﾆつ/　    ＿/ .| .|＿＿＿＿  
　 　　　＼/＿＿＿＿/　（u　⊃  
---------------------------------------------------------------------------------------------------*/






int N, M;
P xy[14];
double dp[1 << 14], dp2[14][1<<14];
//---------------------------------------------------------------------------------------------------
using ld = double;
double smallestEnclosingCircle(vector<P> &v) {
    if (v.size() <= 1) return 0;

    auto c3 = [](const P &a, const P &b, const P &c) {
        ld A = norm(b - c);
        ld B = norm(c - a);
        ld Co = norm(a - b);
        ld S = abs(cross(b - a, c - a));
        P p = (A*(B + Co - A)*a + B * (Co + A - B)*b + Co * (A + B - Co)*c) / (4 * S*S);
        ld r = abs(p - a);
        return C(p, r);
    };

    auto c2 = [](const P &a, const P &b) {
        P c = (a + b) / (ld)2;
        ld r = abs(a - c);
        return C(c, r);
    };

    auto in_circle = [](const P &a, const C &c) {
        return norm(a - c.p) <= c.r *c.r + EPS;
    };

    int n = v.size();
    random_shuffle(v.begin(), v.end());
    C c = c2(v[0], v[1]);
    for (int i = 2; i < n; ++i) {
        if (!in_circle(v[i], c)) {
            c = c2(v[0], v[i]);
            for (int j = 1; j < i; ++j) {
                if (!in_circle(v[j], c)) {
                    c = c2(v[j], v[i]);
                    for (int k = 0; k < j; ++k) {
                        if (!in_circle(v[k], c)) {
                            c = c3(v[i], v[j], v[k]);
                        }
                    }
                }
            }
        }
    }
    return c.r;
}
//---------------------------------------------------------------------------------------------------
void _main() {
    cin >> N >> M;
    rep(i, 0, N) {
        double x, y;
        cin >> x >> y;
        xy[i] = P(x, y);
    }

    rep(msk, 0, 1 << N) {
        vector<P> v;
        rep(i, 0, N) if (msk & (1 << i)) v.push_back(xy[i]);

        dp[msk] = smallestEnclosingCircle(v);
    }

    //rep(msk, 0, 1 << N) printf("dp[%d] = %.10f\n", msk, dp[msk]);

    rep(i, 0, M + 1) rep(msk, 0, 1 << N) dp2[i][msk] = INF;
    dp2[0][0] = 0;
    rep(i, 1, M + 1) rep(msk, 0, 1 << N) {
        for (int j = msk; j>0; j = (j - 1)&msk) {
            chmin(dp2[i][msk], max(dp2[i - 1][msk - j], dp[j]));
        }
    }

    double ans = INF;
    rep(i, 0, M + 1) chmin(ans, dp2[i][(1 << N) - 1]);
    printf("%.12f\n", ans);
}

### Prompt

Construct a CPP code solution to the problem outlined:

Problem

Great Demon King Megumi wants to defeat the $ N $ hero who lives on the ground.
Megumi can cast explosion magic up to $ M $ times.
Explosion magic is magic that extinguishes heroes who exist within a radius of $ r $ around arbitrary coordinates.
The hero is very thin, so you don't have to consider the size.
All $ M $ explosion spells shall be cast with the same radius.
I decided to use the magic of the minimum radius necessary to annihilate the hero.
Minimize the size of the explosive magic radius.

Constraints

The input satisfies the following conditions.

* $ 1 \ leq M \ leq N \ leq 14 $
* $ 0 \ leq x_i $, $ y_i \ leq 10 ^ 5 $

Input


$ N $ $ M $
$ x_1 $ $ y_1 $
...
$ x_N $ $ y_N $


The input is given in the following format.
The number of heroes $ N $ and the number of explosions to cast $ M $ are given as integers on the first line.
From the second line onward, the coordinates $ x_i $ $, $ $ y_i $ of each hero $ i $ are given as integers.

Output

Output the minimum value of the explosive magic radius as a real number.
Must not contain absolute error greater than $ 10 ^ {-3} $.

Examples

Input

5 2
0 0
5 5
10 10
100 100
200 200


Output

70.710678118755


Input

10 5
321 675
4312 6534
312 532
412 6543
21 43
654 321
543 0
32 5
41 76
5 1


Output

169.824909833728


Input

14 3
312 342
4893 432
321 4
389 4
23 543
0 0
1 1
2 2
432 12
435 32
1 5
2 10
100 100
20 50


Output

218.087711712613


Input

5 2
0 0
0 0
0 0
0 0
0 0


Output

0.000000000001

### Response

```cpp
#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=a;i<b;i++)
#define rrep(i,a,b) for(int i=a;i>=b;i--)
#define fore(i,a) for(auto &i:a)
#define all(x) (x).begin(),(x).end()
#pragma GCC optimize ("-O3")
using namespace std; void _main(); int main() { cin.tie(0); ios::sync_with_stdio(false); _main(); }
typedef long long ll; const int inf = INT_MAX / 2; const ll infl = 1LL << 60;
template<class T>bool chmax(T &a, const T &b) { if (a<b) { a = b; return 1; } return 0; }
template<class T>bool chmin(T &a, const T &b) { if (b<a) { a = b; return 1; } return 0; }
//---------------------------------------------------------------------------------------------------
const double EPS = 1e-8, INF = 1e12, PI = 2 * acos(0.0);
typedef complex<double> P;
namespace std {
    bool operator < (const P& a, const P& b) { return real(a) != real(b) ? real(a) < real(b) : imag(a) < imag(b); } 
    bool operator == (const P& a, const P& b) { return abs(real(a) - real(b)) < EPS && abs(imag(a) - imag(b)) < EPS; }
    P operator / (const P& a, const double& b) { return P(real(a) / b, imag(a) / b); }
    P operator * (const P& a, const double& b) { return P(real(a) * b, imag(a) * b); }
}
double cross(const P& a, const P& b) { return imag(conj(a)*b); }
double dot(const P& a, const P& b) { return real(conj(a)*b); }
struct L : public vector<P> { L() {} L(const P &a, const P &b) { push_back(a); push_back(b); } };
typedef vector<P> G;
struct C { P p; double r; C() {} C(const P &p, double r) : p(p), r(r) {} };
int ccw(P a, P b, P c) {
	b -= a; c -= a;
	if (cross(b, c) > 0)   return 0;       // counter clockwise
	if (cross(b, c) < 0)   return 1;       // clockwise
	if (dot(b, c) < 0)     return +2;       // c--a--b on line
	if (norm(b) < norm(c)) return -2;       // a--b--c on line
	return 0;
}
P rotate(P vec, double ang) {
    double x = real(vec), y = imag(vec);
    return P(x * cos(ang) - y * sin(ang), x * sin(ang) + y * cos(ang));
}
bool intersectLL(const L &l, const L &m) {
    return abs(cross(l[1] - l[0], m[1] - m[0])) > EPS || // non-parallel
        abs(cross(l[1] - l[0], m[0] - l[0])) < EPS;   // same line
}
bool intersectLS(const L &l, const L &s) {
    return cross(l[1] - l[0], s[0] - l[0])*       // s[0] is left of l
        cross(l[1] - l[0], s[1] - l[0]) < EPS; // s[1] is right of l
}
bool intersectLP(const L &l, const P &p) {
    return abs(cross(l[1] - p, l[0] - p)) < EPS;
}
bool intersectSS(const L &s, const L &t) {
    return ccw(s[0], s[1], t[0])*ccw(s[0], s[1], t[1]) <= 0 &&
        ccw(t[0], t[1], s[0])*ccw(t[0], t[1], s[1]) <= 0;
}
bool intersectSP(const L &s, const P &p) {
    return abs(s[0] - p) + abs(s[1] - p) - abs(s[1] - s[0]) < EPS; // triangle inequality
}
bool intersectSSwithoutPoint(const L &s, const L &t) { // not verified
    return ccw(s[0], s[1], t[0])*ccw(s[0], s[1], t[1]) < 0 &&
        ccw(t[0], t[1], s[0])*ccw(t[0], t[1], s[1]) < 0;
}
P crosspoint(const L &l, const L &m) {
    double A = cross(l[1] - l[0], m[1] - m[0]);
    double B = cross(l[1] - l[0], l[1] - m[0]);
    if (abs(A) < EPS && abs(B) < EPS) return m[0]; // same line
    if (abs(A) < EPS) assert(false); // !!!PRECONDITION NOT SATISFIED!!!
    return m[0] + B / A * (m[1] - m[0]);
}
double argument(const P &a, const P &b) { // argument for A->B[-PI,PI]
    double ax = real(a), ay = imag(a), bx = real(b), by = imag(b);
    return atan2(by - ay, bx - ax);
}
// A->B->C時計回りの時の外角
// argument for A->B->C
double threePointArgumentTokei(const P &a, const P &b, const P &c) {
    P ba = b - a;
    P cb = c - b;

    double r1 = atan2(real(ba), imag(ba)) - PI / 2.0;
    while (r1<-PI) r1 += 2 * PI;
    while (PI<r1) r1 -= 2 * PI;
    double r2 = atan2(real(cb), imag(cb)) - PI / 2.0;
    while (r2<-PI) r2 += 2 * PI;
    while (PI<r2) r2 -= 2 * PI;

    //cout << id << "\t" << r1 << "\t" << r2 << endl;

    double range = 0;
    if (r1 >= 0 && r2>0) range = r2 - r1;
    else if (r1<0 && r2 >= 0) range = -r1 + r2;
    else if (r1 >= 0 && r2<0) range = (PI - r1) + (PI + r2);
    else range = r2 - r1;

    return range;
}
// A-B-Cの成す各のうち小さい方
double threePointArgument(const P &a, const P &b, const P &c) {
    double x = argument(b, a);
    double y = argument(b, c);

    double res = x - y;
    while (res < 0) res += PI * 2;
    while (PI * 2 <= res) res -= PI * 2;

    return min(res, 2 * PI - res);
}
// 三点より中心導出
P three_point_circle(const P& a, const P& b, const P& c) {
    P x = 1.0 / conj(b - a), y = 1.0 / conj(c - a);
    return (y - x) / (conj(x)*y - x * conj(y)) + a;
}
/*---------------------------------------------------------------------------------------------------
　　　　　　　　　　　 ∧＿∧  
　　　　　 ∧＿∧ 　（´<_｀ ）　 Welcome to My Coding Space!
　　　　 （ ´_ゝ`）　/　 ⌒i     
　　　　／　　　＼　 　  |　|     
　　　 /　　 /￣￣￣￣/　　|  
　 ＿_(__ﾆつ/　    ＿/ .| .|＿＿＿＿  
　 　　　＼/＿＿＿＿/　（u　⊃  
---------------------------------------------------------------------------------------------------*/






int N, M;
P xy[14];
double dp[1 << 14], dp2[14][1<<14];
//---------------------------------------------------------------------------------------------------
using ld = double;
double smallestEnclosingCircle(vector<P> &v) {
    if (v.size() <= 1) return 0;

    auto c3 = [](const P &a, const P &b, const P &c) {
        ld A = norm(b - c);
        ld B = norm(c - a);
        ld Co = norm(a - b);
        ld S = abs(cross(b - a, c - a));
        P p = (A*(B + Co - A)*a + B * (Co + A - B)*b + Co * (A + B - Co)*c) / (4 * S*S);
        ld r = abs(p - a);
        return C(p, r);
    };

    auto c2 = [](const P &a, const P &b) {
        P c = (a + b) / (ld)2;
        ld r = abs(a - c);
        return C(c, r);
    };

    auto in_circle = [](const P &a, const C &c) {
        return norm(a - c.p) <= c.r *c.r + EPS;
    };

    int n = v.size();
    random_shuffle(v.begin(), v.end());
    C c = c2(v[0], v[1]);
    for (int i = 2; i < n; ++i) {
        if (!in_circle(v[i], c)) {
            c = c2(v[0], v[i]);
            for (int j = 1; j < i; ++j) {
                if (!in_circle(v[j], c)) {
                    c = c2(v[j], v[i]);
                    for (int k = 0; k < j; ++k) {
                        if (!in_circle(v[k], c)) {
                            c = c3(v[i], v[j], v[k]);
                        }
                    }
                }
            }
        }
    }
    return c.r;
}
//---------------------------------------------------------------------------------------------------
void _main() {
    cin >> N >> M;
    rep(i, 0, N) {
        double x, y;
        cin >> x >> y;
        xy[i] = P(x, y);
    }

    rep(msk, 0, 1 << N) {
        vector<P> v;
        rep(i, 0, N) if (msk & (1 << i)) v.push_back(xy[i]);

        dp[msk] = smallestEnclosingCircle(v);
    }

    //rep(msk, 0, 1 << N) printf("dp[%d] = %.10f\n", msk, dp[msk]);

    rep(i, 0, M + 1) rep(msk, 0, 1 << N) dp2[i][msk] = INF;
    dp2[0][0] = 0;
    rep(i, 1, M + 1) rep(msk, 0, 1 << N) {
        for (int j = msk; j>0; j = (j - 1)&msk) {
            chmin(dp2[i][msk], max(dp2[i - 1][msk - j], dp[j]));
        }
    }

    double ans = INF;
    rep(i, 0, M + 1) chmin(ans, dp2[i][(1 << N) - 1]);
    printf("%.12f\n", ans);
}

```