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

import sys
input = sys.stdin.readline

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

C = [0]*(10**5+1)
for a in A:
    C[a] += 1

L = 2*10**5+4

CC = rfft(np.array(C), L)
X = (irfft(CC*CC)+0.5).astype(int)
count = 0
ans = 0
for i in reversed(range(L)):
    x = X[i]
    if count + x >= M:
        ans += i*(M-count)
        break
    ans += i*x
    count += x

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 [[]]*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))



from bisect import bisect_left



class Range_BIT():
    def __init__(self, N):
        self.size = N
        self.data0 = [0]*(N+1)
        self.data1 = [0]*(N+1)

    def _add(self, data, k, x):
        while k <= self.size:
            data[k] += x
            k += k & -k

    
    def add(self, l, r, x):
        self._add(self.data0, l, -x*(l-1))
        self._add(self.data0, r, x*(r-1))
        self._add(self.data1, l, x)
        self._add(self.data1, r, -x)

    def _get(self, data, k):
        s = 0
        while k:
            s += data[k]
            k -= k & -k
        return s

    
    def query(self, l, r):
        return self._get(self.data1, r-1)*(r-1)+self._get(self.data0, r-1)\
                -self._get(self.data1, l-1)*(l-1)-self._get(self.data0, l-1)

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

def is_ok(x):
    num = 0
    for a in A:
        p = bisect_left(A,x-a)
        num += N-p
    if num>=M:
        return True
    else:
        return False


l = 0
r = 2*10**5 + 1
while abs(r-l)>1:
    mid = (l+r)//2
    if is_ok(mid):
        l = mid
    else:
        r = mid


bit = Range_BIT(N)
for i in range(1,N+1):
    p = bisect_left(A,l-A[i-1])
    if p!=N:
        bit.add(p+1,N+1,1)
        bit.add(i,i+1,N-p)
x = [0]*N
for i in range(N):
    x[i] = bit.query(i+1,i+2)

ans = 0
for i in range(N):
    ans += A[i]*x[i]


ans -= l*((sum(x)-2*M)//2)

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):
    Y = np.searchsorted(A,x-A)
    return N*N - Y.sum()

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

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

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

Solution (Python):
N, M = map(int, input().split())
a = [2*10**5+1] + sorted(map(int, input().split()), reverse=True)


def solve(x):
   count = 0
   j = N
   for i in range(1, N+1):
      while a[i] + a[j] < x:
         j -= 1
      count += j
   
   return count

left, right = 0, 2*10**5+1
while left < right - 1:
   mid = (left + right) // 2
   if solve(mid) >= M:
      left = mid
   else:
      right = mid

acc = a[:]
acc[0] = 0
for i in range(N):
   acc[i+1] += acc[i]

ans = 0
j = N
for i in range(1, N+1):
   while a[i] + a[j] < left+1:
      j -= 1
   ans += a[i]*j + acc[j]
   M -= j

if M:
   ans += left * M

print(ans)

Solution (Python):
N, M = map(int, input().split())
a = [10**18] + sorted(map(int, input().split()), reverse=True)


def solve(x):
   count = 0
   j = N
   for i in range(1, N+1):
      while a[i] + a[j] < x:
         j -= 1
      count += j
   
   return count


left, right = 0, 10**18
while left < right - 1:
   mid = (left + right) // 2
   if solve(mid) >= M:
      left = mid
   else:
      right = mid

acc = a[:]
acc[0] = 0
for i in range(N):
   acc[i+1] += acc[i]

ans = 0
j = N
for i in range(1, N+1):
   while a[i] + a[j] < left+1:
      j -= 1
   ans += a[i]*j + acc[j]
   M -= j

if M:
   ans += left * M

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 = 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=[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):
import bisect

N, M = map(int, input().split())
A = list(map(int, input().split()))
A.sort()
rev_A = A[::-1]

def check(x, a=A, m=M):
  result = False
  n = len(a)
  count = 0
  for ai in a:
    count += (n - bisect.bisect_left(a, x-ai))
  
  return count>=m


left = 0
right = 2*10**5 + 1
mid = (left + right) // 2
while right - left > 1:
  if check(mid):
    left = mid
  else:
    right = mid
    
  mid = (left + right)//2

