Problem:
 Score : 500 points    Problem Statement Takahashi has come to a party as a special guest.
There are N ordinary guests at the party. The i -th ordinary guest has a power of A_i .  Takahashi has decided to perform M  handshakes to increase the happiness of the party (let the current happiness be 0 ).
A handshake will be performed as follows:   Takahashi chooses one (ordinary) guest x for his left hand and another guest y for his right hand ( x and y can be the same).  Then, he shakes the left hand of Guest x and the right hand of Guest y simultaneously to increase the happiness by A_x+A_y .   However, Takahashi should not perform the same handshake more than once. Formally, the following condition must hold:   Assume that, in the k -th handshake, Takahashi shakes the left hand of Guest x_k and the right hand of Guest y_k . Then, there is no pair p, q  (1 \leq p < q \leq M) such that (x_p,y_p)=(x_q,y_q) .   What is the maximum possible happiness after M handshakes?      Constraints  1 \leq N \leq 10^5  1 \leq M \leq N^2  1 \leq A_i \leq 10^5  All values in input are integers.         Input Input is given from Standard Input in the following format:  N  M  A_1  A_2  ...  A_N       Output Print the maximum possible happiness after M handshakes.        Sample Input 1 5 3
10 14 19 34 33      Sample Output 1 202  Let us say that Takahashi performs the following handshakes:   In the first handshake, Takahashi shakes the left hand of Guest 4 and the right hand of Guest 4 .  In the second handshake, Takahashi shakes the left hand of Guest 4 and the right hand of Guest 5 .  In the third handshake, Takahashi shakes the left hand of Guest 5 and the right hand of Guest 4 .   Then, we will have the happiness of (34+34)+(34+33)+(33+34)=202 .  We cannot achieve the happiness of 203 or greater, so the answer is 202 .       Sample Input 2 9 14
1 3 5 110 24 21 34 5 3      Sample Output 2 1837       Sample Input 3 9 73
67597 52981 5828 66249 75177 64141 40773 79105 16076      Sample Output 3 8128170

Solution (Python):

import sys
input = sys.stdin.readline
import bisect

n, m = map(int, input().split())
a = [int(item) for item in input().split()]
a.sort()
a_rev = sorted(a, reverse=True)

cumsum = [0]
for item in a:
    cumsum.append(cumsum[-1] + item)


l = 0; r = 10**10
while r - l > 1:
    mid = (l + r) // 2
    to_use = 0
    for i, item in enumerate(a_rev):
        useless = bisect.bisect_left(a, mid - item)
        to_use += n - useless
    if to_use >= m:
        l = mid
    else:
        r = mid
ans = 0
total_use = 0
for i, item in enumerate(a_rev):
    useless = bisect.bisect_left(a, l - item)
    to_use = n - useless
    total_use += to_use
    ans += item * to_use + cumsum[n] - cumsum[n - to_use]
print(ans - l * (total_use - m))

Solution (Python):
def modinv(x, mod):
    a,b = x,mod
    u,v = 1,0
    while b:
        t = a//b
        a -= t*b; a,b = b,a
        u -= t*v; u,v = v,u
    return u%mod

def _garner(xs, mods):
    M = len(xs)
    coeffs = [1]*M
    constants = [0]*M
    for i in range(M-1):
        mod_i = mods[i]
        v = (xs[i] - constants[i]) * modinv(coeffs[i], mod_i) % mod_i
        for j in range(i+1,M):
            mod_j = mods[j]
            constants[j] = (constants[j] + coeffs[j] * v) % mod_j
            coeffs[j] = (coeffs[j] * mod_i) % mod_j

    return constants[-1]


def bit_reverse(d):
    n = len(d)
    ns = n>>1; nss = ns>>1
    ns1 = ns + 1
    i = 0
    for j in range(0, ns, 2):
        if j<i:
            d[i], d[j] = d[j], d[i]
            d[i+ns1], d[j+ns1] = d[j+ns1], d[i+ns1]
        d[i+1], d[j+ns] = d[j+ns], d[i+1]
        k = nss; i ^= k
        while k > i:
            k >>= 1; i ^= k
    return d

