/**
 * code generated by JHelper
 * More info: https://github.com/AlexeyDmitriev/JHelper
 * @author cherrypi59
 */

#pragma GCC optimize ("O2")
#pragma GCC optimize ("tree-vectorize")
#pragma GCC target ("sse4")

#include <bits/stdc++.h>
#define FOR(i, a, b) for(int i=(a);i<(b);++i)
#define FFOR(i, a, b) FOR(i, a, b+1)
#define REP(i, n) FOR(i, 0, n)
#define RREP(i, n) FOR(i, 1, n+1)
#define ALL(v) (v).begin(), (v).end()
#define RALL(v) (v).rbegin(), (v).rend()
#define LEN(x) (int)(x).size()
#define DUMP(x) cerr<<__LINE__<<' '<<#x<<"="<<(x)<<endl;
#define pb push_back
#define eb emplace_back
#define mp make_pair

using namespace std;
using lint = long long;
using pii =  pair<int, int>;
using pll =  pair<lint, lint>;
template <typename T> using vc = vector<T>;
template <typename T> using vvc = vector<vector<T>>;
template <typename T> inline bool chmax(T &a, T b){if(a<b){a = b; return true;} return false;}
template <typename T> inline bool chmin(T &a, T b){if(a>b){a = b; return true;} return false;}

constexpr lint ten(int n) {return n==0 ? 1 : ten(n-1)*10;}

vector<complex<double>> dft(vector<complex<double>> a, bool inverse=false){
    // 数列aを離散フーリエ変換(DFT)、または逆離散フーリエ変換(IDFT)した数列を返す
    // inverse=false -> DFT, reverse=true -> IDFT
    // aは長さnが2冪の数列
    // O(n log n)
    const double PI = acos(-1);
    int n = a.size();
    int d = 0;
    for(int i=0;1<<i<n;++i) ++d;
    assert(n==1<<d);
    for(int i=0;i<n;++i){
        int rev = 0;
        for(int j=0;j<d;++j) rev |= (i>>j&1) << (d-j-1);
        if(i<rev) swap(a[i], a[rev]);
    }
    for(int i=1;i<n;i*=2){
        for(int j=0;j<i;++j){
            complex<double> w = polar(1.0, (2.0*PI/(2*i)*j)*(inverse ? 1 : -1));
            for(int k=0;k<n;k+=2*i){
                complex<double> first = a[j+k];
                complex<double> second = w * a[i+j+k];
                a[j+k] = first + second;
                a[i+j+k] = first - second;
            }
        }
    }
    if(inverse){
        for(int i=0;i<n;++i) a[i] /= n;
    }
    return a;
}

vector<complex<double>> square_dft(vector<complex<double>> a, bool inverse=false){
    const double PI = acos(-1);
    int n = a.size();
    vector<complex<double>> res(n);
    for(int i=0;i<n;++i){
        for(int j=0;j<n;++j){
            res[i] += a[j] * polar(1.0, (2.0*PI)/n*i*j*(inverse ? 1 : -1));
        }
    }
    return res;
}

bool f_equal(vector<complex<double>> a, vector<complex<double>> b){
    assert(a.size()==b.size());
    const double EPS = 1e-10;
    bool eq = true;
    for(int i=0;i<int(a.size());++i) eq &= fabs(a[i]-b[i]) < EPS;
    return eq;
}

template <typename T>
vector<double> convolution(vector<T> a, vector<T> b){
    int n = a.size() + b.size() - 1;
    int m = 1;
    while(m<n) m <<= 1;
    vector<complex<double>> a_complex(m), b_complex(m);
    for(int i=0;i<int(a.size());++i) a_complex[i] = complex<double>(a[i], 0);
    for(int i=0;i<int(b.size());++i) b_complex[i] = complex<double>(b[i], 0);
#ifdef DEBUG
    //vector<complex<double>> a_complex_dbg = a_complex;
    //vector<complex<double>> b_complex_dbg = b_complex;
#endif
    a_complex = dft(a_complex);
    b_complex = dft(b_complex);
#ifdef DEBUG
    //vector<complex<double>> a_square = square_dft(a_complex_dbg);
    //vector<complex<double>> b_square = square_dft(b_complex_dbg);
    //assert(f_equal(a_complex, a_square));
    //assert(f_equal(b_complex, b_square));
#endif
    for(int i=0;i<m;++i) a_complex[i] *= b_complex[i];
    vector<complex<double>> c = dft(a_complex, true);
    vector<double> res(n);
    for(int i=0;i<n;++i) res[i] = c[i].real();
    return res;
}

class EHandshake {
public:
    void solve(std::istream& in, std::ostream& out){
        ios_base::sync_with_stdio(false);
        in.tie(nullptr), out.tie(nullptr);
        const int MX = ten(5);
        int N;
        lint M;
        in>>N>>M;
        vc<int> cnt(MX+1, 0);
        REP(i, N){
            int a;
            in>>a;
            cnt[a]++;
        }
        lint ans = 0;
        vc<double> res = convolution(cnt, cnt);
        for(int i=2*MX;0<=i;--i){
            lint c = lint(res[i]+0.5);
            lint t = min(c, M);
            ans += i * t;
            M -= t;
            if(not M) break;
        }
        out<<ans<<'\n';
    }
};


int main() {
	EHandshake solver;
	std::istream& in(std::cin);
	std::ostream& out(std::cout);
	solver.solve(in, out);
	return 0;
}


Measured memory_bytes: 15480