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 numpy as np
N, M = [int(x) for x in input().split()]
A = [int(x) for x in input().split()]
A = np.sort(np.array(A))
def shake_hands(x):
    cnt = 0
    X = np.searchsorted(A, x - A, side="left")
    cnt += N * N - X.sum()
    if cnt >= M:
        return 1
    return 0
left = 0
right = 10 ** 6
while True:
    m = (right + left) // 2
    if shake_hands(m) == 1 and shake_hands(m+1) == 0:
        maximum = m
        break
    elif shake_hands(m) == 0 and shake_hands(m-1) == 1:
        maximum = m - 1
        break
    elif shake_hands(m) == 1:
        left = m + 1
    else:
        right = m - 1
X = np.searchsorted(A,maximum+1-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) * (maximum)
print(happy)

Solution (Python):


import sys, math, itertools, collections, bisect
mans = float('inf')
mod = 10 **9 +7
ans = 0
count = 0
pro = 0

N, M = map(int, input().split())
A = sorted(map(int, input().split()))
B = [0] + A[:]
for i in range(N):
  B [i+1] += B[i]



def solve_binary(mid):
  tmp = 0
  for i, ai in enumerate(A):
    tmp += N - bisect.bisect_left(A, mid-ai)
    
    
  return tmp >= M

def binary_search(N):
  ok = 0
  ng = N
  while abs(ok - ng ) > 1:
    mid = (ok + ng) // 2
    if solve_binary(mid):
      ok = mid
    else:
      ng = mid
  return ok


binresult = binary_search(2*10**5+1)



for i, ai in enumerate(A):
  ans += ai*(N - bisect.bisect_left(A, binresult-ai)) + B[N] - B[bisect.bisect_left(A, binresult-ai)]
  
  count += N - bisect.bisect_left(A, binresult-ai)
ans -= binresult * (count-M)
print(ans)

Solution (Python):
import sys
import cmath
input = sys.stdin.buffer.readline
pi = cmath.pi
exp = cmath.exp

""" https://tjkendev.github.io/procon-library/python/fft/fft.html
"""

def make_exp_t(fn, base):
    exp_t = {0: 1}
    temp = fn
    while temp:
        exp_t[temp] = exp(base / temp)
        temp >>= 1
    return exp_t

def fft_dfs(f, s, fn, st, exp_t):
    if fn==2:
        a = f[s]; b = f[s+st]
        return [a+b, a-b]
    fn2 = fn//2; st2 = st*2
    F0 = fft_dfs(f, s   , fn2, st2, exp_t)
    F1 = fft_dfs(f, s+st, fn2, st2, exp_t)
    w = exp_t[fn]; wk = 1.0
    for k in range(fn2):
        U = F0[k]; V = wk * F1[k]
        F0[k] = U + V
        F1[k] = U - V
        wk *= w
    F0.extend(F1)
    return F0

def fft_init(n: int):
    fn = 2**(n-1).bit_length()
    fft_exp_t = make_exp_t(fn, -2j*pi)
    ifft_exp_t = make_exp_t(fn, 2j*pi)

    def fft(f):
        if fn==1: return f
        F = fft_dfs(f, 0, fn, 1, fft_exp_t)
        return F

    def ifft(F):
        if fn==1: return F
        f = fft_dfs(F, 0, fn, 1, ifft_exp_t)
        f = [x/fn for x in f]
        return f

    return fn, fft, ifft

def pad(f, fn):
    return f + [0,]*max(0, fn-len(f))

def roll(f, shift):
    return f[-shift:] + f[:-shift]

def main():
    N, M = map(int, input().split())
    vA = list(map(int, input().split()))
    nH = 10**5 + 1
    vH = [0,]*nH
    for a in vA:
        vH[a] += 1

    fn, fft, ifft = fft_init(nH*2)
    vH = pad(vH, fn)
    
    

    fA = fft(vH)
    
    fC = [a*a for a in fA]
    vC = ifft(fC)
    rC = [int(x.real+.5) for x in vC]
    

    res = 0
    i = len(rC)
    for c in rC[::-1]:
        i -= 1
        if M>=c:
            res += i*c
            M -= c
        else:
            res += i*M
            break
    print(res)

