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):
from numpy.fft import rfft, irfft
import numpy as np

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

a = np.zeros(2 ** 18)
for Ai in A:
    a[Ai] += 1

fa = rfft(a)
a2 = irfft(fa ** 2)
a2 = np.round(a2).astype(np.int)

cum_a2 = np.cumsum(a2)

index = np.searchsorted(cum_a2, N ** 2 - M)

a2[index] -= cum_a2[index] - (N ** 2 - M)

ans = sum(A) * 2 * N - (a2[: index + 1] * np.arange(index + 1)).sum()
print(ans)

Solution (Python):
import sys
from bisect import bisect, bisect_left
import heapq

def solve():
    input = sys.stdin.readline
    N, M = map(int, input().split())
    A = [int(a) for a in input().split()]
    A.sort()
    S = [0] * (N + 1)
    for i, a in enumerate(A): S[i + 1] = a + S[i]
 
    l, h = 1, A[N-1] * 2 + 1
    while h - l > 1: 
        m = (h + l) // 2
        count = 0
        for i, a in enumerate(A): count += N - bisect_left(A, m - a)
        if count >= M: 
            l = m
        else: h = m

    total = 0
    count = 0
    Q = []
    heapq.heapify(Q)
    min_out = 10000000
    for i, a in enumerate(A):
        j = bisect_left(A, l - a)  
        total += a * (N - j) + (S[N] - S[j])
        count += N - j
        if 0 <= j < N: min_out = min(min_out, a + A[j])

    if count > M:
        total -= min_out * (count - M)

    print(total)

    return 0

if __name__ == "__main__":
    solve()

Solution (Python):
from bisect import bisect as B
LI=lambda:map(int,input().split())
A=0
n,m=LI()
a=sorted(LI())
R=a[::-1]
def Q(x,M):
  w=0
  for i in range(n):
    y=x-R[i]
    c=B(a,y-1)
    w+=n-c
    if w>=M:return 1
  return 0
N=a[-1]*3
O=-1
while N-O>1:
  f=(O+N)//2
  if Q(f,m):O=f
  else:N=f
P=0
for i in range(n):
  e=n-B(a,O-a[i]-1)
  P+=e
  A+=e*a[i]*2
print(A-O*(P-m))

Solution (Python):
N,M=map(int,input().split())
A=list(map(int,input().split()))
A.sort()
A=[-10**9]+A+[10**9]


import bisect

H=A[-2]*2
L=0

while abs(H-L)>2:
    Mid=(H+L)//2
    
    cnt=0
    for i in range(1,N+1):
        cnt+=N+1-bisect.bisect_left(A,Mid-A[i])
        
    if cnt>M:
        L=Mid
    else:
        H=Mid

X=0
Y=0
Yc=0
for j in range(15):
    Mid2=L-4+j
    cnt=0
    for i in range(1,N+1):
        cnt+=N+1-bisect.bisect_left(A,Mid2-A[i])
        
    if cnt>M:
        X=Mid2
        Xc=cnt
    elif cnt==M:
        X=Mid2
        Y=Mid2
        Yc=cnt
    else:
        Y=Mid2
        Yc=cnt
        break



R=[0]
for i in range(1,len(A)):
    R.append(R[-1]+A[i])

ans=0
for i in range(1,N+1):
    cnt=N+1-bisect.bisect_left(A,Y-A[i])
    
    if cnt>0:
        ans+=cnt*A[i]
        ans+=R[N]-R[N-cnt]
        
    
ans+=X*(M-Yc)
print(ans)

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

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

    def bisect_min(self, reva, func,M):
        ok = 10**5 
        ng = 0 
        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)
    for i in range(n):
        tmp = bisect.bisect_left(reva,tgt-reva[i])
        tmp = n - tmp
        res += tmp
    if M <= res:
        return True
    else:
        return False

def solve2(tgt):
    print(tgt)
    if tgt > 12:
        return True
    else:
        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=[0]
for i in range(N):
    r.append(r[i]+a[i])

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):
from bisect import *
n, m = map(int, input().split())
a = list(map(int, input().split()))
a.sort()