x = left
if check(right):
  x = right

accum_A = [rev_A[0]]
for i, a in enumerate(rev_A):
  if i == 0:
    continue
  accum_A.append(accum_A[-1] + a)

ans = 0
count = 0
surplus = 2*10**5
for ai in rev_A:
  idx = bisect.bisect_left(A, x-ai)
  if idx == N:
    break
  ans += (N-idx)*ai + accum_A[N-idx-1]
  count += N - idx
  surplus = min(surplus, ai+A[idx])


print(ans - (count-M)*surplus)

Solution (Python):
from itertools import accumulate
n, m, *a = map(int, open(0).read().split())
a = [-x for x in a]
a.sort()
acc = list(accumulate(a))
l, r = 0, 2 * 10 ** 5 + 1
while r - l > 1:
	mid = (l + r) // 2
	border = n - 1
	cnt = 0
	for x in a:
		while x + a[border] > -mid:
			border -= 1
			if border < 0:
				break
		if border < 0:
			break
		cnt += border + 1
		if cnt >= m:
			break
	if cnt >= m:
		l = mid
	else:
		r = mid
ans = tot = 0
border = n - 1
for x in a:
	while x + a[border] > -l:
		border -= 1
		if border < 0:
			break
	if border < 0:
		break
	ans += acc[border] + x * (border + 1)
	tot += border + 1
ans += l * (tot - m)
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 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

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):
n,m = map(int, input().split())
A = sorted(list(map(int, input().split())),reverse=True)

def cnt(x):
    c = 0
    i = 0
    for j in reversed(range(len(A))):
        while i <n and A[j]+A[i]>=x:
            i += 1
        c += i
    return c

def happiness(x):
    c = 0
    i = 0
    ci = 0
    for j in reversed(range(len(A))):
        while i <n and A[j]+A[i]>=x:
            ci += A[i]
            i += 1
        c += ci + A[j]*i
    return c

def binary_search(l, r):
    mid = (l+r)//2
    c = cnt(mid)
    if c < m:
        return binary_search(l, mid)   
    else:
        if l == mid:
            return(l, c-m)
        return binary_search(mid, r)

x, dm = binary_search(0, A[0]*2+1)
print(happiness(x)-dm*x)

Solution (Python):
from itertools import accumulate
from operator import mul
N, M = map(int, input().split())
A = list(map(int, input().split()))
B = [0] * (10 ** 5 * 2 + 1)
for a in A:
    B[a] += 1
C = list(accumulate(B[::-1]))[::-1]
l, r = 1, 10 ** 5 * 2 + 1
while r - l > 1:
    m = l + (r - l + 1) // 2
    if sum(C[max(0, m - a)] for a in A) >= M:
        l = m
    else:
        r = m
D = list(accumulate(map(mul, range(10 ** 5 * 2, -1, -1), B[::-1])))[::-1]
ans = sum(D[max(0, l - a)] - (l - a) * C[max(0, l - a)] for a in A) + l * M
print(ans)

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

nA, sA, i = [N, N], [0],  N - 1
for a in range(1, A[0] * 2):
    di = 0
    while i - di >= 0 and A[i - di] == a:
        di += 1
    nA.append(nA[-1] - di)
    i -= di
for a in A:
    sA.append(sA[-1] + a)

L, R, ans = 0, A[0] * 2 + 1, 0
while R - L > 1:
    mid = L + (R - L) // 2
    mA = [nA[max(mid - A[i], 0)] for i in range(N)]
    sM = sum(mA)
    if M <= sM:
        L = mid
        ans = sum([A[i] * mA[i] + sA[mA[i]]
                   for i in range(N)]) - mid * (sM - M)
    else:
        R = mid

print(ans)

Solution (Python):
from itertools import accumulate
from bisect import *
n, m, *a = map(int, open(0).read().split())
a = [-x for x in a]
a.sort()
acc = list(accumulate(a)) + [0]
l = 0
r = 2 * 10 ** 5 + 1
while r - l > 1:
	mid = (l + r) // 2
	border = n - 1
	cnt = 0
	for x in a:
		while x + a[border] > -mid:
			border -= 1
			if border < 0:
				break
		if border < 0:
			break
		cnt += border + 1
		if cnt >= m:
			break
	if cnt >= m:
		l = mid
	else:
		r = mid
