#pragma GCC optimize ("fast-math")
#include<bits/stdc++.h>
using namespace std;
using LL=long long;
#define FOR(i,l,r)for(int i=(l);i<=(r);++i)
#define REP(i,n)FOR(i,0,(n)-1)
#define ssize(x)int(x.size())
#ifdef DEBUG
auto operator<<(auto&o,auto x)->decltype(x.end(),o);
auto&operator<<(auto&o,pair<auto,auto>p){return o<<"("<<p.first<<", "<<p.second<<")";}
auto&operator<<(auto&o,tuple<auto,auto,auto>t){return o<<"("<<get<0>(t)<<", "<<get<1>(t)<<", "<<get<2>(t)<<")";}
auto&operator<<(auto&o,tuple<auto,auto,auto,auto>t){return o<<"("<<get<0>(t)<<", "<<get<1>(t)<<", "<<get<2>(t)<<", "<<get<3>(t)<<")";}
auto operator<<(auto&o,auto x)->decltype(x.end(),o){o<<"{";int i=0;for(auto e:x)o<<","+!i++<<e;return o<<"}";}
#define debug(X...)cerr<<"["#X"]: ",[](auto...$){((cerr<<$<<"; "),...)<<endl;}(X)
#else
#define debug(...){}
#endif

using T = std::chrono::steady_clock::time_point;

const T start = std::chrono::steady_clock::now();

bool finish() {
	T end = std::chrono::steady_clock::now();
	auto value = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
	const int hard_limit = 2800;
	return value > hard_limit;
}

using D = double;
using V = vector<D>;
const D INF = 1e18;

using Complex = complex<D>;
void fft(vector<Complex> &a) {
	int n = ssize(a), L = 31 - __builtin_clz(n);
	static vector<complex<long double>> R(2, 1);
	static vector<Complex> rt(2, 1);
	for(static int k = 2; k < n; k *= 2) {
		R.resize(n), rt.resize(n);
		auto x = polar(1.0L, acosl(-1) / k);
		FOR(i, k, 2 * k - 1)
			rt[i] = R[i] = i & 1 ? R[i / 2] * x : R[i / 2];
	}
	vector<int> rev(n);
	REP(i, n) rev[i] = (rev[i / 2] | (i & 1) << L) / 2;
	REP(i, n) if(i < rev[i]) swap(a[i], a[rev[i]]);
	for(int k = 1; k < n; k *= 2) {
		for(int i = 0; i < n; i += 2 * k) REP(j, k) {
			Complex z = rt[j + k] * a[i + j + k]; // mozna zoptowac rozpisujac
			a[i + j + k] = a[i + j] - z;
			a[i + j] += z;
		}
	}
}
V conv(V a, V b) {
	if(a.empty() || b.empty()) return {};
	V res(ssize(a) + ssize(b) - 1);
	if (min(ssize(a), ssize(b)) <= 300) {
		REP(i, ssize(a))
			REP(j, ssize(b))
				res[i + j] += a[i] * b[j];
		return res;
	}
	int L = 32 - __builtin_clz(ssize(res)), n = (1 << L);
	vector<Complex> in(n), out(n);
	copy(a.begin(), a.end(), in.begin());
	REP(i, ssize(b)) in[i].imag(b[i]);
	fft(in);
	for(auto &x : in) x *= x;
	REP(i, n) out[i] = in[-i & (n - 1)] - conj(in[i]);
	fft(out);
	REP(i, ssize(res)) res[i] = imag(out[i]) / (4 * n);
	return res;
}
V conv_many(vector<V> v) {
	assert(ssize(v));
	v.reserve(2 * ssize(v));
	multiset<pair<int, int>> s;
	REP(i, ssize(v)) {
		s.emplace(ssize(v[i]), i);
	}
	while (ssize(s) > 1) {
		auto [_, a] = *s.begin();
		s.erase(s.begin());
		auto [__, b] = *s.begin();
		s.erase(s.begin());
		v.emplace_back(conv(v[a], v[b]));
		s.emplace(ssize(v.back()), ssize(v) - 1);
	}
	auto ret = v[s.begin()->second];
	return ret;
}