left = 0
right = 200001

while right - left > 1:
    mid = left + (right - left) // 2
    count = 0
    for i in range(n-1, -1, -1):
        count += n - bisect_left(a, mid - a[i])
        if count >= m:
            break


    if m <= count:
        left = mid
    else:
        right = mid

    

ans = 0
a_accum = [0]
for i in range(n):
    a_accum.append(a_accum[i] + a[-1-i])


for i in range(n-1, -1, -1):
    count = n - bisect_left(a, left - a[i])
    ans += a_accum[count] + a[i] * count
    m -= count


ans += left * m
print(ans)

Solution (Python):
import bisect
LI=lambda:map(int,input().split())
A=0
n,m=LI()
a=sorted(LI())
R=a[::-1]
def Q(x,M):
    w=0
    for i in range(n):
        y=x-R[i]
        c=bisect.bisect(a,y-1)
        w+=n-c
        if w>=M:return True
    return False
N=a[-1]*3
O=-1
while N-O>1:
    f=(O+N)//2
    if Q(f,m):
        O=f
    else:
        N=f
P=0
for i in range(n):
    e=n-bisect.bisect(a,O-a[i]-1)
    P+=e
    A+=e*a[i]*2
print(A-O*(P-m))

Solution (Python):
from numpy.fft import rfft, irfft
import numpy as np
from bisect import bisect_left, bisect_right

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

a = np.zeros(2 ** 18)
for Ai in A:
    a[Ai] += 1

fa = rfft(a)
a2 = irfft(fa ** 2)
a2 = np.round(a2).astype(np.int)

cum_a2 = np.cumsum(a2)

index = bisect_left(cum_a2, N ** 2 - M)
eli = cum_a2[index] - (N ** 2 - M)

a2[index] -= eli

ans = sum(A) * 2 * N - (a2[: index + 1] * np.arange(index + 1)).sum()
print(ans)

Solution (Python):
from bisect import bisect as B
LI=lambda:map(int,input().split())
A=0
n,m=LI()
a=sorted(LI())
R=a[::-1]
def Q(x,M):
    w=0
    for i in range(n):
        y=x-R[i]
        c=B(a,y-1)
        w+=n-c
        if w>=M:return True
    return False
N=a[-1]*3
O=-1
while N-O>1:
    f=(O+N)//2
    if Q(f,m):
        O=f
    else:
        N=f
P=0
for i in range(n):
    e=n-B(a,O-a[i]-1)
    P+=e
    A+=e*a[i]*2
print(A-O*(P-m))

Solution (Python):
import bisect as bi

def S(T): 
    global N
    res=0
    for i in range(N):
        res+=k(i,T)
    return res

def k(i,T): 
    global A
    N=len(A)
    res=N-bi.bisect_left(A,T-A[i])
    return res

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


l=A[0]*2;r=A[N-1]*2
while r-l>1:    
    m=(l+r)//2
    
    
    if S(m)>=M: 
        l=m
    else:
        r=m



b=[0] 
for i in range(N):
    b.append(b[-1]+A[N-1-i])


ans=0
cnt=0
for i in range(N):
    cnt+=k(i,l)
    ans+=A[i]*k(i,l)+b[k(i,l)]


if cnt>M:
    ans-=l*(cnt-M)

print(ans)

Solution (Python):
from bisect import bisect as B
LI=lambda:map(int,input().split())
A=0
n,m=LI()
a=sorted(LI())
ra=a[::-1]
def Q(x,M):
    w=0
    for i in range(n):
        y=x-ra[i]
        c=B(a,y-1)
        w+=n-c
        if w>=M:return True
    return False
N=a[-1]*3
O=-1
while abs(O-N)>1:
    f=(O+N)//2
    if Q(f,m):
        O=f
    else:
        N=f
O
P=0
for i in range(n):
    e=n-B(a,O-a[i]-1)
    P+=e
    A+=e*a[i]*2
print(A-O*(P-m))