if __name__ == "__main__":
    main()

Solution (Python):
import sys
import cmath
input = sys.stdin.buffer.readline
pi = cmath.pi
exp = cmath.exp

""" https://tjkendev.github.io/procon-library/python/fft/fft.html
"""

def make_exp_t(fn, base):
    exp_t = {0: 1}
    temp = fn
    while temp:
        exp_t[temp] = exp(base / temp)
        temp >>= 1
    return exp_t

def fft_dfs(f, s, fn, st, exp_t):
    if fn==2:
        a = f[s]; b = f[s+st]
        return [a+b, a-b]
    fn2 = fn//2; st2 = st*2
    F0 = fft_dfs(f, s   , fn2, st2, exp_t)
    F1 = fft_dfs(f, s+st, fn2, st2, exp_t)
    w = exp_t[fn]; wk = 1.0
    for k in range(fn2):
        U = F0[k]; V = wk * F1[k]
        F0[k] = U + V
        F1[k] = U - V
        wk *= w
    F0.extend(F1)
    return F0

def fft_init(n: int):
    fn = 2**(n-1).bit_length()
    fft_exp_t = make_exp_t(fn, -2j*pi)
    ifft_exp_t = make_exp_t(fn, 2j*pi)

    def fft(f):
        if fn==1: return f
        F = fft_dfs(f, 0, fn, 1, fft_exp_t)
        return F

    def ifft(F):
        if fn==1: return F
        f = fft_dfs(F, 0, fn, 1, ifft_exp_t)
        f = [x/fn for x in f]
        return f

    return fn, fft, ifft

def pad(f, fn):
    return f + [0,]*max(0, fn-len(f))

def roll(f, shift):
    return f[-shift:] + f[:-shift]

def main():
    N, M = map(int, input().split())
    vA = list(map(int, input().split()))
    nH = 10**5 + 1
    vH = [0,]*nH
    for a in vA:
        vH[a] += 1

    fn, fft, ifft = fft_init(nH*2)
    vH = pad(vH, fn)
    
    vB = roll(vH, 1)

    fA = fft(vH)
    fB = fft(vB)
    fC = [a*b for a,b in zip(fA,fB)]
    vC = ifft(fC)
    rC = [int(x.real+.5) for x in vC]
    print(*rC[50:70], sep=" ", file=sys.stderr)

    res = 0
    for i,c in zip(reversed(range(len(rC))),rC[::-1]):
        if M>=c:
            res += (i-1)*c
            M -= c
        else:
            res += (i-1)*M
            break
    print(res)

if __name__ == "__main__":
    main()

Solution (Python):
from sys import stdout
import bisect
printn = lambda x: stdout.write(x)
inn = lambda : int(input())
inl   = lambda: list(map(int, input().split()))
inm   = lambda:      map(int, input().split())
DBG = True 
BIG = 999999999
R = 10**9 + 7

def ddprint(x):
  if DBG:
    print(x)

def f(x):
    sm = 0
    for i in range(n-1,-1,-1):
        j = bisect.bisect_left(a,x-a[i])
        sm += n-j
        if sm>=m:
            return True
    return False

n,m = inm()
a = inl()
a.sort()
acc = [0]*(n+1)
for i in range(n-1,-1,-1):
    acc[i] = acc[i+1]+a[i]
mn = 2*min(a)
mx = 2*max(a)+1
while mx-mn>=2:
    mid = (mx+mn)//2
    if f(mid):
        mn = mid
    else:
        mx = mid


sm = cnt = 0
for i in range(n):
    j = bisect.bisect_left(a,mx-a[i])
    sm += acc[j]+(n-j)*a[i]
    cnt += n-j
print(sm+mn*(m-cnt))

Solution (Python):

import sys
from bisect import bisect_left
from itertools import accumulate

INF = float("inf")