ans = tot = 0
for x in a:
	b = bisect(a, -l - x)
	tot += b
	ans += acc[b - 1] + x * b
ans += l * (tot - m)
print(-ans)

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


def combinations(x):
  s = 0
  i = 0
  for j in reversed(range(N)):
    while i<N and A[i]+A[j]>=x:
      i += 1
    s += i
  return s
 

def koufukudo(x):
  s = 0
  si = 0
  i = 0
  for j in reversed(range(N)):
    while i<N and A[i]+A[j]>=x:
      si += A[i]
      i += 1
    s += si + A[j]*i  
  return s
 

def bsearch(l, u):
  m = (l+u)//2
  c = combinations(m)

  if c<M:
    return bsearch(l, m)
  else:
    if l==m:
      return (l, c-M)
    return bsearch(m, u)
 
x, dm = bsearch(0, A[0]*2+1)
print(koufukudo(x)-dm*x)

Solution (Python):
import sys
from itertools import accumulate
n, m = [int(i) for i in sys.stdin.readline().split()]
a_ls = [int(i) for i in sys.stdin.readline().split()]
a_ls.sort(reverse=True)
ls = [0 for i in range(10**6)]
for a in a_ls:
    ls[a+1] += 1
ls = [i for i in accumulate(ls)]

def check(x):
    res = 0
    for a in a_ls:
        res += n - ls[max(0,x - a)]
    return res

up = 10**6
down = 0
while True:
    mid = (up + down) // 2
    if check(mid) >= m:
        if down == mid:
            break
        down = mid
    else:
        if up == mid:
            if check(up) == m:
                mid = up
            break
        up = mid

acc_ls = [i for i in accumulate([0] + a_ls)]
ans = 0
for a in a_ls:
    x = n - ls[max(0, mid - a)]
    ans += x * a + acc_ls[x]
print(ans - (check(mid) - m) * mid)

Solution (Python):
from bisect import bisect_left

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

hi = 2*10**5+1
lo = 0



while hi - lo > 1:
    mid = (hi+lo)//2
    
    cnt_hs = 0
    for a in A:
        
        right = mid - a
        idx = bisect_left(A, right)
        
        cnt_hs += N - idx
    
    
    if cnt_hs >= M:
        lo = mid
    else:
        hi = mid
        



B = sorted(A, reverse=True)
accum_B = [0]*N
accum_B[0] = B[0]
for i in range(1,N):
    accum_B[i] = accum_B[i-1] + B[i]

ans = 0
cnt_hs = 0
for a in A:
    
    right = lo - a
    
    if A[-1] < right:
        continue
    idx = bisect_left(A, right)
    
    
    ans += a*(N-idx) + accum_B[N-1 - idx]
    
    cnt_hs += N-idx




if cnt_hs > M:
    ans -= (cnt_hs - M) * lo

print(ans)

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

C=[0]*(10**5*2+1)
for a in A:
    C[a]+=1
for i in range(1,len(C)):
    C[i]+=C[i-1]
 

def pattern_num_ge_x(x):
    p=0
    for a in A:
        p+=(N-(C[x-a-1] if x-a-1>=0 else 0))
    return(p)
 


l=-1
r=10**5*2+1
while r-l>1:
    m=(r+l)//2
    if pattern_num_ge_x(m)<M:
        r=m
    else:
        l=m
 

S=[0]+list(reversed(A[:]))
for i in range(1,N+1):
    S[i]+=S[i-1]
ans=0
for a in A:
    idx=bisect.bisect_left(A,r-a)
    ans+=a*(N-idx) + S[N-idx]
    M-=(N-idx)
 

ans+=l*M
 
print(ans)

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


C=[0]*(10**5*2+1)
for a in A:
    C[a]+=1
for i in range(1,len(C)):
    C[i]+=C[i-1]