Solution (Python):
import sys
from itertools import accumulate
from bisect import bisect_left as bi_l, bisect_right as bisect_right

n, m, *a = map(int, sys.stdin.read().split())
a.sort()

def possible(border):
    pair_cnt = 0
    for i in range(n):
        aj = border - a[i]
        j = bi_l(a, aj)
        pair_cnt += n - j
    return pair_cnt >= m

def main():
    s = [0] + list(accumulate(a))
    
    lo = 1; hi = a[-1] * 2 + 1
    while lo + 1 < hi:
        border = (lo + hi) // 2
        if possible(border):
            lo = border
        else:
            hi = border
    res = lo
    ans = 0
    minimum = float('inf')
    cnt = 0
    for i in range(n):
        aj = res - a[i]
        j = bi_l(a, aj)
        ans += s[n] - s[j] + a[i] * (n - j)
        cnt += n - j
        if j < n:
            minimum = min(minimum, a[i] + a[j])

    ans -= (cnt - m) * minimum
    return ans


if __name__ == '__main__':
    ans = main()
    print(ans)

Solution (Python):
import bisect
import sys
N,M = map(int,input().split())
A = list(map(int,input().split()))
A.sort()
L = A[:]
R = A[:]

for i in range(N-1):
    
    L[i+1] += L[i]
    R[-i-2] += R[-i-1]

def func(l,r):
    if r-l == 1:
        return l
    m = (l+r)//2
    count = 0
    for i in A:
        count += N-bisect.bisect_left(A,m-i)
    if count >= M:
        return func(m,r)
    return func(l,m)

f = func(0,10**6)
ans = 0
count = 0
for p in range(N):
    x = N-bisect.bisect_left(A,f-A[p])
    if x == 0:
        continue
    ans += A[p]*x + R[-x]
    count += x
ans -= (count-M)*f
print(ans)

Solution (Python):
N, M = map(int, input().split())
A = sorted([int(a) for a in input().split()])
C = [0] * 101010
S = [0] * 101010

for a in A[::-1]:
    C[a] += 1
    S[a] += a
for i in range(101000)[::-1]:
    C[i] += C[i+1]
    S[i] += S[i+1]

def cnt(n):
    re = 0
    for a in A:
        b = min(max(n-a, 0), 101000)
        re += C[b]
    return re
def calc(n):
    re = 0
    for a in A:
        b = min(max(n-a, 0), 101000)
        re += S[b]  + C[b] * a
    return re

l, r = 0, 200001
while r - l > 1:
    m = (l+r) // 2
    if cnt(m) >= M:
        l = m
    else:
        r = m

print(calc(l) - (cnt(l) - M) * l)

Solution (Python):
from bisect import bisect_left as bl
from itertools import accumulate as ac
n,m=map(int,input().split())
a=sorted(list(map(int,input().split())))
s=[0]+list(ac(a[::-1]))
l,r=2*a[-1]*(n**2)+1,0
while l-r>1:
    t=(l+r)//2
    p=0
    for i in range(n-1,-1,-1):
        if a[i]+a[-1]>=t:
            x=bl(a,t-a[i])
            p+=n-x
        else:
            break
    if p<m:
        l=t
    else:
        r=t
c,d=0,0
w=0
for i in range(n-1,-1,-1):
    if a[i]+a[-1]>=r:
        g=bl(a,r-a[i])
        x=n-g
        c+=s[x]+x*a[i]
        w+=x
    else:
        break
print(c-r*max(0,(w-m)))

Solution (Python):
from bisect import bisect_left as bl
from itertools import accumulate as ac
n,m=map(int,input().split())
a=sorted(list(map(int,input().split())))
s=[0]+list(ac(a[::-1]))
l,r=2*a[-1]*(n**2)+1,0
while l-r>1:
    t=(l+r)//2
    p=0
    for i in range(n-1,-1,-1):
        if a[i]+a[-1]>=t:
            x=bl(a,t-a[i])
            p+=n-x
        else:
            break
    if p<m:
        l=t
    else:
        r=t