def solve(N: int, M: int, A: "List[int]"):
    A.sort()

    
    
    left = 0
    right = 3*(10**5)
    while right - left > 1:
        X = (right + left)//2

        
        tot = 0
        for left_hand in A:
            tot += N-bisect_left(A, X-left_hand)

        if tot < M:
            right = X
        else:
            left = X

    
    
    

    
    
    happy = 0

    acc = [0]+list(accumulate(A))
    tot = 0
    for left_hand in A:
        i = bisect_left(A, right-left_hand)
        tot += N-i
        happy += (acc[N]-acc[i]) + left_hand*(N-i)

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

    return


def main():

    def iterate_tokens():
        for line in sys.stdin:
            for word in line.split():
                yield word
    tokens = iterate_tokens()
    N = int(next(tokens))  
    M = int(next(tokens))  
    A = [int(next(tokens)) for _ in range(N)]  
    solve(N, M, A)


if __name__ == '__main__':
    main()

Solution (Python):
import sys
import math
import heapq
import bisect

sys.setrecursionlimit(10**7)
INTMAX = 9223372036854775807
INTMIN = -9223372036854775808
DVSR = 1000000007
def POW(x, y): return pow(x, y, DVSR)
def INV(x, m=DVSR): return pow(x, m - 2, m)
def DIV(x, y, m=DVSR): return (x * INV(y, m)) % m
def LI(): return [int(x) for x in input().split()]
def LF(): return [float(x) for x in input().split()]
def LS(): return input().split()
def II(): return int(input())
def FLIST(n):
    res = [1]
    for i in range(1, n+1): res.append(res[i-1]*i%DVSR)
    return res
def gcd(x, y):
    if x < y: x, y = y, x
    div = x % y
    while div != 0:
        x, y = y, div
        div = x % y
    return y

N,M=LI()
As=LI()
As.sort()
len_As = len(As)
Acum=[As[0]]
for i in range(1, len_As):
    Acum.append(Acum[i-1]+As[i])



ok = 0
ng = 2*max(As)+1

pat = 0
sm = 0
while ng - ok > 1:
    mid_sum = (ok + ng)//2
    pat = 0
    sm = 0
    for A in As:
        A2 = mid_sum - A
        
        lft = bisect.bisect_left(As, A2)
        pat += len_As - lft
        sm += Acum[len_As-1] - (lft and Acum[lft-1]) + (len_As - lft)*A
    

    if pat < M: ng = mid_sum
    else: ok = mid_sum

pat = 0
sm = 0
for A in As:
    A2 = ok - A
    
    lft = bisect.bisect_left(As, A2)
    pat += len_As - lft
    sm += Acum[len_As-1] - (lft and Acum[lft-1]) + (len_As - lft)*A



print(sm - (pat-M)*ok)

Solution (Python):

N, M = map(int, input().split())
tmp = input().split()

A = [int(tmp[i]) for i in range(N)]
A.sort()
cum_A = [None]*N
cum_A[N-1] = A[N-1]
for i in range(1,N):
    cum_A[N-1-i] = cum_A[N-i] + A[N-1-i]
cum_A.append(0)


def geq(A, cum_A, X, N):  
    sum = 0
    val = 0
    index = 0
    for i in range(N):
        while True:
            if index != N and A[i] + A[N - 1 - index] >= X:
                index += 1
            else:
                sum += index
                val += A[i] * index + cum_A[N-index]
                break
    return sum, val


low = 2
high = 2*A[N-1]
while low <= high:
    X = (low + high)//2
    sum, val = geq(A, cum_A, X, N)
    if sum == M:
        break
    elif sum > M:
        low = X+1
    else:
        high = X-1
if sum < M:
    X -= 1
sum, val = geq(A, cum_A, X, N)
ans = val - X*(sum-M)
print(ans)

Solution (Python):
import bisect
import copy

def check(a, b, bar):
    sm = 0
    cnt = 0
    n = len(a)
    for x in a:
        i = bisect.bisect_left(a, bar - x)
        if i == n:
            continue

        cnt += n - i
        sm += b[i] + x * (n-i)
        
    return cnt, sm

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