def pattern_num_ge_x(x):
    p=0
    for a in A:
        p+=(N-(C[x-a-1] if x-a-1>=0 else 0))
    return(p)



l=-1
r=10**5*2+1
while r-l>1:
    m=(r+l)//2
    if pattern_num_ge_x(m)<M:
        r=m
    else:
        l=m


S=[0]+list(reversed(A[:]))
for i in range(1,N+1):
    S[i]+=S[i-1]
ans=0
for a in A:
    idx=bisect.bisect_left(A,r-a)
    ans+=a*(N-idx) + S[N-idx]
    M-=(N-idx)


ans+=l*M

print(ans)

Solution (Python):
import bisect

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

ans = 0
ok, ng = 0, 1e6
while abs(ng-ok) > 1:
    mid = (ok+ng) // 2
    t = 0
    for a in A:
        t += n - bisect.bisect_left(A, mid-a)
    if t >= m: ok = mid
    else: ng = mid

ans = 0
t = 0
for i in range(n):
    k = bisect.bisect_left(A, ng - A[i])
    t += n - k
    ans += cumA[n] - cumA[k] + (n-k)*A[i]
ans += ok*(m-t)
print(int(ans))

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


def pattern_num_ge_x(x):
    p=0
    for a in A:
        p+=(N-bisect.bisect_left(A, x-a))
    return(p)



l=-1
r=10**5*2+1
while r-l>1:
    m=(r+l)//2
    if pattern_num_ge_x(m)<M:
        r=m
    else:
        l=m


S=[0]+list(reversed(A[:]))
for i in range(1,N+1):
    S[i]+=S[i-1]
ans=0
for a in A:
    idx=bisect.bisect_left(A,r-a)
    ans+=a*(N-idx) + S[N-idx]
    M-=(N-idx)


ans+=l*M

print(ans)

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


def pattern_num_ge_x(x):
    p=0
    for a in A:
        p+=(N-bisect.bisect_left(A, x-a))
    return(p)



l=-1
r=10**10+1
while r-l>1:
    m=(r+l)//2
    if pattern_num_ge_x(m)<M:
        r=m
    else:
        l=m


S=[0]+list(reversed(A[:]))
for i in range(1,N+1):
    S[i]+=S[i-1]
ans=0
for a in A:
    idx=bisect.bisect_left(A,r-a)
    ans+=a*(N-idx) + S[N-idx]
    M-=(N-idx)


ans+=l*M

print(ans)

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


def pattern_num_ge_x(x):
    p=0
    for a in A:
        p+=(N-bisect.bisect_right(A, x-1-a))
    return(p)



l=-1
r=10**10+1
while r-l>1:
    m=(r+l)//2
    if pattern_num_ge_x(m)<M:
        r=m
    else:
        l=m


S=[0]+list(reversed(A[:]))
for i in range(1,N+1):
    S[i]+=S[i-1]
ans=0
for a in A:
    idx=bisect.bisect_right(A,r-a-1)
    ans+=a*(N-idx) + S[N-idx]
    M-=(N-idx)


ans+=l*M

print(ans)

Solution (Python):
import cmath
pi = cmath.pi
exp = cmath.exp

import sys
readline = sys.stdin.readline
write = sys.stdout.write

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

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

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

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

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

l = 10**5+1
L = 2**(2*l-1).bit_length()

power = list(map(int,input().split()))
count = [0.0 for _ in range(l)]

for p in power:
    count[p] += 1

fft_exp_t = make_exp_t(L, -2j*pi)
ifft_exp_t = make_exp_t(L, 2j*pi)

count.extend([0]*(L-l))

C = fft(count, L)

CC =[C[i]*C[i] for i in range(L)]

cc = ifft(CC,L)

shake = 0
ans = 0

for i in reversed(range(len(cc))):
    if shake+int(cc[i].real+0.5)<=M:
        ans += i*int(cc[i].real+0.5)
        shake += int(cc[i].real+0.5)
    else:
        ans += i*(M-shake)
        break
        
print(ans)

Solution (Python):