c,d=0,0
w=0
for i in range(n-1,-1,-1):
    if a[i]+a[-1]>=r:
        g=bl(a,r-a[i])
        x=n-g
        c+=s[x]+x*a[i]
        w+=x
    else:
        break
print(c-r*max(0,(w-m)))

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

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

A.sort()

Col = [0]*(3*10**5+1)
ind = 0
for i, a in enumerate(A):
    while ind <= a:
        Col[ind] = i
        ind += 1

while ind <= 3*10**5:
    Col[ind] = N
    ind += 1

l = 0
r = 2*10**5+1
while r-l > 1:
    m = (l+r)//2
    count = 0
    for i, a in enumerate(A):
        rem = m - a
        if rem < 0:
            count += N
        else:
            count += N - Col[rem]
    if count >= M:
        l = m
    else:
        r = m

b = 0
B = [0]
for a in A:
    b += a
    B.append(b)



C = l
ans = 0
Count = 0
for i, a in enumerate(A):
    rem = C-a
    if rem < 0:
        Count += N
        ans += B[N] + a*N
    else:
        c = Col[rem]
        pittari = Col[rem+1] - c
        Count += N-c
        ans += B[-1] - B[-(N-c)-1] + a*(N-c)

ans += (M-Count)*C
print(ans)

Solution (Python):
from heapq import heappush,heappop,heapify
from collections import deque,defaultdict,Counter
import itertools
from itertools import permutations,combinations
import sys
import bisect
import string
import math
import time
import random
def I():
    return int(input())
def MI():
    return map(int,input().split())
def LI():
    return [int(i) for i in input().split()]
def LI_():
    return [int(i)-1 for i in input().split()]
def StoI():
    return [ord(i)-97 for i in input()]
def ItoS(nn):
    return chr(nn+97)
def GI(V,E,Directed=False,index=0):
    org_inp=[]
    g=[[] for i in range(n)]
    for i in range(E):
        inp=LI()
        org_inp.append(inp)
        if index==0:
            inp[0]-=1
            inp[1]-=1
        if len(inp)==2:
            a,b=inp
            g[a].append(b)
            if not Directed:
                g[b].append(a)
        elif len(inp)==3:
            a,b,c=inp
            aa=(inp[0],inp[2])
            bb=(inp[1],inp[2])
            g[a].append(bb)
            if not Directed:
                g[b].append(aa)
    return g,org_inp
def bit_combination(k,n=2):
    rt=[]
    for tb in range(n**k):
        s=[tb//(n**bt)%n for bt in range(k)]
        rt+=[s]
    return rt
def show(*inp,end='\n'):
    if show_flg:
        print(*inp,end=end)
YN=['Yes','No']
mo=10**9+7
inf=float('inf')
l_alp=string.ascii_lowercase
u_alp=string.ascii_uppercase
ts=time.time()

input=lambda: sys.stdin.readline().rstrip()

def ran_input():
    import random
    n=random.randint(4,16)
    rmin,rmax=1,10
    a=[random.randint(rmin,rmax) for _ in range(n)]
    return n,a

show_flg=False


ans=0

n,m=LI()
a=sorted(LI())
ra=a[::-1]


def Bigger(x,M): 
    cnt=0
    for i in range(n):
        y=x-ra[i]
        c=bisect.bisect(a,y-1)
        cnt+=n-c

        if cnt>=M:
            return True
    return False


def check(x):
    if x==condition:
        rt=True
    else:
        rt=False
    return rt


ng=1+a[-1]+a[-1]
ok=-1

while abs(ok-ng)>1:
    mid=(ok+ng)//2
    if Bigger(mid,m):
        ok=mid
    else:
        ng=mid

mean=ok





choice=0
for i in range(n):
    num=n-bisect.bisect(a,mean-a[i]-1)
    choice+=num
    ans+=num*a[i]*2



print(ans-mean*(choice-m))

Solution (Python):

from collections import defaultdict,deque
from heapq import heappush, heappop
from itertools import permutations
import sys
import math
import bisect
import random
def LI(): return [int(x) for x in sys.stdin.readline().split()]
def I(): return int(sys.stdin.readline())
def LS():return [list(x) for x in sys.stdin.readline().split()]
def S():
    res = list(sys.stdin.readline())
    if res[-1] == "\n":
        return res[:-1]
    return res
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)]