b = copy.deepcopy(a)
for i,_ in enumerate(b[:-1]):
    b[n-i-2] += b[n-i-1]
    
left = 0
right = b[0] * 2

while right - left > 1:
    middle = ( right + left ) // 2    
    if check(a, b, middle)[0] < m:
        right = middle
    else:
        left = middle
             
print(check(a, b, left)[1] - left * (check(a, b, left)[0]-m) )

Solution (Python):
import numpy as np
from numpy.fft import rfft, irfft
import sys
input=sys.stdin.readline

N, M = map(int, input().split())
fft_len = 1<<18
MAX = 2*10**5
A = np.zeros(fft_len)
for a in [int(x) for x in input().split()]:
    A[a] += 1
F = rfft(A, fft_len)
f = np.rint(irfft(F*F))
n, ans = 0, 0
for i, c in enumerate(f[:MAX+1][::-1]):
    if M < n+c:
        ans += (M-n)*(MAX-i)
        break
    ans += c*(MAX-i)
    n += c
print(int(ans))

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

input = sys.stdin.readline

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

a.sort(reverse = True)

li = [0] * (2*10**5+1)

for i in a:
    li[i] += 1

acm = list(accumulate(li))
acm_a = list(accumulate(a))

l = 0
r = 2 * 10 ** 5

while l <= r:
    mid = (l + r) // 2
    count = 0
    for i in a:
        if mid - i >= 0:
            count += n - acm[mid-i]
        else:
            count += n
    if count >= m:
        l = mid + 1
    else:
        r = mid - 1

sum_idx = 0
ans = 0
b = a[::-1]
for left in a:
    right_min = l - left
    idx = bisect.bisect_left(b,right_min)
    if n - idx > 0:
        sum_idx += n - idx
        ans += left * (n - idx) + acm_a[n - idx - 1]

ans -= l * (sum_idx - m)

print(ans)

Solution (Python):
import sys
sys.setrecursionlimit(10000000)

from bisect import *





INF = 500000
MOD = 10**9+7

N, M = map(int, input().split())
A = [int(x) for x in input().split()]
A.sort()
l, r = 0, 10**6
while l+1 < r:
    m = (l+r)//2
    cnt = 0
    for a in A:
        cnt += N-bisect_left(A, m-a)
    if cnt < M:
        r = m
    else:
        l = m
f = [0]+A[::-1]
for i in range(N):
    f[i+1] += f[i]
cnt, ans = 0, 0
for a in A:
    i = bisect_right(A, l-a)
    ans += f[N-i]+(N-i)*a
    cnt += N-i
print(ans+(M-cnt)*l)

Solution (Python):
import sys
input=sys.stdin.readline
import numpy as np
from numpy.fft import rfft,irfft

n,m=[int(j) for j in input().split()]
l=np.array([int(j) for j in input().split()])

a=np.bincount(l)

fft_len=1<<18
fft = np.fft.rfft
ifft = np.fft.irfft
Ff = fft(a,fft_len)
x=np.rint(ifft(Ff * Ff,fft_len)).astype(np.int64)
p=x.cumsum()
i=np.searchsorted(p,n*n-m)

q=n*n-m-p[i-1]

ans=(x[:i]*np.arange(i,dtype=np.int64)).sum()+i*q

print(int(l.sum()*2*n-ans))

Solution (Python):
from bisect import bisect_left
from bisect import bisect_right



n,m=map(int, input().split())
A=list(map(int, input().split()))
A.sort()
S=[0]*(n+1)
for i in range(1,n+1):
    S[i]=S[i-1]+A[i-1]

def solve(x):
    case=0
    for i in range(n):
        case+=n-bisect_right(A,x-A[i]-1)
    return case

right=2*A[-1]+1
left=0
         
while abs(right-left)>1:
    mid=(right+left)//2
    if solve(mid)>=m:
        left=mid
    else:
        right=mid
    
Kijun=left
ans=0
for i in range(n):
    case=n-bisect_right(A,Kijun-A[i]-1)
    ans+=case*A[i]+(S[n]-S[n-case])
    