def count(a, k, n):
    l = []
    i, j, count = 0, n-1, 0
    while True:
        tmp = a[i] + a[j]
        if tmp >= k: 
            if j > 0:
                j -= 1
        elif tmp < k:
            l.append(n-j-1)
            count += n - j - 1
            i += 1
        if i >= n:
            break
        tmp = a[i] + a[j]
        if j == 0 and tmp >= k:
            l.append(n-j)
            count += n - j
            i += 1
        if i >= n:
            break
    return [count, l]

n, m = map(int, input().split())
a = list(map(int, input().split()))
a.sort()
l, r = 10**6, 1


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

ans = 0
if l == r:
    tmp = count(a, l, n)
    for i in range(n):
        ans += a[i] * tmp[1][i] * 2
else:
    p = count(a, l, n)
    q = count(a, r, n)
    check = -1
    for i in range(n):
        if p[1][i] != q[1][i]:
            check = i
            break
    ans += (a[i] + a[n - p[1][check] - 1]) * (m - p[0])
    tmp = count(a, l, n)
    for i in range(n):
        ans += a[i] * tmp[1][i] * 2
print(ans)

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


def pattern_num_ge_x(x):
    p=0
    for a in A:
        ng=-1
        ok=N
        while (ok-ng)>1:
            m=(ok+ng)//2
            if a+A[m]>=x:
                ok=m
            else:
                ng=m
        p+=N-ok
    return(p)



l=-1
r=10**10+1
while r-l>1:
    m=(r+l)//2
    if pattern_num_ge_x(m)<M:
        r=m
    else:
        l=m





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

ans=0
for a in A:
    ng=-1
    ok=N
    while (ok-ng)>1:
        m=(ok+ng)//2
        if a+A[m]>=r:
            ok=m
        else:
            ng=m
    ans+=a*(N-ok) + (S[N-ok-1] if N-ok-1>=0 else 0)
    M-=(N-ok)
    


ans+=l*M

print(ans)

Solution (Python):
N,M = map(int,input().split())
A = sorted(list(map(int,input().split())),reverse=True)
bi = [0]*(A[0]+1)
for a in A:
  bi[a]+=1
  
cum = [0]*(A[0]+1)
for i in range(1,A[0]+1):
  cum[i]=cum[i-1]+bi[i]

cum_max=cum[-1]
left = 0
right = A[0]*2+1
while right-1!=left:
  x = (right+left)//2
  num = 0
  for a in A:
    b = max(x-a,1)
    if b<=A[0]:
      num += cum_max-cum[b-1]
    
  if num>=M:
    left=x
  else:
    right=x

ans = 0
cum_A = [0]*(N+1)
for i in range(1,N+1):
  cum_A[i]=cum_A[i-1]+A[i-1]
for a in A:
  b = max(right-a,1)
  if b>A[0]:
    break
  c = cum_max-cum[b-1]
  M-=c
  ans+=(a*c)+cum_A[c]
  

ans+=M*left
print(ans)

Solution (Python):
N,M = map(int,input().split())
A = sorted(list(map(int,input().split())),reverse=True)
bi = [0]*(A[0]+1)
for a in A:
  bi[a]+=1
  
cum = [0]*(A[0]+1)
for i in range(1,A[0]+1):
  cum[i]=cum[i-1]+bi[i]

cum_max=cum[-1]
left = 0
right = A[0]*2+1
while right-1!=left:
  x = (right+left)//2
  num = 0
  for a in A:
    b = max(x-a,1)
    if b<=A[0]:
      num += cum_max-cum[b-1]
    
  if num>=M:
    left=x
  else:
    right=x

ans = 0
cum_A = [0]*(N+1)
for i in range(1,N+1):
  cum_A[i]=cum_A[i-1]+A[i-1]
for a in A:
  b = max(right-a,1)
  if b>A[0]:
    break
  c = cum_max-cum[b-1]
  M-=c
  ans+=(a*c)+cum_A[c]
  

ans+=M*left
print(ans)

Solution (Python):