sys.setrecursionlimit(1000000)
mod = 1000000007

def solve():
    n,m = LI()
    a = LI()
    a.sort()
    l = 0
    r = 2*a[-1]
    while r-l > 1:
        x = (l+r) >> 1
        s = 0
        for i in a:
            j = bisect.bisect_left(a,x-i)
            s += n-j
        if s <= m:
            r = x
        else:
            l = x
    s = [0]*n
    s[-1] = a[-1]
    for i in range(n-1)[::-1]:
        s[i] = s[i+1]+a[i]
    res = m
    ans = 0
    for i in a:
        j = bisect.bisect_left(a,r-i)
        if res < n-j:
            j = n-res
        if j < n:
            ans += min(res,(n-j))*i
            ans += s[j]
            res -= n-j
            if res < 0:
                res = 0
    print(ans+res*l)
    return


if __name__ == "__main__":
    solve()

Solution (Python):
from bisect import bisect_right
N,M = map(int,input().split())
A = list(map(int,input().split()))
B = [-a for a in A]
B.sort()
C = [0]
for b in B:
    C.append(C[-1] + b)

def is_ok(m):
    p = 0
    for b in B:
        i = bisect_right(B, -b-m)
        p += i
        if p >= M: return True
    return False

ok = 0
ng = 2*10**5+1

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

def calc(m):
    p = s = 0
    for b in B:
        i = bisect_right(B, -b-m)
        p += i
        s += C[i] + i*b
    return (p,-s)

p1,s1 = calc(ok)
p2,s2 = calc(ng)
ans = (s1-s2) // (p1-p2) * (M-p2) + s2
print(ans)

Solution (Python):
n, m = [int(i) for i in input().split()]
a = [int(i) for i in input().split()]
a.sort()
b = [0]
for i in reversed(a):
    b.append(b[-1] + i)


def search(i):
    if a[-1] < i:
        return 0
    mi = 0
    ma = n - 1
    while mi != ma:
        m = (mi + ma) // 2
        if a[m] >= i:
            ma = m
        else:
            mi = m + 1
    return n-mi

def countmax(i):
    ans = 0
    for j in a:
        ans += search(i - j)
    return ans

mi = 0
ma = a[-1] * 2
while mi != ma:
    md = (mi + ma) // 2 + 1
    
    if countmax(md) > m:
        mi = md
    else:
        ma = md - 1

ans = 0
for j in a:
    i = search(mi + 1 - j)
    ans+=b[i]+i*j
ans += (m-countmax(mi+1))*mi
print(ans)

Solution (Python):
from bisect import bisect_left

N, M = map(int, input().split())
A = [int(i) for i in input().split()]

A.sort()

def calc(m) :
  c = 0
  for i in range(N) :
    c += N - bisect_left(A, m - A[i])
  return c

ng, ok = 0, A[-1] * 2 + 1
while ok - ng > 1 :
  m = (ok + ng) // 2

  if calc(m) <= M :
    ok = m
  else :
    ng = m
  
S = [A[0]]
for i in range(1, N) :
  S.append(S[-1] + A[i])
S.append(0)
S.append(0)

all = sum(A)
ret = 0
for i in range(N) :
  ret += all - S[bisect_left(A, ng - A[i])-1]
  ret += A[i] * (N - bisect_left(A, ng - A[i]))
ret -= (calc(ng) - M) * (ng)
print(ret)

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

A.sort()

OK=0
NG=max(A)*2+1

import bisect

while NG-OK>1: 
    mid=(OK+NG)//2

    score=0

    for a in A:
        x=bisect.bisect_left(A,mid-a)
        score+=N-x
        if score>=M:
            OK=mid
            break
    else:
        NG=mid
        
from itertools import accumulate
S=list(accumulate(A))