class NTT:
    def __init__(self, mod, primitive_root):
        self.mod = mod
        self.root = primitive_root

    def _ntt(self, a, sign):
        n = len(a)
        mod, g = self.mod, self.root
        tmp = (mod - 1) * modinv(n, mod) % mod 
        h = pow(g, tmp, mod) 
        if sign < 0: h = modinv(h, mod)

        a = bit_reverse(a)

        m = 1
        while m < n:
            m2 = m<<1
            _base = pow(h, n//m2, mod)
            _w = 1
            for x in range(m):
                for s in range(x,n,m2):
                    u = a[s]
                    d = a[s+m]* _w % mod
                    a[s] = (u+d)%mod
                    a[s+m] = (u-d)%mod
                _w = _w * _base % mod
            m <<= 1
        return a

    def ntt(self, a):
        return self._ntt(a,1)

    def intt(self, a):
        mod = self.mod
        n = len(a)
        n_inv = modinv(n, mod)
        a = self._ntt(a,-1)
        for i in range(n):
            a[i] = a[i]*n_inv%mod
        return a

    def convolution(self, a, b):
        mod = self.mod
        ret_size = len(a)+len(b)-1
        n = 1 << (ret_size-1).bit_length()
        _a = a + [0] * (n - len(a))
        _b = b + [0] * (n - len(b))
        _a = self.ntt(_a)
        _b = self.ntt(_b)
        _a = [x*y%mod for x,y in zip(_a,_b)]
        _a = self.intt(_a)
        _a = _a[:ret_size]
        return _a


def convolution_ntt(a, b, mod):
    a = [x%mod for x in a]
    b = [x%mod for x in b]

    mods = (167772161, 469762049, 1224736769, mod)

    ntt1 = NTT(mods[0], 3)
    ntt2 = NTT(mods[1], 3)
    ntt3 = NTT(mods[2], 3)

    x1 = ntt1.convolution(a, b)
    x2 = ntt2.convolution(a, b)
    x3 = ntt3.convolution(a, b)

    n = len(x1)
    ret = [0]*n
    for i in range(n):
        xs = [x1[i], x2[i], x3[i], 0]
        ret[i] = _garner(xs, mods)

    return ret

m,k = map(int,input().split())
f = list(map(int,input().split()))
n = max(f)+1
g = [0]*n
for x in f: g[x] += 1
h = convolution_ntt(g,g,10**11)
ans = 0
cnt = 0
for i in range(len(h)-1,-1,-1):
    if cnt + h[i] < k:
        ans += h[i]*i
        cnt += h[i]
    else:
        ans += (k-cnt)*i
        break
print(ans)

Solution (Python):
import sys
import bisect
input = sys.stdin.readline
import numpy as np

class Bisect(object):
    def bisect_max(self, reva, func,M):
        ok = 0 
        ng = 2*(10**5)+1 
        while abs(ok-ng) > 1:
            cnt = (ok + ng) // 2
            if func(cnt,reva,M):
                ok = cnt
            else:
                ng = cnt
        return ok

def solve1(tgt,reva,M):
    res=0
    n = len(reva)
    if (n - bisect.bisect_left(reva,tgt-reva[n-1]))*n < M:
        return False
    for i in range(n):
        res += n - bisect.bisect_left(reva,tgt-reva[n-i-1])
        if M <= res:
            return True
    return False

N,M = map(int,input().split())
a = list(map(int,input().split()))
a.sort(reverse=True)
reva = a[::-1]
bs = Bisect()
Kmax = (bs.bisect_max(reva,solve1,M))

r=np.zeros((N+1),dtype='int64')
np.cumsum(a,out=r[1:])
r = list(r)
res = 0
cnt = 0
t = 0
for i in range(N):
    tmp = bisect.bisect_left(reva,Kmax-reva[N-i-1])
    if t != 1:
        tmp2 = bisect.bisect_right(reva,Kmax-reva[N-i-1])
        if tmp!=tmp2:
            t = 1
    tmp = N - tmp
    cnt += tmp
    res += tmp*a[i]+r[tmp]
if t==1:
    res -= (cnt-M)*Kmax
print(res)

Solution (Python):
import sys
import bisect
input = sys.stdin.readline
import numpy as np

class Bisect(object):
    def bisect_max(self, reva, func,M):
        ok = 0 
        ng = 2*(10**5)+1 
        while abs(ok-ng) > 1:
            cnt = (ok + ng) // 2
            if func(cnt,reva,M):
                ok = cnt
            else:
                ng = cnt
        return ok

def solve1(tgt,reva,M):
    res=0
    n = len(reva)
    if (n - bisect.bisect_left(reva,tgt-reva[n-1]))*n < M:
        return False
    for i in range(n):
        res += n - bisect.bisect_left(reva,tgt-reva[n-i-1])
        if M <= res:
            return True
    return False

N,M = map(int,input().split())
a = list(map(int,input().split()))
a.sort(reverse=True)
reva = a[::-1]
bs = Bisect()
Kmax = (bs.bisect_max(reva,solve1,M))

r=np.zeros((N+1),dtype='int64')
np.cumsum(a,out=r[1:])

res = 0
cnt = 0
t = 0
for i in range(N):
    tmp = bisect.bisect_left(reva,Kmax-reva[N-i-1])
    tmp2 = bisect.bisect_right(reva,Kmax-reva[N-i-1])
    if tmp!=tmp2:
        t = 1
    tmp = N - tmp
    cnt += tmp
    res += tmp*a[i]+r[tmp]
if t==1:
    res -= (cnt-M)*Kmax
print(res)

Solution (Python):
def main():
  import bisect
  import sys
  from itertools import accumulate
  input = sys.stdin.readline
  def hsnum(x):
    num = 0
    for i in a:
      y = x-i
      num += n-bisect.bisect_left(a,y)
    return num
  n,m = map(int, input().split())
  a = list(map(int, input().split()))
  a.sort()
  suma = [0]+list(accumulate(a))
  lx = 0
  rx = 2*10**5+2
  error = 0
  while lx+1 < rx:
    mx = (lx+rx)//2
    if hsnum(mx) < m:
      rx = mx
    elif hsnum(mx) > m:
      lx = mx + 1
    else:
      lx = mx
      break
  if hsnum(lx)<m:
    lx -= 1
  error = hsnum(lx)-m
  ans = 0
  for i in a:
    rema = lx - i
    j = bisect.bisect_left(a,rema)
    ans += (n-j)*i-suma[j]
  ans += suma[n]*n
  print(ans-error*lx)
if __name__ == "__main__":
    main()

Solution (Python):
def main():
  import bisect
  import sys
  from itertools import accumulate
  input = sys.stdin.readline
  def hsnum(x):
    num = 0
    for i in a:
      y = x-i
      num += n-bisect.bisect_left(a,y)
    return num
  n,m = map(int, input().split())
  a = list(map(int, input().split()))
  a.sort()
  suma = [0]+list(accumulate(a))
  lx = 0
  rx = 2*10**5+2
  error = 0
  while lx+1 < rx:
    mx = (lx+rx)//2
    if hsnum(mx) < m:
      rx = mx
    elif hsnum(mx) > m:
      lx = mx + 1
    else:
      lx = mx
      break
  if hsnum(lx)<m:
    lx -= 1
  error = hsnum(lx)-m
  ans = 0
  for i in a:
    rema = lx - i
    j = bisect.bisect_left(a,rema)
    ans += (n-j)*i-suma[j]
  ans += suma[n]*n
  print(ans-error*lx)
if __name__ == "__main__":
    main()

Solution (Python):
from bisect import bisect_left

N, M=map(int, input().split())
A=list(map(int, input().split()))

A.sort()
left=0
right=max(A) * 2 + 1

while left+1<right:
  mid=(left + right) // 2
  
  cnt=0
  for i in range(N):
    cnt+=N - bisect_left(A, mid - A[i])
  if cnt>=M:
    left=mid
  else:
    right=mid

mon=[0]
for a in A:
  mon.append(mon[-1] + a)

ans=0
cnt=0
for i in range(N):
  bi = bisect_left(A, left - A[i])
  ans+=(mon[-1] - mon[bi]) + (N - bi) * A[i]
  cnt+=(N - bi)
ans-=(cnt - M) * left

print(ans)

Solution (Python):
n, m = map(int, input().split())
l = list(map(int, input().split()))
l.sort()

import bisect 

def func(x):
  C = 0
  for p in l:
    q = x -p
    j = bisect.bisect_left(l, q)
    C += n-j
  if C >= m:
    return True
  else:
    return False
  
l_ = 0
r_ = 2*10**5 +1
while l_+1 < r_:
  c_ = (l_+r_)//2
  if func(c_):
    l_ = c_
  else:
    r_ = c_
    
ans = 0
cnt = 0
lr = sorted(l, reverse=True)
from itertools import accumulate
cum = [0] + list(accumulate(lr))
for i in lr:
  j = bisect.bisect_left(l, l_-i)
  ans += i*(n-j) + cum[n-j]
  cnt += n -j
ans -= (cnt-m)*l_
print(ans)

Solution (Python):
n, m = map(int, input().split())
l = list(map(int, input().split()))
l.sort()

import bisect 

def func(x):
  C = 0
  for p in l:
    q = x -p
    j = bisect.bisect_left(l, q)
    C += n-j
  if C >= m:
    return True
  else:
    return False
  
l_ = 0
r_ = 2*10**5 +1
while l_+1 < r_:
  c_ = (l_+r_)//2
  if func(c_):
    l_ = c_
  else:
    r_ = c_
    
ans = 0
cnt = 0
lr = sorted(l, reverse=True)
from itertools import accumulate
cum = [0] + list(accumulate(lr))
for i in lr:
  j = bisect.bisect_left(l, l_-i)
  ans += i*(n-j) + cum[n-j]
  cnt += n -j
ans -= (cnt-m)*l_
print(ans)

Solution (Python):
import bisect
from itertools import accumulate

def func(x):
    cnt = 0
    for i in a:
        index = bisect.bisect_left(a, x - i)
        cnt += N - index
    if cnt >= M:
        return True
    else:
        return False

N, M = map(int, input().split())
a = list(map(int, input().split()))
a.sort()
ar = sorted(a, reverse=True)
b = [0] + list(accumulate(ar))

MIN = 0
MAX = 2 * 10 ** 5 + 1
while MAX - MIN > 1:
    MID = (MIN + MAX) // 2
    if func(MID):
        MIN = MID
    else:
        MAX = MID

ans = 0
cnt = 0
for i in ar:
    index = bisect.bisect_left(a, MIN - i)
    ans += i * (N - index) + b[N - index]
    cnt += N - index
print(ans - (cnt - M) * MIN)

Solution (Python):
N, M = map(int, input().split())
A = list(map(int, input().split()))

A.sort(reverse=True)
LB = A[-1] * 2  
UB = A[0] * 2   


cgs = [0] * (UB + 1)
for a in A:
    cgs[a] += 1
for i in range(UB, 0, -1):
    cgs[i - 1] += cgs[i]



def is_ok(n):
    m = 0
    for a in A:
        m += cgs[max(n - a, 0)]
    return m >= M


ok = LB
ng = UB + 1
while ng - ok > 1:
    m = (ok + ng) // 2
    if is_ok(m):
        ok = m
    else:
        ng = m


cps = A[:]
for i in range(1, N):
    cps[i] += cps[i - 1]

result = 0
for a in A:

    t = cgs[max(ok - a, 0)] 
    if t != 0:
        result += a * t + cps[t - 1] 
        M -= t
result += ok * M 
print(result)

Solution (Python):
from sys import stdin
f_i = stdin

N, M = map(int, f_i.readline().split())
A = list(map(int, f_i.readline().split()))

A.sort()

from bisect import bisect_left

def shake_count(power):
    return N ** 2 - sum(bisect_left(A, power - left_hand) for left_hand in A)

lower = A[0] * 2 - 1
upper = A[-1] * 2 + 1
while lower + 1 < upper:
    mid = (lower + upper) // 2
    if shake_count(mid) >= M:
        lower = mid
    else:
        upper = mid

X = tuple(bisect_left(A, upper - left_hand) for left_hand in A)

from itertools import accumulate
acc = tuple(accumulate([0] + A))

s = acc[-1]
ans = sum(s - acc[x] for x in X) + sum(a * (N - x) for a, x in zip(A, X))

rest_shake_count = M - (N ** 2 - sum(X))
ans += lower * rest_shake_count

print(ans)

Solution (Python):
from bisect import bisect_left
from itertools import accumulate


n, m = map(int, input().split())
a = list(map(int, input().split()))
a.sort()

b = [0] + a
b.sort(reverse=True)
b = [0] +b
b = list(accumulate(b))

def isok(sum_):
    n_pattern = 0
    for i in range(n):
        n_pattern += n - bisect_left(a, sum_ - a[i])
        if n_pattern >= m:
            return True

    return False

ok = min(a)*2
ng = max(a)*2 +1

while ng - ok > 1:
    mid = (ok+ng)//2
    if isok(mid):
        ok = mid
    else:
        ng = mid

sum_ng = ng
sum_ok = ok

m_counter_ng = 0

ans = 0
for i in range(n):
    id_ng = bisect_left(a, sum_ng - a[i])
    ans += b[n-id_ng]
    ans += a[i]*(n - id_ng)
    m_counter_ng += n - id_ng

ans += sum_ok * (m - m_counter_ng)

print(ans)

Solution (Python):
from bisect import bisect_left


def calc(N, A, cumsum, x):
    num = N ** 2
    total = cumsum[N] * N * 2
    for a in A:
        idx = bisect_left(A, x - a)
        num -= idx
        total -= cumsum[idx] + idx * a
    return num, total


def main():
    N, M = map(int, input().split())
    A = list(map(int, input().split()))
    A.sort()
    cumsum = [0] * (N+1)
    for i in range(N):
        cumsum[i+1] = cumsum[i] + A[i]
    lo = 0
    hi = int(2e5+1)
    while lo < hi:
        mid = (lo + hi) // 2
        num, total = calc(N, A, cumsum, mid)
        if num <= M:
            hi = mid
        else:
            lo = mid + 1
    num, total = calc(N, A, cumsum, lo)
    ans = total - (num - M) * (lo - 1)
    print(ans)


if __name__ == "__main__":
    main()

Solution (Python):
import sys
import math
from collections import defaultdict

sys.setrecursionlimit(10**7)
def input():
    return sys.stdin.readline()[:-1]

mod = 10**9 + 7

def I(): return int(input())
def II(): return map(int, input().split())
def III(): return list(map(int, input().split()))
def Line(N,num):
    if N<=0:
        return [[] for _ in range(num)]
    elif num==1:
        return [I() for _ in range(N)]
    else:
        read_all = [tuple(II()) for _ in range(N)]
        return map(list, zip(*read_all))



import cmath
pi = cmath.pi
exp = cmath.exp

def convolution(a,b):
    
    def fft(a,sz,inv=False):
        tmp = [0]*sz
        mask = sz-1
        p = 0
        i = sz>>1
        sign = 1 if inv else -1

        while i:
            if p&1:
                cur,nex = tmp,a
            else:
                cur,nex = a,tmp
            ei = exp(2j*pi*i*sign/sz)
            w = 1
            for j in range(0,sz,i):
                for k in range(i):
                    nex[j+k] = cur[((j<<1)&mask)+k] + w*cur[(((j<<1)+i)&mask)+k]
                w *= ei
            p += 1
            i >>= 1

        if p&1:
            a,tmp = tmp,a
        if inv:
            a = list(map(lambda x: x/sz, a))
        return a
    
    sz = 1<<(len(a)+len(b)-2).bit_length()
    a = a + [0]*(sz-len(a))
    b = b + [0]*(sz-len(b))

    fa = fft(a,sz)
    fb = fft(b,sz)
    fafb = [fai*fbi for fai,fbi in zip(fa,fb)]
    ab = fft(fafb,sz,inv=True)
    return [round(x.real) for x in ab[:len(a)+len(b)-1]]

N,M = II()
A = III()

h = [0]*(max(A))
for a in A:
    h[a-1] += 1

conv = [0,0] + convolution(h,h)

ans = 0
for k in range(2,2*max(A)+1)[::-1]:
    if conv[k]<M:
        ans += k*conv[k]
        M -= conv[k]
    else:
        ans += k*M
        break
        
print(ans)

Solution (Python):
import sys
import math
from collections import defaultdict

sys.setrecursionlimit(10**7)
def input():
    return sys.stdin.readline()[:-1]

mod = 10**9 + 7

def I(): return int(input())
def II(): return map(int, input().split())
def III(): return list(map(int, input().split()))
def Line(N,num):
    if N<=0:
        return [[] for _ in range(num)]
    elif num==1:
        return [I() for _ in range(N)]
    else:
        read_all = [tuple(II()) for _ in range(N)]
        return map(list, zip(*read_all))







import numpy as np
def Convolution(a,b):
    bit = (len(a)+len(b)-2).bit_length()
    L = 2**bit
    fa,fb = np.fft.rfft(a,L), np.fft.rfft(b,L)
    c = np.rint(np.fft.irfft(fa*fb,L)).astype(np.int64)
    return c[:len(a)+len(b)-1]

N,M = II()
A = III()

h = [0]*(max(A))
for a in A:
    h[a-1] += 1

conv = np.append([0,0],Convolution(h,h))
ans = 0
for k in range(2,2*max(A)+1)[::-1]:
    if conv[k]<M:
        ans += k*conv[k]
        M -= conv[k]
    else:
        ans += k*M
        break
        
print(ans)

Solution (Python):
import sys
import math
from collections import defaultdict

sys.setrecursionlimit(10**7)
def input():
    return sys.stdin.readline()[:-1]

mod = 10**9 + 7

def I(): return int(input())
def II(): return map(int, input().split())
def III(): return list(map(int, input().split()))
def Line(N,num):
    if N<=0:
        return [[] for _ in range(num)]
    elif num==1:
        return [I() for _ in range(N)]
    else:
        read_all = [tuple(II()) for _ in range(N)]
        return map(list, zip(*read_all))







import numpy as np
def FFT(a,b):
    bit = (len(a)+len(b)).bit_length()+1
    L = 2**bit
    fa,fb = np.fft.rfft(a,L), np.fft.rfft(b,L)
    c = np.rint(np.fft.irfft(fa*fb,L)).astype(np.int64)
    return c[:len(a)+len(b)-1]

N,M = II()
A = III()

h = [0]*(max(A))
for a in A:
    h[a-1] += 1

conv = np.append([0,0],FFT(h,h))
ans = 0
count = 0

for k in range(2,2*max(A)+1)[::-1]:
    if conv[k]:
        num = min(M-count,conv[k])
        count += num
        ans += k*num
    if count==M:
        break

print(ans)

Solution (Python):
N, M = map(int, input().split())
A = list(map(int, input().split()))

A.sort(reverse=True)
LB = A[-1] * 2  
UB = A[0] * 2   


cgs = [0] * (UB + 1)
for a in A:
    cgs[a] += 1
for i in range(UB, 0, -1):
    cgs[i - 1] += cgs[i]



def is_ok(n):
    m = 0
    for a in A:
        m += cgs[max(n - a, 0)]
    return m >= M


ok = LB
ng = UB + 1
while ng - ok > 1:
    m = (ok + ng) // 2
    if is_ok(m):
        ok = m
    else:
        ng = m


cps = A[:]
for i in range(1, N):
    cps[i] += cps[i - 1]

result = 0
for a in A:
    t = cgs[max(ok - a, 0)]
    if t != 0:
        result += a * t + cps[t - 1]
        M -= t
result += ok * M
print(result)

Solution (Python):
from itertools import accumulate
n,m=map(int,input().split())
a=list(map(int,input().split()))
a.sort(reverse=True)
b=[0]*(2*a[0]+1)
for i in range(n):
  b[a[i]]+=1
b=list(accumulate(b[::-1]))[::-1]
s=list(accumulate([0]+a))
ok=0
ng=2*a[0]+1
while ng-ok>1:
  mid=(ng+ok)//2
  cnt=0
  for i in range(n):
    cnt+=b[max(0,mid-a[i])]
  if cnt>=m:
    ok=mid
  else:
    ng=mid
cnt=0
ans=0
for i in range(n):
  t=min(b[max(0,ng-a[i])],m-cnt)
  ans+=s[t]+a[i]*t
  cnt+=t
ans+=(m-cnt)*ok
print(ans)

Solution (Python):

N, M = map(int, list(input().split()))
A = list(map(int, input().split()))


A.sort(reverse=True)


def calc(x):
    count, sum = 0, 0

    j, t = 0, 0
    for i in reversed(range(N)):
        while j < N and A[i]+A[j] >= x:
            t += A[j]
            j += 1
        count += j
        sum += A[i]*j + t
    
    return (count, sum)


def binary_search(x, y):
    mid = (x+y)//2

    count, sum = calc(mid)
    
    if count < M:
        return binary_search(x, mid)
    else:
        if x == mid:
            print(sum-(count-M)*mid)
        else:
            return binary_search(mid, y)


binary_search(0, A[0]*2+1)

Solution (Python):
from bisect import bisect_left

N, M = map(int, input().split())
A = list(map(int, input().split()))

A.sort()
left = 0
right = max(A) * 2 + 1

while left + 1 < right:
  mid = (left + right) // 2
  
  cnt = 0
  for i in range(N):
    cnt += N - bisect_left(A, mid - A[i])
  if cnt >= M:
    left = mid
  else:
    right = mid

Acum = [0]
for a in A:
  Acum.append(Acum[-1] + a)

answer = 0
cnt = 0
for i in range(N):
  bi = bisect_left(A, left - A[i])
  answer += (Acum[-1] - Acum[bi]) + (N - bi) * A[i]
  cnt += (N - bi)
answer -= (cnt - M) * left

print(answer)

Solution (Python):
from bisect import bisect_left

n, m = map(int, input().split())
a = list(map(int, input().split()))

a.sort(reverse=True)
b = list(reversed(a))

rs = [0] * n
rs[0] = a[0]
for i, x in enumerate(a[1:]):
    rs[i+1] = rs[i] + x

L = a[-1] * 2
H = a[0] * 2

while H - L > 1:
    M = (L + H) // 2
    count = 0
    for x in a:
       count += n - bisect_left(b, M - x)
    if count == m:
        L = M
        break
    elif count > m:
        L = M
    else:
        H = M

ans = 0
count = 0
for x in a:
    c = n - 1 -  bisect_left(b, L - x)
    count += c + 1
    if c >= 0:
        ans += rs[c] + x * (c+1)

print(ans - (count-m) * L)

Solution (Python):
from bisect import bisect_left

N,M= [int(i) for i in input().split()]
n = sorted([int(x) for x in input().split()])

s =[0]*(N+1)
s[0]=n[0]
for i in range(1,N+1):
        s[i] = s[i-1]+n[i-1]
MM=n[-1]*2
mm=n[0]*2

while True:
        c=0        
        r=(mm+MM)//2
        for i in range(N):
                c += N-bisect_left(n,r-n[i])
        if (MM-mm) < 2:
                break
        elif c>=M:
                mm=r
        elif c<M:
                MM=r

v=0
for i in range(N):
        j=bisect_left(n,r-n[i])
        v+= (s[N] - s[j]) + (N-j)*n[i]

v = v - (c-M)*(r)

print(v)

Solution (Python):
import numpy

n, m = [int(i) for i in input().split()]
A = [int(i) for i in input().split()]
sup_A = 10 ** 5
B = [0] * (1 << 18)
for a in A:
    B[a] += 1
C = numpy.fft.fft(B)
D = numpy.fft.ifft(C * C)

ans = 0
cnt = 0
for i in range(2 * sup_A, -1, -1):
    d = int(D[i].real + 0.5)
    if cnt + d >= m:
        ans += (m - cnt) * i
        break
    cnt += d
    ans += i * d

print(ans)

Solution (Python):
import numpy

n, m = [int(i) for i in input().split()]
A = [int(i) for i in input().split()]
sup_A = 10 ** 5
B = [0] * (2 * sup_A + 1)
for a in A:
    B[a] += 1
C = numpy.fft.fft(B)
D = numpy.fft.ifft(C * C)

ans = 0
cnt = 0
for i in range(2 * sup_A, -1, -1):
    d = int(D[i].real + 0.5)
    if cnt + d >= m:
        ans += (m - cnt) * i
        break
    cnt += d
    ans += i * d

print(ans)

Solution (Python):
import sys
input = sys.stdin.readline
import numpy as np
N, M = map(int, input().split())
L = list(map(int, input().split()))
A = np.zeros(1<<18)
for i in L:
  A[i] += 1
C = (np.fft.irfft(np.fft.rfft(A) * np.fft.rfft(A)) + .5).astype(np.int64)
G = C.cumsum()
count = np.searchsorted(G, N*N-M)
rest = N*N-M - G[count-1]
temp = (C[:count]*np.arange(count, dtype=np.int64)).sum() + count*rest
ans = sum(L)*2*N - temp
print(ans)

Solution (Python):
import bisect
N, M = map(int, input().split())
A = [int(x) for x in input().split()]
A.sort()
INF = 10 ** 9
A.append(INF)

left, right = 0, 2 * max(A)
while left + 1 < right:
    mid = (left + right) >> 1
    cnt = sum(N - bisect.bisect_left(A, mid-A[i]) for i in range(N))
    if cnt >= M:
        left = mid
    else:
        right = mid
acc = [0]
for i in range(N):
    acc.append(acc[-1] + A[i])
ans = 0
cnt = 0
for i in range(N):
    idx = bisect.bisect_left(A, left-A[i])
    cnt += N - idx
    ans += (N - idx) * A[i] + acc[N] - acc[idx]
ans -= (cnt - M) * left
print(ans)

Solution (Python):
import bisect
N, M = map(int, input().split())
A = [int(x) for x in input().split()]
A.sort()
INF = 10 ** 9
A.append(INF)

left, right = 0, 2 * max(A)
while left + 1 < right:
    mid = (left + right) >> 1
    cnt = sum(N - bisect.bisect_left(A, mid-A[i]) for i in range(N))
    if cnt >= M:
        left = mid
    else:
        right = mid
acc = [0]
for i in range(N):
    acc.append(acc[-1] + A[i])
ans = 0
cnt, min_val = 0, INF
for i in range(N):
    idx = bisect.bisect_left(A, left-A[i])
    cnt += N - idx
    min_val = min(min_val, A[i] + A[idx])
    ans += (N - idx) * A[i] + acc[N] - acc[idx]
ans -= (cnt - M) * min_val
print(ans)

Solution (Python):
import numpy as np

n,m = [int(i) for i in input().split()]
a = [int(i) for i in input().split()]
d = 2**18

f = np.zeros(d,dtype=int)
for i in a:
    f[i]+=1

tf = np.fft.rfft(f)
f = np.fft.irfft(tf*tf)
f = [int(i+0.5) for i in f]

ans=0
for i in range(len(f)-1,0,-1):
    if f[i]<=m:
        ans+=i*f[i]
        m-=f[i]
    elif f[i]>m:
        ans+=i*m
        break
print(ans)

Solution (Python):
import numpy as np

n,m = [int(i) for i in input().split()]
a = [int(i) for i in input().split()]
d = 2**18

f = np.zeros(d,dtype=int)
for i in a:
    f[i]+=1

tf = np.fft.fft(f)
f = np.fft.ifft(tf*tf)
f = [int(i+0.5) for i in f]

ans=0
for i in range(len(f)-1,0,-1):
    if f[i]<=m:
        ans+=i*f[i]
        m-=f[i]
    elif f[i]>m:
        ans+=i*m
        break
print(ans)

Solution (Python):
import numpy as np

n,m = [int(i) for i in input().split()]
a = [int(i) for i in input().split()]
d = 2**18

f = np.array([0]*d)
for i in a:
    f[i]+=1

tf = np.fft.fft(f)
f = np.fft.ifft(tf*tf)
f = [int(i+0.5) for i in f]

ans=0
for i in range(len(f)-1,0,-1):
    if f[i]<=m:
        ans+=i*f[i]
        m-=f[i]
    elif f[i]>m:
        ans+=i*m
        break
print(ans)

Solution (Python):
import bisect

N,M = map(int,input().split())
A = sorted(list(map(int,input().split())))



csum = [A[0] for _ in range(N+1)]
for i in range(1,N):
  csum[i] = csum[i-1]+A[i]




l = 0
r = 10**6
while (l+1 < r):
  m = (l+r)//2
  p = 0 
  for a in A:
    q = bisect.bisect_left(A,m-a)
    
    p += N-q
  
  if p >= M:
    l = m
  else:
    r = m





ans = 0
cnt = 0
for a in A:
  q = bisect.bisect_right(A,l-1-a)
  if q == 0:
    ans += csum[N-1]+a*N
  elif q < N:
    ans += csum[N-1]-csum[q-1]+a*(N-q)
  cnt += N-q
  
ans += l*(M-cnt)
print(ans)

Solution (Python):
import bisect

N,M = map(int,input().split())
A = sorted(list(map(int,input().split())))



csum = [A[0] for _ in range(N+1)]
for i in range(1,N):
  csum[i] = csum[i-1]+A[i]



l = 0
r = 10**6
while (l+1 < r):
  m = (l+r)//2
  p = 0 
  for a in A:
    q = bisect.bisect_left(A,m-a)
    
    p += N-q
  
  if p >= M:
    l = m
  else:
    r = m





ans = 0
cnt = 0
for a in A:
  q = bisect.bisect_right(A,l-1-a)
  if q == 0:
    ans += csum[N-1]+a*N
    cnt += N
  elif q < N:
    ans += csum[N-1]-csum[q-1]+a*(N-q)
    cnt += N-q
  
ans +=l*(M-cnt)
print(ans)

Solution (Python):
import bisect

N,M = map(int,input().split())
A = sorted(list(map(int,input().split())))



csum = [A[0] for _ in range(N+1)]
for i in range(1,N):
  csum[i] = csum[i-1]+A[i]




l = 0
r = 10**6
while (l+1 < r):
  m = (l+r)//2
  p = 0 
  for a in A:
    q = bisect.bisect_left(A,m-a)
    
    p += N-q
  
  if p >= M:
    l = m
  else:
    r = m





ans = 0
cnt = 0
for a in A:
  q = bisect.bisect_right(A,l-a)
  if q == 0:
    ans += csum[N-1]+a*N
    cnt += N
  elif q < N:
    ans += csum[N-1]-csum[q-1]+a*(N-q)
    cnt += N-q
  
ans +=l*(M-cnt)
print(ans)

Solution (Python):
from collections import defaultdict, deque, Counter
from heapq import heappush, heappop, heapify
from bisect import bisect_right, bisect_left
import random
from itertools import permutations, accumulate, combinations, product
import sys
import string
from bisect import bisect_left, bisect_right
from math import factorial, ceil, floor, gamma, log
from operator import mul
from functools import reduce


sys.setrecursionlimit(2147483647)
INF = 1 << 100
def LI(): return list(map(int, sys.stdin.buffer.readline().split()))
def I(): return int(sys.stdin.buffer.readline())
def LS(): return sys.stdin.buffer.readline().rstrip().decode('utf-8').split()
def S(): return sys.stdin.buffer.readline().rstrip().decode('utf-8')
def IR(n): return [I() for i in range(n)]
def LIR(n): return [LI() for i in range(n)]
def SR(n): return [S() for i in range(n)]
def LSR(n): return [LS() for i in range(n)]
def SRL(n): return [list(S()) for i in range(n)]
def MSRL(n): return [[int(j) for j in list(S())] for i in range(n)]
mod = 1000000007


n, m = LI()
A = sorted(LI(), reverse=True)
max_num = 10 ** 5 * 2 + 1
D = [0] * max_num
acc = [0] + A
for i in range(1, n):
    acc[i + 1] += acc[i]

for i in A:
    D[i] += 1

for i in range(max_num - 1, 0, -1):
    D[i - 1] += D[i]

ok = 1
ng = max(A) * 2 + 1
while ng > ok + 1:
    mid = (ng + ok) // 2
    ret = 0
    happiness = 0
    for i in range(n):
        d = D[max(mid - A[i], 0)]
        ret += d
    if ret >= m:
        ok = mid
    else:
        ng = mid


ret2 = 0
ans = 0
for i in range(n):
    d2 = D[max(ng - A[i], 0)]
    ret2 += d2
    ans = ans + d2 * A[i] + acc[d2]

ans += (m - ret2) * ok
print(ans)

Solution (Python):
import sys
input = sys.stdin.readline

from bisect import bisect_left

def isOver(A,n,x):
	cnt = 0
	for a in A:
		cnt += n-bisect_left(A,x-a)
	return cnt

def calOut(A,s,n,x):
	out = 0
	for a in A:
		idx = bisect_left(A,x-a)
		out += s[n]-s[idx]+a*(n-idx)
	return out - (isOver(A,n,x)-m)*x

def calX(A,n,m,l,r):
	while r-l>1:
		
		x = (l+r)//2
		if isOver(A,n,x) >= m:
			l = x
		else: 
			r = x 
	return l

n,m = map(int, input().split())
A = list(map(int, input().split()))


A.sort()
S = [0]*(n+1)
for i in range(n):
	S[i+1] = S[i]+A[i]

x = calX(A,n,m,A[0]*2,A[n-1]*2)

print(calOut(A,S,n,x))

Solution (Python):
import sys
input = sys.stdin.readline

from bisect import bisect_left

def isOver(a,n,x):
	cnt = 0
	for i in range(n):
		cnt += n-bisect_left(a,x-a[i])
	return cnt

def calX(a,n,m,l,r):
	while r-l>1:
		
		x = (l+r)//2
		if isOver(a,n,x) >= m:
			l = x
		else: 
			r = x 
	return l

n,m = map(int, input().split())
a = list(map(int, input().split()))


a.sort()
s = [0]*(n+1)
for i in range(n):
	s[i+1] = s[i]+a[i]

x = calX(a,n,m,a[0]*2,a[n-1]*2)

out = 0
for i in range(n):
	idx = bisect_left(a,x-a[i])
	out += s[n]-s[idx]+a[i]*(n-idx)

out -= (isOver(a,n,x)-m)*x

print(out)

Solution (Python):
from bisect import bisect_left
n, m = map(int, input().split())
A = list(map(int, input().split()))
A.sort()
A_r = sorted(A, reverse = True)
cumsumA = [0]
total = 0
for i in range(n):
    total += A_r[i]
    cumsumA.append(total)

def func(x):
    ans = 0
    for i in range(n):
        kazu = x - A[i]
        idx = bisect_left(A, kazu)
        cou = n - idx
        ans += cou
    if ans >= m:
        return True
    else:
        return False

mi = 0
ma = 2 * 10 ** 5 + 1
while ma - mi > 1:
    mid = (mi + ma) // 2
    if func(mid):
        mi = mid
    else:
        ma = mid

ans = 0
count = 0

for ai in A_r:
    idx = bisect_left(A, mi - ai)
    ans += ai * (n - idx) + cumsumA[n - idx]
    count += n - idx
print(ans - (count-m) * mi)