D solve(V v, int t) {
	sort(v.rbegin(), v.rend());
	const int n = ssize(v);
	debug(n, t, v);

	auto pref = v;
	REP(i, n - 1)
		pref[i + 1] += pref[i];

	const auto estimate = [&](int k) -> D {
		D expected = pref[k - 1];
		const int minimum = (t + k + 1) / 2;
		return (expected - minimum) / sqrt(k);
	};
	const auto find_peak = [&]() -> int {
		pair<D, int> best = {-INF, INF};
		FOR(k, 1, n) {
			best = max(best, pair(estimate(k), -k));
		}
		return -best.second;
	};

	const auto adjust_range = [&](int l, int r) -> pair<int, int> {
		if (l < t) {
			const int diff = t - l;
			l += diff;
			r += diff;
		}
		if (r > n) {
			const int diff = r - n;
			r -= diff;
			l -= diff;
		}
		l = max(l, t);
		r = min(r, n);
		return {l, r};
	};

	D ret = 0;

	const auto eval = [&](const V& poly, int k) {
		if (k % 2 != t % 2)
			return;
		const int minimum = (t + k + 1) / 2;
		D value = 0;
		FOR(d, minimum, k) {
			value += poly[d];
		}
		ret = max(ret, value);
	};

	const auto calc = [&](int k) -> V {
		vector<V> h;
		REP(i, k) {
			h.emplace_back(V{1 - v[i], v[i]});
		}
		auto poly = conv_many(h);
		eval(poly, k);
		return poly;
	};

	const auto consider_interval = [&](int l, int r) {
		const auto [left, right] = adjust_range(l, r);

		auto poly = calc(left);
		poly.resize(n + 1);

		FOR(k, left + 1, right) {
			const D p = v[k - 1];
			const D q = 1 - p;
			for (int i = k; i >= 1; --i) {
				poly[i] = poly[i] * q + poly[i - 1] * p;
			}
			poly[0] *= q;
			eval(poly, k);
			if (k % 10 == 0) {
				if (finish())
					return;
			}
		}
	};

	const int pref_length = 1000;
	consider_interval(t, t + pref_length);

	const int length = 10000;
	const int peak = find_peak();
	consider_interval(peak - length / 2, peak + length / 2);

	return ret;
}

D solve_quadratic(vector<D> v, int t) {
	sort(v.rbegin(), v.rend());
	const int n = ssize(v);
	debug(n, t, v);
	D ret = 0;
	V poly(n + 1);
	poly[0] = 1;
	FOR(k, 1, n) {
		const D p = v[k - 1];
		const D q = 1 - p;
		for (int i = k; i >= 1; --i) {
			poly[i] = poly[i] * q + poly[i - 1] * p;
		}
		poly[0] *= q;
		D current = 0;
		const int minimum = (t + k + 1) / 2;
		FOR(d, minimum, k) {
			current += poly[d];
		}
		debug(k, current);
		ret = max(ret, current);
		if (k % 10 == 0) {
			if (finish())
				return ret;
		}
	}
	return ret;
}

int main() {
	cin.tie(0)->sync_with_stdio(0);

	int n, t;
	cin >> n >> t;
	vector<D> v(n);
	for (auto& x : v)
		cin >> x;
	const int quadratic_limit = 17'000;
	D answer;
	if (n <= quadratic_limit) {
		answer = solve_quadratic(v, t);
	}
	else {
		answer = solve(v, t);
	}
	answer = max(answer, D(0));
	answer = min(answer, D(1));
#ifdef TESTS
	const int precision = 7;
#else
	const int precision = 9;
#endif
	cout << setprecision(precision) << fixed;
	cout << answer << '\n';
}