ANS=0
ups=0
for a in A:
    x=bisect.bisect_right(A,OK-a)
    if x!=0:
        ANS+=(N-x)*a+(S[-1]-S[x-1])
    else:
        ANS+=(N-x)*a+S[-1]
        
    
    ups+=N-x

ANS+=(M-ups)*OK
print(ANS)

Solution (Python):
N,M = map(int,input().split())
A = list(map(int,input().split()))
A.sort()
sumA = [0]
for i in range(N):
	sumA.append(sumA[-1] + A[i])

import bisect
def num_shakes(border):
	ret = 0
	for i in range(N):
		if A[i] + A[-1] < border:
			continue
		rest = border - A[i]
		j = bisect.bisect_left(A, rest)
		ret += N - j
	return ret
l,r = A[0]*2, A[-1]*2 + 1 
while r - l > 1:
	p = (l + r)//2
	if num_shakes(p) >= M:
		l = p
	else:
		r = p
m = num_shakes(r)
ans = l * (M-m)
for i in range(N):
	if A[i] + A[-1] < r:
		continue
	rest = r - A[i]
	j = bisect.bisect_left(A, rest)
	ans += A[i] * (N - j) + (sumA[-1] - sumA[j])
print(ans)

Solution (Python):
def f(k):
    cnt=0
    left=N-1
    
    for i in range(N):
        while A[i]+A[left]<k:
            left-=1
            if left<0:
                break
        else:
            cnt+=left+1
            continue
        break
        
    if cnt>=M:
        return True
    else:
        return False

N,M=map(int,input().split())
A=list(map(int,input().split()))
A.sort(reverse=True)
S=[A[0]]
for i in range(1,N):
    S.append(S[-1]+A[i])
low=0
high=2*10**5+100
while high-low>=1:
    mid=(high+low)//2
    if f(mid):
        if low==mid:
            if f(high):
                low=high
            break
        low=mid
    else:
        high=mid
left=N-1
ans=0
cnt=M
for i in range(N):
    while A[i]+A[left]<low:
        left-=1
        if left<0:
            break
    else:
        ans+=A[i]*(left+1)+S[left]
        cnt-=left+1
        continue
    break
print(ans+cnt*low)

Solution (Python):
import sys
input = sys.stdin.readline
n, m = map(int,input().split())
A = list(map(int,input().split()))
A.sort(reverse=True)

import bisect
S = [0] * n
S[0]=A[0]
for i in range(1,n):
    S[i] = S[i-1]+A[i]


l = -1
r = A[0]*2+1
mi = (l+r)//2
while r-l > 1:
    t = n-1
    kai = 0
    for i in range(n):
        while t >= 0 and A[t] + A[i] < mi:
            t -= 1
        
        kai += t+1
    
    if kai < m:
        r = mi
    else:
        l = mi
    mi = (l+r)//2


t = n-1
ans = 0
kai = 0
for i in range(n):
    while t >= 0 and A[t] + A[i] < l:
        t -= 1
    kai += t+1
    
    if t >= 0:
        ans += S[t] + (t+1)*A[i]


if kai > m:
    ans -= (kai-m)*l
print(ans)

Solution (Python):
N,M = list(map(int,input().split()))
A = list(map(int,input().split()))
A.sort(reverse=True)
B = [0]
b=  0
for i in range(N):
    b+=A[i]
    B.append(b)