import sys
from itertools import accumulate
from bisect import bisect_left

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

    right = 2*A[-1]+1
    left = 0

    while right - left > 1:
        mid = (left+right)//2

        
        
        m = 0
        for i in range(N):
            a = A[i]
            res = mid - a
            
            index = bisect_left(A,res)
            m += N-index
        
        if m >= M:
            left = mid
        else:
            right = mid

    
    B = sorted(A,reverse=True)
    accum_B = list(accumulate(B))
    answer = 0
    handshake_count = 0
    for i in range(N):
        res = left - A[i]
        if res > A[-1]:
            continue
        index = bisect_left(A,res)
        answer += (N-index)*A[i]+accum_B[N-index-1]
        handshake_count += (N-index)

    
    if handshake_count>M:
        answer -= left*(handshake_count-M)
    print(answer)
    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 numpy as np 

I = np.array(sys.stdin.read().split(), dtype=np.int64)
n, m = I[:2]
a = I[2:]
a.sort()
ma = a[-1]
cnt = np.bincount(a, minlength=ma+2)
cnt_not_less_than = np.cumsum(cnt[::-1])[::-1]

right = a[:, None]
def possible(border):
    left = border - right
    left = np.maximum(left, 0)
    left = np.minimum(left, ma+1)
    pair_cnt = np.sum(cnt_not_less_than[left])
    return pair_cnt >= m

def main():
    lo = 1; hi = ma * 2 + 1
    while lo + 1 < hi:
        border = (lo + hi) // 2
        if possible(border):
            lo = border
        else:
            hi = border

    res = lo
    left = res - right
    left = np.maximum(left, 0)
    left = np.minimum(left, ma+1)
    c = cnt_not_less_than[left[left != ma+1]]
    pair_cnt = np.sum(c)
    minimum = np.amin(right[left != ma + 1] + a[-c])

    s = np.cumsum(a[::-1])
    ans = np.sum(s[c-1] + right[left != ma+1] * c) - minimum * (pair_cnt - m)
    return ans

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

Solution (Python):


import sys
import math
from bisect import bisect_right as br
from bisect import bisect_left as bl
sys.setrecursionlimit(2147483647)
from heapq import heappush, heappop,heappushpop
from collections import defaultdict
from itertools import accumulate
from collections import Counter
from collections import deque
from operator import itemgetter
from itertools import permutations
mod = 10**9 + 7
inf = float('inf')
def I(): return int(sys.stdin.readline())
def LI(): return list(map(int,sys.stdin.readline().split()))


def f(x):
    cnt = 0 
    
    for i in a:
        idx = bl(a,x-i)
        cnt += n - idx
        
        if cnt >= m:
            return True
    return False

n,m = LI()
a = LI()
a.sort()
ng,ok = 2*max(a)+1,2

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

ans = 0
tmp = 0
b = [0] + list(accumulate(a))

for i in a:
    idx = bl(a,ok-i)
    
    
    ans += i*(n - idx) + (b[-1] - b[idx])
    
    tmp += n - idx


ans = ans - (tmp - m)*ok
print(ans)

Solution (Python):
import sys
from itertools import accumulate

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

ma = a[-1]
cnt = [0] * (ma + 1)
cur = ma
for i in range(n-1, -1, -1):
    if a[i] == cur:
        cnt[cur] += 1
    else:
        while True:
            cur -= 1
            cnt[cur] = cnt[cur+1]
            if cur == a[i]:
                cnt[cur] += 1
                break
for i in range(cur-1, -1, -1):
    cnt[i] = cnt[i+1]


def possible(border):
    pair_cnt = 0
    for i in range(n):
        aj = border - a[i]
        if aj < 0:
            aj = 0
        elif aj > ma:
            continue
        pair_cnt += cnt[aj]
    return pair_cnt >= m

def main():
    s = [0] + list(accumulate(a[::-1]))
    
    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')
    pair_cnt = 0
    for i in range(n):
        aj = res - a[i]
        if aj < 0:
            aj = 0
        elif aj > ma:
            continue
        c = cnt[aj]
        ans += s[c] + a[i] * c
        pair_cnt += c
        minimum = min(minimum, a[i] + a[-c])

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

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

Solution (Python):


