1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
#pragma GCC optimize("O3")
#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define ld double
#define pb push_back
#define ff first
#define ss second
#define MOD 1000000009
#define INF 1000000019
#define INFL 1000000000000000099LL
#define IP 26




#define rep(i, a, b) for (int i = (a); i < (b); ++i)
#define sz(x) ((int)(x).size())
#define all(x) (x).begin(), (x).end()
typedef vector<int> vi;


/**
 * Author: Ludo Pulles, chilli, Simon Lindholm
 * Date: 2019-01-09
 * License: CC0
 * Source: http://neerc.ifmo.ru/trains/toulouse/2017/fft2.pdf (do read, it's excellent)
   Accuracy bound from http://www.daemonology.net/papers/fft.pdf
 * Description: fft(a) computes $\hat f(k) = \sum_x a[x] \exp(2\pi i \cdot k x / N)$ for all $k$. N must be a power of 2.
   Useful for convolution:
   \texttt{conv(a, b) = c}, where $c[x] = \sum a[i]b[x-i]$.
   For convolution of complex numbers or more than two vectors: FFT, multiply
   pointwise, divide by n, reverse(start+1, end), FFT back.
   Rounding is safe if $(\sum a_i^2 + \sum b_i^2)\log_2{N} < 9\cdot10^{14}$
   (in practice $10^{16}$; higher for random inputs).
   Otherwise, use NTT/FFTMod.
 * Time: O(N \log N) with $N = |A|+|B|$ ($\tilde 1s$ for $N=2^{22}$)
 * Status: somewhat tested
 * Details: An in-depth examination of precision for both FFT and FFTMod can be found
 * here (https://github.com/simonlindholm/fft-precision/blob/master/fft-precision.md)
 */
//#pragma once

typedef complex<double> C;
typedef vector<double> vd;
void fft(vector<C>& a) {
	int n = sz(a), L = 31 - __builtin_clz(n);
	static vector<complex<long double>> R(2, 1);
	static vector<C> rt(2, 1);  // (^ 10% faster if double)
	for (static int k = 2; k < n; k *= 2) {
		R.resize(n); rt.resize(n);
		auto x = polar(1.0L, acos(-1.0L) / k);
		rep(i,k,2*k) rt[i] = R[i] = i&1 ? R[i/2] * x : R[i/2];
	}
	vi rev(n);
	rep(i,0,n) rev[i] = (rev[i / 2] | (i & 1) << L) / 2;
	rep(i,0,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,0,k) {
			// C z = rt[j+k] * a[i+j+k]; // (25% faster if hand-rolled)  /// include-line
			auto x = (double *)&rt[j+k], y = (double *)&a[i+j+k];        /// exclude-line
			C z(x[0]*y[0] - x[1]*y[1], x[0]*y[1] + x[1]*y[0]);           /// exclude-line
			a[i + j + k] = a[i + j] - z;
			a[i + j] += z;
		}
}
vd multiply(vd a, vd b) {
	if (a.empty() || b.empty()) return {};
	vd res(sz(a) + sz(b) - 1);
	int L = 32 - __builtin_clz(sz(res)), n = 1 << L;
	vector<C> in(n), out(n);
	copy(all(a), begin(in));
	rep(i,0,sz(b)) in[i].imag(b[i]);
	fft(in);
	for (C& x : in) x *= x;
	rep(i,0,n) out[i] = in[-i & (n - 1)] - conj(in[i]);
	fft(out);
	rep(i,0,sz(res)) res[i] = imag(out[i]) / (4 * n);
	return res;
}

ld a,b,c,d;
ll n,K;
vector<ld>liczby;
vector<ld>prob[51007],pref[IP+4];

int main()
{
    ios_base::sync_with_stdio(0);cin.tie(0);
    cin>>n>>K;
    cout<<setprecision(20);
    for(ll i=0;i<n;i++){
        cin>>a;
        liczby.pb(a);
    }
    sort(liczby.begin(),liczby.end());
    reverse(liczby.begin(),liczby.end());
    n+=IP+4;
    for(ll i=0;i<IP+4;i++)liczby.pb(0);
    vector<ll>konce={0};
    for(ll i=0;i<IP;i++){
        konce.pb(((i+1)*n)/IP);
        prob[konce[i]+1]={(ld)1-liczby[konce[i]],0,liczby[konce[i]]};
        for(ll j=konce[i]+2;j<=konce.back();j++){
  prob[j].resize(prob[j-1].size() + 2, 0);
            for(ll k=0;k<prob[j-1].size();k++){
                prob[j][k]+=prob[j-1][k]*((ld)1-liczby[j-1]);
                prob[j][k+2]+=prob[j-1][k]*liczby[j-1];
            }

        }
    }
//cout<<"faza 1: "<<flush;
    pref[0]={1};
    for(int i=1;i<=IP;i++){
        pref[i]=multiply(pref[i-1],prob[konce[i]]);
        pref[i].resize(1 + 2 * konce[i]);
    }

    for(ll i=1;i<=IP;i++){
        for(ll j=pref[i].size()-2;j>=0;j--)
            pref[i][j]+=pref[i][j+1];
    }
    for(ll i=n;i>0;i--){
        prob[i]=prob[i-1];
    }
    prob[1]={1};
    ld wyn=0;
    konce[0]=-1;
    //cout<<"faza 2:   "<<flush;
    for(ll i=0;i<IP;i++){
        for(ll j=konce[i]+2;j<=konce[i+1]+1;j++){
            ld akwyn=0;
            for(ll k=max(0LL,(ll)((prob[j].size()-1)/2+K-(pref[i].size()-1)/2));k<prob[j].size() && (pref[i].size()-1)/2+max((ll)-(pref[i].size()-1)/2,(ll)(K-(k-(prob[j].size()-1)/2)))<pref[i].size();k++){
                akwyn+=prob[j][k]*pref[i][(pref[i].size()-1)/2+max((ll)-(pref[i].size()-1)/2,(ll)(K-(k-(prob[j].size()-1)/2)))];
            }
            wyn=max(wyn,akwyn);
        }
    }
    cout<<fixed<<wyn;

    return 0;
}