ans-=(solve(Kijun)-m)*Kijun
print(ans)

Solution (Python):
import sys
import numpy as np

stdin = sys.stdin

ri = lambda: int(rs())
rl = lambda: list(map(int, stdin.readline().split())) 
rs = lambda: stdin.readline().rstrip()  

N, M = rl()
A = np.array(rl())
A.sort()

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

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


Acum = np.zeros(N+1, dtype=np.int64)
Acum[1:] = A.cumsum()
Y = np.searchsorted(A, right-A)
happy = (Acum[-1] - Acum[Y]).sum() + (A * (N-Y)).sum()
happy += left * (M-shake_cnt(right))
print(happy)

Solution (Python):
from bisect import bisect_left

N, M, *A = map(int, open(0).read().split())
A.sort()


def calc(t):
    cnt = 0
    for ai in A:
        cnt += N - bisect_left(A, t - ai)
    if M <= cnt:
        return True
    else:
        return False


def solve():
    left, right = 0, 2 * 10 ** 5 + 1
    while 1 < right - left:
        mid = (left + right) // 2
        if calc(mid):
            left = mid
        else:
            right = mid
    
    sum_A = [0] * N
    sum_A[-1] = A[-1]
    for i in range(N - 2, -1, -1):
        sum_A[i] = sum_A[i + 1] + A[i]
    
    ans = 0
    add_cnt = 0
    for ai in A:
        index = bisect_left(A, left - ai)
        if index == N:
            continue
        add_cnt += N - index
        ans += (N - index) * ai + sum_A[index]
    ans -= left * (add_cnt - M)
    print(ans)


if __name__ == '__main__':
    solve()

Solution (Python):

n,m=map(int,input().split())
a=sorted(list(map(int,input().split())),reverse=True)
cs=[]
sum=0
for i in range(n):
	sum+=a[i]
	cs.append(sum)


def overX(x):
	count=0
	
	if a[0]*2<x:
		return 0
	else:
		for i in range(n):
			if a[i]+a[0]<x:
				count+=0
			else:
				lb, ub = 0, n
				
				while ub-lb>1:
					mid=(lb+ub)//2
					if a[i]+a[mid]>=x:
						lb=mid
					else:
						ub=mid
				count+=lb+1
		return count



lb=1 
ub=2*10**5+1
while ub-lb>1:
	mid=(ub+lb)//2
	if overX(mid)>=m:
		lb=mid
	else:
		ub=mid

maxx=lb

hp=0


if a[0]*2>=maxx+1:
	for i in range(n):
		if a[i]+a[0]<maxx+1:
			break
		else:
			lb, ub = 0, n
			
			while ub-lb>1:
				mid=(lb+ub)//2
				if a[i]+a[mid]>=maxx+1:
					lb=mid
				else:
					ub=mid
			
			hp+=a[i]*(lb+1)+cs[lb]
hp+=(m-overX(maxx+1))*maxx
print(hp)

Solution (Python):
import bisect

N, M = map(int, input().split())
A = list(map(int, input().split()))
A.sort()
B = [0] * N
B[-1] = A[-1]
for i in range(N - 2, -1, -1):
    B[i] = B[i+1]+A[i]


def C(mid):
    tmp = 0
    for a in A:
        pos = bisect.bisect_right(A, mid - a)
        tmp += N-pos
    return tmp > M


lb = 0
rb = 10**6
while rb - lb > 1:
    happy = (lb + rb) // 2
    if C(happy):
        lb = happy
    else:
        rb = happy

ans = 0
cnt = 0
for a in A:
    pos = bisect.bisect_right(A, rb - a)
    if pos == N:
        continue
    ans += B[pos]+(N-pos)*a
    cnt += N - pos
ans += (M-cnt)*rb
print(ans)

Solution (Python):
import numpy as np

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

n0 = 2**int(np.ceil(np.log2(2*Amax+1)))

Afre=np.zeros(n0).astype(int)
for i in range(N):
    Afre[A[i]]+=1