import sys
import math
from bisect import bisect_right as br
from bisect import bisect_left as bl
sys.setrecursionlimit(2147483647)
from heapq import heappush, heappop,heappushpop
from collections import defaultdict
from itertools import accumulate
from collections import Counter
from collections import deque
from operator import itemgetter
from itertools import permutations
mod = 10**9 + 7
inf = float('inf')
def I(): return int(sys.stdin.readline())
def LI(): return list(map(int,sys.stdin.readline().split()))



def f(x):
    cnt = 0
    for i in a:
        idx = bl(a,x-i)
        cnt += n - idx
        if cnt >= m:
            return True
    return False

n,m = LI()
a = LI()
a.sort()
ng,ok = 2*max(a)+1,2

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

ans = 0
tmp = 0
b = [0] + list(accumulate(a))

for i in a:
    idx = bl(a,ok-i)
    
    
    ans += i*(n - idx) + (b[-1] - b[idx])
    
    tmp += n - idx


ans = ans - (tmp - m)*ok
print(ans)

Solution (Python):


import sys
import math
from bisect import bisect_right as br
from bisect import bisect_left as bl
sys.setrecursionlimit(2147483647)
from heapq import heappush, heappop,heappushpop
from collections import defaultdict
from itertools import accumulate
from collections import Counter
from collections import deque
from operator import itemgetter
from itertools import permutations
mod = 10**9 + 7
inf = float('inf')
def I(): return int(sys.stdin.readline())
def LI(): return list(map(int,sys.stdin.readline().split()))


def f(x):
    cnt = 0
    for i in a:
        idx = bl(a,x-i)
        cnt += n - idx
        if cnt >= m:
            return True
    return False

n,m = LI()
a = LI()
a.sort()
ok,ng = 2*max(a)+1,2
while ok - ng > 1:
    mid = (ok+ng) // 2
    if f(mid):
        ng = mid
    else:
        ok = mid

ans = 0
tmp = 0
b = [0] + list(accumulate(a))
for i in a:
    idx = bl(a,ng-i)
    ans += i*(n - idx) + (b[-1] - b[idx])
    tmp += n - idx
ans = ans - (tmp - m)*ng
print(ans)

Solution (Python):
def count(pre, a, x):
  cnt = 0
  val = 0
  for i in reversed(range(n)):
    if 2*a[i] < x:
      break
    lo, hi = 0, n-1
    while lo < hi:
      mid = (lo + hi)//2
      if a[mid] + a[i] < x:
        lo = mid + 1
      else:
        hi = mid
    cnt += 1 + 2*(i - lo)
    val += 2*(i-lo+1)*a[i] + (pre[i] - pre[lo])*2
  return cnt, val

n, m = map(int, input().split())
a = list(sorted(map(int, input().split())))
pre = [0]
for i in range(n):
  pre.append(a[i] + pre[-1])
lo, hi = 0, 2*a[-1]
while lo < hi:
  mid = (lo + hi + 1) // 2
  cnt, _ = count(pre, a, mid)
  if cnt < m:
    hi = mid - 1
  else:
    lo = mid
cnt, val = count(pre, a, lo)
print(val - (cnt - m)*lo)

Solution (Python):
import bisect

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

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

def solve(min_power):
    ans = 0
    for pos in range(n):
        tmp = min_power - a[pos]
        next_pos = bisect.bisect_left(a, tmp, pos) - pos
        times = n - pos - next_pos
        if times == n - pos:
            ans += 2 * (times - 1) + 1 
        else:
            ans += 2 * times
    if ans >= m:
        return True
    else:
        return False

ok = 0
ng = 10** 9
while abs(ok - ng) > 1:
    mid =  (ok + ng) // 2 
    if solve(mid):
        ok = mid
    else:
        ng = mid

cnt = 0
ans = 0

for pos in range(n):
    tmp = ok - a[pos]
    next_pos = bisect.bisect_left(a, tmp, pos) - pos
    times = n - pos - next_pos
    if times == n - pos:
        cnt += 2 * (times - 1) + 1 
        ans += 2 * (ruiseki[n] - ruiseki[n-times+1] + (times-1)*a[pos]) + 2*a[pos]
    else:
        cnt += 2 * times
        ans += 2 * (ruiseki[n] - ruiseki[n-times] + (times)*a[pos])

print(ans - (cnt-m)*ok)