def bisect1(x,i):
    l = 0
    r = N-1
    if A[i]+A[l]<=x:
        return 0
    if A[i]+A[r]>x:
        return N
    while r-l>1:
        if A[i]+A[(l+r)//2]>x:
            l = (l+r)//2
        else:
            r = (l+r)//2
    return r

def bisect2():
    l=-1
    r = 2*A[0]
    while r-l>1:
        s = 0
        for i in range(N):
            s += bisect1((l+r)//2,i)
        if s>M:
            l = (l+r)//2
        else:
            r = (l+r)//2
    return r

minv = bisect2()
s = 0
ans = 0
for i in range(N):
    j = bisect1(minv,i)
    ans += j*A[i] + B[j]
    s += j

ans += (M-s)*minv

print(ans)

Solution (Python):
from bisect import bisect_left, bisect_right

maxA = 10**5

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

As.sort()

def isOK(v):
    num = 0
    for A in As:
        i = bisect_left(As, v-A)
        num += N-i
    return num >= M

ng, ok = 2*maxA+1, 0
while abs(ok-ng) > 1:
    mid = (ng+ok) // 2
    if isOK(mid):
        ok = mid
    else:
        ng = mid

accRevAs = [As[-1]]
for A in reversed(As[:-1]):
    accRevAs.append(accRevAs[-1]+A)
accRevAs.reverse()

v = ok
ans = 0
numGt = 0

for A in As:

    iR = bisect_right(As, v-A)
    if iR < N:
        ans += accRevAs[iR] + A*(N-iR)
    numGt += N-iR


ans += v * (M-numGt)
print(ans)

Solution (Python):
from bisect import bisect_left


def solve(n, m, a_list):

    a_list_s = list(sorted(a_list))
    a_list_rs_cum = []
    cum = 0
    for v in reversed(a_list_s):
        cum += v
        a_list_rs_cum.append(cum)
    
    left = 0
    right = 2 * 10 ** 5 + 1

    while right > left + 1:
        mid = (left + right) // 2
        res_mid = 0
        for i in range(n):
            res_mid += n - bisect_left(a_list_s, mid - a_list[i])
        if res_mid >= m:
            left = mid
        else:
            right = mid
    
    res = 0
    res_cnt = 0
    for i in range(n):
        c = n - bisect_left(a_list_s, left - a_list[i])
        if c > 0:
            res += a_list_rs_cum[c - 1] + c * a_list[i]
            res_cnt += c
            

    if res_cnt > m:
        res -= left * (res_cnt - m)
    
    return res


def main():
    n, m = map(int, input().split())
    a_list = list(map(int, input().split()))
    res = solve(n, m, a_list)
    print(res)


def test():
    assert solve(5, 3, [10, 14, 19, 34, 33]) == 202
    assert solve(9, 14, [1, 3, 5, 110, 24, 21, 34, 5, 3]) == 1837
    assert solve(9, 73, [67597, 52981, 5828, 66249, 75177, 64141, 40773, 79105, 16076]) == 8128170


if __name__ == "__main__":
    test()
    main()

Solution (Python):
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,copy,functools
import time,random

sys.setrecursionlimit(10**7)
inf = 10**20
eps = 1.0 / 10**10
mod = 10**9+7
dd = [(-1,0),(0,1),(1,0),(0,-1)]
ddn = [(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1)]

def LI(): return list(map(int, sys.stdin.readline().split()))
def LLI(): return [list(map(int, l.split())) for l in sys.stdin.readlines()]
def LI_(): return [int(x)-1 for x in sys.stdin.readline().split()]
def LF(): return [float(x) for x in sys.stdin.readline().split()]
def LS(): return sys.stdin.readline().split()
def I(): return int(sys.stdin.readline())
def F(): return float(sys.stdin.readline())
def S(): return input()
def pf(s): return print(s, flush=True)
def pe(s): return print(str(s), file=sys.stderr)
def JA(a, sep): return sep.join(map(str, a))
def JAA(a, s, t): return s.join(t.join(map(str, b)) for b in a)

def bs(f, mi, ma):
    mm = -1
    while ma > mi:
        mm = (ma+mi) // 2
        if f(mm):
            mi = mm + 1
        else:
            ma = mm
    if f(mm):
        return mm + 1
    return mm

def main():
    n,m = LI()
    a = sorted(LI(), reverse=True)
    t = a[0]
    r = [0] * (t*2+1)
    u = [0] * (t*2+1)
    for c in a:
        r[c] += 1
        u[c] += c

    for i in range(t-1,-1,-1):
        r[i] += r[i+1]
        u[i] += u[i+1]

    def f(i):
        t = 0
        for c in a:
            if c >= i:
                t += n
            else:
                t += r[i-c]
        return t >= m

    i1 = bs(f,0,t*2+1) - 1

    rr = 0
    k = 0
    for c in a:
        ci = max(i1-c,0)
        ac = r[ci]
        if ac == 0:
            continue
        rr += u[ci]
        rr += c * ac
        k += ac

    rr -= (k-m) * i1

    return rr


print(main())

Solution (Python):
import sys
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines




import numpy as np

N,M = map(int,readline().split())
A = np.array(read().split(),np.int64)

A.sort()

def shake_cnt(x):
    
    
    X = np.searchsorted(A,x-A)
    return N * N - X.sum()

left = 0 
right = 10 ** 6 
while left + 1 < right:
    x = (left + right) // 2
    if shake_cnt(x) >= M:
        left = x
    else:
        right = x

left,right



X = np.searchsorted(A,right-A) 
Acum = np.zeros(N+1,np.int64) 
Acum[1:] = np.cumsum(A)

shake = N * N - X.sum()
happy = (Acum[-1] - Acum[X]).sum() + (A * (N - X)).sum()

happy += (M - shake) *left
print(happy)

Solution (Python):
import os
import sys

import numpy as np

if os.getenv("LOCAL"):
    sys.stdin = open("_in.txt", "r")

sys.setrecursionlimit(10 ** 9)
INF = float("inf")
IINF = 10 ** 18
MOD = 10 ** 9 + 7



N, M = list(map(int, sys.stdin.buffer.readline().split()))
A = list(map(int, sys.stdin.buffer.readline().split()))

A.sort()
A = np.array(A, dtype=int)
cs = A[::-1].cumsum()[::-1]


def test(k):
    
    counts = N - np.searchsorted(A, k - A)
    return counts.sum() >= M


def count(k):
    counts = N - np.searchsorted(A, k - A)
    cnt = counts.sum()

    pos = np.searchsorted(A, k - A)
    ret = ((A * counts)[pos < N] + cs[pos[pos < N]]).sum()
    ret -= abs(M - cnt) * k
    return ret


ok = 0
ng = max(A) * 2 + 1
while abs(ok - ng) > 1:
    mid = (ok + ng) // 2
    if test(mid):
        ok = mid
    else:
        ng = mid

ans = count(ok)
print(ans)

Solution (Python):
import numpy as np
N, M = map(int, input().split())
A = list(map(int, input().split()))
B = [0]*(10**5+1)
for ai in A:
    B[ai]+=1
from numpy.fft import rfft, irfft
import numpy as np
B = np.array(B)
fft_len = 2**20
_B = _A = rfft(B, fft_len)
C = list(np.round(irfft(_A*_A, fft_len)).astype('int'))
cnt = 0
ans = 0
for i in range(10**5*2+1)[::-1]:
    if cnt+C[i] > M:
        ans+=i*(M-cnt)
        break
    cnt+=C[i]
    ans+=C[i]*i
print(ans)

Solution (Python):
from bisect import bisect_left

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

def check(mid):
    cnt = 0
    for i in range(n):
        s = mid-a[i]
        cnt += n-bisect_left(a, s)
        if cnt >= m:
            return True
    return False

left = 0
right = 10**11
while right-left>1:
    mid = (right+left)//2
    if check(mid):
        left = mid
    else:
        right = mid
cs = [0]
for i in reversed(range(n)):
    cs.append(cs[-1]+a[i])
ans = 0

cnt = 0
for i in range(n):
    j = n-bisect_left(a, right-a[i])
    ans += j*a[i] + cs[j]
    cnt += j
ans += (m-cnt)*left
print(ans)

Solution (Python):
import numpy as np
from numpy.fft import rfft, irfft

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

fft_len = 1<<20
PA = np.bincount(A)
P = (irfft(rfft(PA, fft_len) * rfft(PA, fft_len)) + .5).astype(int)

ans = 0
for i in range(len(P)-1, -1, -1):
    if P[i] > 0:
        if M <= P[i]:
            ans += i * M
            print(ans)
            exit()
        else:
            M -= P[i]
            ans += i * P[i]