S = np.rint(np.fft.irfft(np.fft.rfft(Afre)*np.fft.rfft(Afre))).astype(int)
Scum =S.cumsum()
bd = N*N-M
i=np.searchsorted(Scum,bd)

remove=((Scum[-1]-Scum[i])-M)*i
ans=0
for j in range(i+1,n0):
    ans+=S[j]*j
print(ans-remove)

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()))

counts = np.bincount(A)
size = 1 << len(counts).bit_length() + 1
counts2 = np.rint(np.fft.irfft(np.fft.rfft(counts, size) ** 2))

cnt = 0
ans = 0
for i in np.where(counts2 > 0)[0][::-1]:
    cnt += counts2[i]
    ans += i * counts2[i]
    if cnt > M:
        ans -= i * (cnt - M)
        break
print(int(ans))

Solution (Python):


import sys

def input(): return sys.stdin.readline().strip()
def list2d(a, b, c): return [[c] * b for i in range(a)]
def list3d(a, b, c, d): return [[[d] * c for j in range(b)] for i in range(a)]
def list4d(a, b, c, d, e): return [[[[e] * d for j in range(c)] for j in range(b)] for i in range(a)]
def ceil(x, y=1): return int(-(-x // y))
def INT(): return int(input())
def MAP(): return map(int, input().split())
def LIST(N=None): return list(MAP()) if N is None else [INT() for i in range(N)]
def Yes(): print('Yes')
def No(): print('No')
def YES(): print('YES')
def NO(): print('NO')
sys.setrecursionlimit(10 ** 9)
INF = 10 ** 18
MOD = 10 ** 9 + 7

def fft(A, B, L):
    import numpy as np
    from numpy.fft import rfft, irfft

    
    res = irfft(rfft(A, L) * rfft(B, L), L)
    
    res = np.rint(res).astype(np.int64)
    return res

MAX = 10 ** 5
N, M = MAP()
A1 = [0] * (MAX+1)
A2 = [0] * (MAX+1)
for a in LIST():
    A1[a] += 1
    A2[a] += 1

L = 2 ** 18
res = fft(A1, A2, L)

ans = cnt = 0
for a in range(L-1, -1, -1):
    c = res[a]
    if cnt + c <= M:
        cnt += c
        ans += a * c
    else:
        
        ans += a * (M - cnt)
        break
print(ans)

Solution (Python):
from itertools import accumulate

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

a.sort(reverse=True)
a_acm = list(accumulate(a))

a_cnt = [0] * (10 ** 6 + 1)
for e in a:
    a_cnt[e] += 1

acm = list(accumulate(a_cnt))


def count_handshake(x):
    ret = 0
    for e in a:
        if x - e - 1 >= 0:
            cnt = acm[-1] - acm[x-e-1]
        else:
            cnt = n

        ret += cnt

    return ret


left = 0
right = 10 ** 6

while right - left > 1:
    mid = (left + right) // 2
    if count_handshake(mid) >= m:
        left = mid
    else:
        right = mid

ans = 0
handshake_cnt = 0
for e in a:
    if left - e - 1 >= 0:
        cnt = n - acm[left-e-1]
    else:
        cnt = n

    if cnt > 0:
        ans += e * cnt + a_acm[cnt-1]
        handshake_cnt += cnt

if handshake_cnt > m:
    ans -= left * (handshake_cnt - m)

print(ans)

Solution (Python):

import bisect
import heapq
import math
import random
import sys
from collections import Counter, defaultdict, deque
from decimal import ROUND_CEILING, ROUND_HALF_UP, Decimal
from functools import lru_cache, reduce
from itertools import combinations, combinations_with_replacement, product, permutations
from operator import add, mul, sub

sys.setrecursionlimit(100000)
input = sys.stdin.readline
INF = 2**62-1

def read_int():
    return int(input())


def read_int_n():
    return list(map(int, input().split()))


def read_float():
    return float(input())


def read_float_n():
    return list(map(float, input().split()))


def read_str():
    return input().strip()


def read_str_n():
    return list(map(str, input().split()))


def error_print(*args):
    print(*args, file=sys.stderr)


def mt(f):
    import time

    def wrap(*args, **kwargs):
        s = time.time()
        ret = f(*args, **kwargs)
        e = time.time()

        error_print(e - s, 'sec')
        return ret

    return wrap


def f(n):
    r = 2
    ret = 1
    for i in range(n, n-r, -1):
        ret *= i
    for i in range(1, r+1):
        ret //= i
    return ret * 2 + n


class Bisect:
    def __init__(self, func):
        self.__func = func

    def bisect_left(self, x, lo, hi):
        while lo < hi:
            mid = (lo+hi)//2
            if self.__func(mid) < x:
                lo = mid+1
            else:
                hi = mid
        return lo

    def bisect_right(self, x, lo, hi):
        while lo < hi:
            mid = (lo+hi)//2
            if x < self.__func(mid):
                hi = mid
            else:
                lo = mid+1
        return lo

@mt
def slv(N, M, A):
    A.sort(reverse=True)
    print(A)
    b = Bisect(f)

    i = b.bisect_left(M, 1, N)
    if f(i) != M:
        i -= 1
    l = f(i)

    ans = 0
    for j in range(i):
        ans += A[j] * 2
        ans += A[j] * (i-1)

    rem = M - l
    print(i, l, rem)
    j = 0
    while rem != 0:
        ans += A[j] + A[i]
        rem -= 1
        if rem == 0:
            break
        ans += A[j] + A[i]
        rem -= 1
        if rem == 0:
            break
        j += 1


    return ans

import numpy as np
@mt
def slv2(N, M, A):
    C = Counter(A)
    L = 1 << (max(A).bit_length() + 1)
    B = [C[i] for i in range(L)]
    D = np.fft.rfft(B)
    D *= D
    E = list(map(int, np.round(np.fft.irfft(D))))
    ans = 0
    c = 0
    i = L - 1

    while c != M:
        n = E[i]
        if c + n > M:
            n = M-c
        c += n
        ans += i * n
        i -= 1

    return ans

def main():
    N, M = read_int_n()
    A = read_int_n()
    print(slv2(N, M,A))

    
    
    
    


if __name__ == '__main__':
    main()

Solution (Python):

import bisect
import heapq
import math
import random
import sys
from collections import Counter, defaultdict, deque
from decimal import ROUND_CEILING, ROUND_HALF_UP, Decimal
from functools import lru_cache, reduce
from itertools import combinations, combinations_with_replacement, product, permutations
from operator import add, mul, sub

sys.setrecursionlimit(100000)
input = sys.stdin.readline
INF = 2**62-1

def read_int():
    return int(input())


def read_int_n():
    return list(map(int, input().split()))


def read_float():
    return float(input())


def read_float_n():
    return list(map(float, input().split()))


def read_str():
    return input().strip()


def read_str_n():
    return list(map(str, input().split()))


def error_print(*args):
    print(*args, file=sys.stderr)


def mt(f):
    import time

    def wrap(*args, **kwargs):
        s = time.time()
        ret = f(*args, **kwargs)
        e = time.time()

        error_print(e - s, 'sec')
        return ret

    return wrap


def f(n):
    r = 2
    ret = 1
    for i in range(n, n-r, -1):
        ret *= i
    for i in range(1, r+1):
        ret //= i
    return ret * 2 + n


class Bisect:
    def __init__(self, func):
        self.__func = func

    def bisect_left(self, x, lo, hi):
        while lo < hi:
            mid = (lo+hi)//2
            if self.__func(mid) < x:
                lo = mid+1
            else:
                hi = mid
        return lo

    def bisect_right(self, x, lo, hi):
        while lo < hi:
            mid = (lo+hi)//2
            if x < self.__func(mid):
                hi = mid
            else:
                lo = mid+1
        return lo

@mt
def slv(N, M, A):
    A.sort(reverse=True)
    print(A)
    b = Bisect(f)

    i = b.bisect_left(M, 1, N)
    if f(i) != M:
        i -= 1
    l = f(i)

    ans = 0
    for j in range(i):
        ans += A[j] * 2
        ans += A[j] * (i-1)

    rem = M - l
    print(i, l, rem)
    j = 0
    while rem != 0:
        ans += A[j] + A[i]
        rem -= 1
        if rem == 0:
            break
        ans += A[j] + A[i]
        rem -= 1
        if rem == 0:
            break
        j += 1


    return ans

import numpy as np
@mt
def slv2(N, M, A):
    C = Counter(A)
    L = 1 << (max(A).bit_length() + 1)
    B = [C[i] for i in range(L)]
    D = np.fft.rfft(B)
    D *= D
    E = list(map(int, np.round(np.fft.irfft(D).real)))
    ans = 0
    c = 0
    i = L - 1

    while c != M:
        n = E[i]
        if c + n > M:
            n = M-c
        c += n
        ans += i * n
        i -= 1

    return ans

def main():
    N, M = read_int_n()
    A = read_int_n()
    print(slv2(N, M,A))

    
    
    
    


if __name__ == '__main__':
    main()

Solution (Python):

import bisect
import heapq
import math
import random
import sys
from collections import Counter, defaultdict, deque
from decimal import ROUND_CEILING, ROUND_HALF_UP, Decimal
from functools import lru_cache, reduce
from itertools import combinations, combinations_with_replacement, product, permutations
from operator import add, mul, sub

sys.setrecursionlimit(100000)
input = sys.stdin.readline
INF = 2**62-1

def read_int():
    return int(input())


def read_int_n():
    return list(map(int, input().split()))


def read_float():
    return float(input())


def read_float_n():
    return list(map(float, input().split()))


def read_str():
    return input().strip()


def read_str_n():
    return list(map(str, input().split()))


def error_print(*args):
    print(*args, file=sys.stderr)


def mt(f):
    import time

    def wrap(*args, **kwargs):
        s = time.time()
        ret = f(*args, **kwargs)
        e = time.time()

        error_print(e - s, 'sec')
        return ret

    return wrap


def f(n):
    r = 2
    ret = 1
    for i in range(n, n-r, -1):
        ret *= i
    for i in range(1, r+1):
        ret //= i
    return ret * 2 + n


class Bisect:
    def __init__(self, func):
        self.__func = func

    def bisect_left(self, x, lo, hi):
        while lo < hi:
            mid = (lo+hi)//2
            if self.__func(mid) < x:
                lo = mid+1
            else:
                hi = mid
        return lo

    def bisect_right(self, x, lo, hi):
        while lo < hi:
            mid = (lo+hi)//2
            if x < self.__func(mid):
                hi = mid
            else:
                lo = mid+1
        return lo

@mt
def slv(N, M, A):
    A.sort(reverse=True)
    print(A)
    b = Bisect(f)

    i = b.bisect_left(M, 1, N)
    if f(i) != M:
        i -= 1
    l = f(i)

    ans = 0
    for j in range(i):
        ans += A[j] * 2
        ans += A[j] * (i-1)

    rem = M - l
    print(i, l, rem)
    j = 0
    while rem != 0:
        ans += A[j] + A[i]
        rem -= 1
        if rem == 0:
            break
        ans += A[j] + A[i]
        rem -= 1
        if rem == 0:
            break
        j += 1


    return ans

import numpy as np
@mt
def slv2(N, M, A):
    C = Counter(A)
    L = 1 << (max(A).bit_length() + 1)
    B = [C[i] for i in range(L)]
    D = np.fft.rfft(B)
    D *= D
    E = list(map(int, np.round(np.fft.irfft(D).real)))
    ans = 0
    c = 0
    i = L - 1
    while c != M:
        n = E[i]
        if c + n > M:
            n = M-c
        c += n
        ans += i * n
        i -= 1

    return ans

def main():
    N, M = read_int_n()
    A = read_int_n()
    print(slv2(N, M,A))

    
    
    
    


if __name__ == '__main__':
    main()