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
from bisect import *

def judge(x):
    cnt = 0
    
    for i in range(N):
        idx = bisect_left(A, x-A[i])
        cnt += N-idx
    
    return cnt>=M
    
def binary_search():
    l, r = 0, 10**11
    
    while l<=r:
        mid = (l+r)//2
        
        if judge(mid):
            l = mid+1
        else:
            r = mid-1
    
    return r

N, M = map(int, input().split())
A = list(map(int, input().split()))
A.sort()
border = binary_search()
acc = [0]

for Ai in A:
    acc.append(acc[-1]+Ai)

ans = 0
cnt = 0

for i in range(N):
    idx = bisect_left(A, border-A[i]+1)
    ans += A[i]*(N-idx)+acc[N]-acc[idx]
    cnt += N-idx

ans += (M-cnt)*border

print(ans)

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

def cumsum(s):
    n = len(s)
    cs = [0] * (n+1)
    for i in range(n):
        cs[i+1] = cs[i] + s[i]
    return cs

def bs_list(a, f):
    l, r = -1, len(a)
    while r - l > 1:
        x = (l + r) // 2
        if f(a[x]): r = x
        else: l = x
    return None if r == len(a) else r

a.sort()
ca = cumsum(a)

def detect(x):
    num = 0
    for b in a[::-1]:
        res = bs_list(a, lambda y: y >= x - b)
        if res is None: break
        num += n - res
    return num <= m

l, r = -1, 10**5*2+10
while r - l > 1:
    x = (l+r) // 2
    if detect(x): r = x
    else: l = x

s, c = 0, 0
for b in a[::-1]:
    res = bs_list(a, lambda x: x >= r - b)
    if res is None: break
    c += (n - res)
    s += b * (n - res) + (ca[n] - ca[res])

print(s + (m - c) * l)

Solution (Python):
N,M = map(int,input().split())
A = list(map(int,input().split()))
A.sort()
B = [0 for i in range(10**5+1)]
S = [0]
s = 0
for i in range(N):
    B[A[i]] += 1
    s += A[i]
    S.append(s)
for i in range(2,10**5+2):
    B[-i] += B[-i+1]
d = 0
u = 2*10**5
for i in range(30):
    x = (u+d)//2
    if i == 28:
        x += 1
    m = 0
    for i in range(N):
        a = A[i]
        y = max(x - a,0)
        if y <= 10**5:
            m += B[y]
    if m >= M:
        d = x
    else:
        u = x
ans = 0
for i in range(N):
    a = A[i]
    x = max(d - a,0)
    if x <= 10**5:
        y = B[x]
    else:
        y = 0
    ans += (S[N] - S[N-y]) + a*y
ans -= (m-M)*d
print(ans)

Solution (Python):

import sys
from itertools import accumulate
from bisect import bisect_left

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

A2 = sorted(A, reverse=True)
A = sorted(A)




def f1(m):
    cnt = 0
    for a1 in A2:
        tmp = bisect_left(A, m - a1)
        if tmp == N: break
        cnt += N - tmp

    return cnt >= M


l = 0
r = A2[0] * 2 + 10

while r - l > 1:
    m = (l + r) // 2
    if f1(m):
        l = m
    else:
        r = m

X = l


acc = list(accumulate(A2))  

total_cnt = 0
ans = 0
for a in A2:
    
    min_right = X - a
    if min_right > A[-1]:
        continue

    
    idx = bisect_left(A, min_right)
    cnt = N - idx

    
    ans += a * cnt + acc[cnt - 1]
    total_cnt += cnt

    
    

if total_cnt > M:
    ans -= X * (total_cnt - M)
print(ans)

Solution (Python):
import numpy as np
from numpy.fft import fft
from numpy.fft import ifft
class FFT:
  def exe(s, A, B):
    N, M = 1, len(A) + len(B) - 1
    while N < M: N <<= 1
    A, B = s.arrN(N, A), s.arrN(N, B)
    A = fft(A) * fft(B)
    A = ifft(A).real[:M] + 0.5
    return list(map(int, A))
  def arrN(s, N, L):
    return np.zeros(N) + (L + [0] * (N - len(L)))

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

MAX = max(A)
L = [0] * (MAX + 1)
for i in A:
  L[i] += 1

f = FFT()
L = f.exe(L, L)

ans = 0
for i in range(len(L) - 1, -1, -1):
  if L[i] < M:
    ans += i * L[i]
    M -= L[i]
  else:
    ans += i * M
    break

print(ans)

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

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

a.sort()

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

ok = 0
ng = 200001

while ng-ok > 1:
    mid = (ok+ng)//2
    cnt = 0
    for e in a:
        k = bisect.bisect_left(a,mid-e)
        cnt += n-k
    if cnt >= m:
        ok = mid
    else:
        ng = mid

c = 0
res = 0

for e in a:
    k = bisect.bisect_left(a,ok+1-e)
    res += b[k] + e*(n-k)
    c += n-k

res += ok*(m-c)
print(res)

Solution (Python):











import sys
sys.setrecursionlimit(10**6)
input=sys.stdin.readline
from math import floor,ceil,sqrt,factorial,hypot,log 
from heapq import heappop, heappush, heappushpop
from collections import Counter,defaultdict,deque
from itertools import accumulate,permutations,combinations,product,combinations_with_replacement
from bisect import bisect_left,bisect_right
from copy import deepcopy
inf=float('inf')
mod = 10**9+7
def pprint(*A): 
    for a in A:     print(*a,sep='\n')
def INT_(n): return int(n)-1
def MI(): return map(int,input().split())
def MF(): return map(float, input().split())
def MI_(): return map(INT_,input().split())
def LI(): return list(MI())
def LI_(): return [int(x) - 1 for x in input().split()]
def LF(): return list(MF())
def LIN(n:int): return [I() for _ in range(n)]
def LLIN(n: int): return [LI() for _ in range(n)]
def LLIN_(n: int): return [LI_() for _ in range(n)]
def LLI(): return [list(map(int, l.split() )) for l in input()]
def I(): return int(input())
def F(): return float(input())
def ST(): return input().replace('\n', '')
def main():
    N,M=MI()
    A=sorted(LI())  
    
    *cumsum_rev, = accumulate(A[::-1]) 

    def is_ok(x): 
        cnt = 0
        for a in A:
            cnt += N-bisect_left(A, x-a)
        return cnt >= M 

    ok = 0
    ng = 10**7
    while ng-ok>1:
        mid = (ok+ng)//2
        if is_ok(mid):
            ok = mid
        else:
            ng = mid
    
    ans = 0
    cnt_sum = 0
    for a in A:
        cnt = N-bisect_left(A,ok-a)
        cnt_sum += cnt
        if cnt:
            ans += (a * cnt) + cumsum_rev[cnt-1]
        
    ans -= (cnt_sum-M)*ok
    print(ans)
if __name__ == '__main__':
    main()

Solution (Python):


from itertools import accumulate

def my_bisect(a, x, key, lo=0, hi=None, side="left"):
    if hi is None:
        hi = len(a)

    lov = key(a[lo])
    hiv = key(a[hi])

    if side == "left":
        while lo < hi:
            mid = (lo + hi) // 2
            miv = key(a[mid])

            if miv < val: lo, lov = mid+1, miv
            else: hi, hiv = mid, miv
    else:
        while lo < hi:
            mid = (lo + hi) // 2
            miv = key(a[mid])

            if val < miv: hi, hiv = mid, miv
            else: lo, lov = mid+1, miv

    return lo

from bisect import bisect_right

def resolve():
    N, M = map(int, input().rstrip().split())
    a = list(map(int, input().rstrip().split()))

    a.sort(reverse=True)
    b = list(accumulate(a))
    c = [-i for i in a]

    def count_shake(z):
        n = 0
        for i, x in enumerate(a):
            y = x - z
            
            j = bisect_right(c, y, i)
            if i == j:
                break
            n += 2 * (j - i) - 1
        return n



    lo = a[-1] * 2
    hi = a[0] * 2

    while True:
        mid = (lo + hi) // 2
        x = count_shake(mid)

        if hi - lo <= 1:
            break
        if x < M: hi = mid
        else: lo = mid

    s = 0
    d = [-1] * N
    m = x
    for i, x in enumerate(a):
        y = x - mid
        j = bisect_right(c, x - mid)
        if j == 0: break
        d[i] = j-1
        s += b[j-1] + x * j

    print(s - mid * (m-M))




if __name__ == "__main__":
    resolve()

Solution (Python):

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


import numpy as np
A = np.array(A)
A.sort()

def number_of_shake(K):
    p = np.searchsorted(A, K-A)
    return N*N - p.sum()

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

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

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

print(happy)

Solution (Python):
import bisect 
from itertools import accumulate

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

xOK = A[0]*2
xNG = A[-1]*2 + 1
nx = 0
while xOK+1 < xNG:
    x = (xOK + xNG)//2
    nx = 0
    for a in A:
        nx += N-(bisect.bisect_left(A,x-a))
        if nx >= M:
            xOK = x
            break
    else:
        xNG = x

S=[0]+list(accumulate(A))

ans = 0
nnx = 0        
for a in A:
    nx = bisect.bisect_right(A,xOK-a)
    ans += a*(N-nx) + (S[N]-S[nx])
    nnx += (N-nx)
ans += xOK*(M-nnx)
        
    
print(ans)

Solution (Python):
def nibutan(n, M):
  ok = 0
  ng = M + 1
  while abs(ok - ng) > 1:
    mid = (ok + ng) // 2
    
    if nasu(mid) >= n:
      ok = mid
    else:
      ng = mid
  return ok

def nasu(x):
  ans = 0
  for i in A:
    
    t = N - hogetan(x - i)
    
    ans += t
  return ans

def hogetan(n):
  ok = N
  ng = -1
  while abs(ok - ng) > 1:
    mid = (ok + ng) // 2
    if A[mid] >= n:
      ok = mid
    else:
      ng = mid
  return ok

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


A.sort()
MAX = A[-1] + A[-1]
inf = MAX + 10

t = nibutan(M, MAX)
mi = inf
cnt = 0
ans = 0
for i in A:
  k = hogetan(t - i)
  n = N - k
  cnt += n
  if n != 0:
    mi = min(mi, i + A[k])
  ans += n * 2 * i

ans -= mi * (cnt - M)

print(ans)

Solution (Python):
from bisect import bisect_left

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

def doesOver(k):
    cnt = 0
    for a in A:
        cnt += N - bisect_left(A, k - a)
    return cnt >= M

overEq = 0
less = A[-1] * 2 + 100

while less - overEq > 1:
    mid = (less + overEq) // 2
    if doesOver(mid):
        overEq = mid
    else:
        less = mid

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

ans = 0
cnt = 0
for a in A:
    left = bisect_left(A, overEq - a)
    ans += accA[N] - accA[left]
    ans += (N - left) * a
    cnt += N - left

ans -= (cnt - M) * overEq
print(ans)

Solution (Python):
import sys
import bisect
sys.setrecursionlimit(10**7)

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

a.sort()

def hantei(t):
    global n,m,a
    ret = 0
    for sa in a:
        ret += n - bisect.bisect_left(a,t-sa)
    if ret<=m:return True
    else:return False

def nibun(l,r):
    if r == l:
        return l
    f = (l+r)//2
    if hantei(f):
        return nibun(l,f)
    else:
        return nibun(f+1,r)


p=nibun(0,3*(10**5))

a.reverse()
na=[]
sna = 0
for sa in a:
    sna += sa
    na.append(sna)
na.reverse()
na.append(0)
a.reverse()

res = 0
f=0
for sa in a:
    q = bisect.bisect_left(a,p-sa)
    f += n-q
    res += na[q]+sa*(n-q)
res += (p-1)*(m-f)

print(res)

Solution (Python):
from bisect import bisect_left
import numpy


def main():
    def isOK(score):
        count = 0
        for a in reversed(A):
            border_score = score - a
            index = bleft(A, border_score)
            count += N - index
            if count >= M:
                return True
        return False

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

    ng = A[-1] + A[-1] + 1
    ok = A[0] + A[0]
    bleft = bisect_left

    while(abs(ok-ng) > 1):
        
        mid = (ok+ng) // 2
        f = isOK(mid)
        if f:
            ok = mid
        else:
            ng = mid
        
        
    C = [0,] + A
    C = numpy.cumsum(C[::-1])
    
    min_score = float('inf')
    count = 0
    for a in reversed(A):
        index = bleft(A, ok-a)
        if index < N:
            min_score = min(a+A[index], min_score)
            ans += C[N-index-1] + (N-index) * a
            count += N-index

    ans -= min_score * (count-M)
    print(ans)


if __name__ == '__main__':
    main()

Solution (Python):
from numpy import*
n,m,*a=int_(open(0).read().split())
a=int_(fft.irfft(fft.rfft(bincount(a,[1]*n,2**18))**2)+.5)
c=0
for i in where(a>0)[0][::-1]:t=min(m,a[i]);c+=i*t;m-=t
print(c)

Solution (Python):
from numpy import*
n,m,*a=int_(open(0).read().split())
a=int_(fft.irfft(fft.rfft(bincount(a,[1]*n,2**18))**2)+.5)
c=0
for i in where(a>0)[0][::-1]:
  t=min(m,a[i])
  c+=i*t
  m-=t
  if m<1:break
print(c)

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

N=2**18

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

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

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,*a=map(int,open(0).read().split())
b=[0]*N
for i in a:b[i]+=1
i=N
c=0
for a in ifft([t*t for t in fft(b,N)],N)[::-1]:
  a=int(a.real+.5)
  i-=1
  if a:
    t=min(m,a)
    c+=i*t
    m-=t
    if not m:break
print(c)

Solution (Python):
from scipy.signal import*
n,m,*a=map(int,open(0).read().split())
b=[0]*2**17
for i in a:b[i]+=1
i=2**18-1
c=0
for a in(fftconvolve(b,b)+.5).astype(int)[::-1]:
  i-=1
  if a:
    t=min(m,a)
    c+=i*t
    m-=t
    if not m:break
print(c)

Solution (Python):
from numpy import*
n,m,*a=int_(open(0).read().split())
b=[0]*2**18
for i in a:b[i]+=1
i=2**18
c=0
for a in int_(fft.irfft(fft.rfft(b)**2)+.5)[::-1]:
  i-=1
  if a:
    t=min(m,a)
    c+=i*t
    m-=t
    if not m:break
print(c)

Solution (Python):
from itertools import accumulate

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

A.sort(reverse=True)
d = [0] * (A[0] * 2 + 1)
for i in range(N):
    d[A[i]] += 1
d = list(accumulate(d[::-1]))[::-1]
S = list(accumulate([0] + A))

left = 0
right = A[0] * 2 + 1
while right - left > 1:
    mid = (left + right) // 2
    m = 0
    for i in range(N):
        m += d[max(0, mid - A[i])]
    if m >= M:
        left = mid
    else:
        right = mid

ans = 0
m = 0
for i in range(N):
    t = min(d[max(0, right - A[i])], M - m)
    ans += S[t] + A[i] * t
    m += t
ans += left * (M - m)
print(ans)

Solution (Python):






















N, M = [int(x) for x in input().split()]
A = [int(x) for x in input().split()]
A.sort()
A.reverse()

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

def binary_search(l, u):
  m = (l+u)//2
  c, k = PairCount_and_TotalSum_over(m)
  if c<M:
    return binary_search(l, m)
  else:
    if l==m:
      print(k-(c-M)*l)
    else:
      return binary_search(m, u)

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

Solution (Python):
from cmath import pi, exp


def convolve(a, b):
    """配列a, bの畳み込みをFFT(fast Fourier transform)によって計算する
    res[k] = sum{i=0...n}a[i]*b[k-i]
    """
    def fft(a, inverse=False):
        
        for i in range(n):
            j = 0
            for k in range(h):
                j |= (i >> k & 1) << (h - 1 - k)
            if i < j:
                a[i], a[j] = a[j], a[i]

        
        sign = 1.0 if inverse else -1.0
        b = 1
        while b < n:
            for j in range(b):
                w = exp(sign * 2.0j * pi / (2.0 * b) * j)
                for k in range(0, n, b*2):
                    s = a[j + k]
                    t = a[j + k + b] * w
                    a[j + k] = s + t
                    a[j + k + b] = s - t
            b *= 2
 
        
        if inverse:
            for i in range(n):
                a[i] /= n
        return a

    def ifft(a):
        return fft(a, inverse=True)

    n = 1 << (len(a) + len(b) - 1).bit_length()
    h = n.bit_length() - 1
    a += [0] * (n - len(a))
    b += [0] * (n - len(b))
    
    
    fft(a), fft(b)
    for i in range(n):
        a[i] *= b[i]
    ifft(a)

    res = [int(abs(v) + 0.5) for i, v in enumerate(a)]
    return res[1:]
            

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

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

cnt_pair = convolve(cnt[0:], cnt[0:])
len_cnt_pair = len(cnt_pair)
ans = 0
c = 0
for i, num in enumerate(cnt_pair[::-1]):
    val = len_cnt_pair - i
    if num + c <= m:
        ans += num * val
        c += num
    else:
        ans += (m - c) * val
        break
print(ans)

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

n, m = map(int, input().split())
hands = [int(x) for x in input().split()]
increasing = sorted(hands)
hands = list(reversed(increasing))
xs = [0] + list(accumulate(hands))

max_x = 2 * 10**5 + 1
min_x = 0


def canGet(x):
    count = 0
    for hand in increasing:
        idx = bisect_left(increasing, x-hand)
        count += n - idx
    return count >= m


while max_x - min_x > 1:
    
    mid = (max_x + min_x) // 2
    if canGet(mid):
        min_x = mid
    else:
        max_x = mid

ans = 0
count = 0

for Ai in hands:
    idx = bisect_left(increasing, min_x-Ai)
    ans += Ai*(n-idx) + xs[n-idx]
    count += n-idx
print(ans-(count-m)*min_x)

Solution (Python):
n,m=map(int,input().split())
A=[int(i) for i in input().split()]
A.sort()
s=sum(A)
SA=[0]
for a in A:
  SA.append(SA[-1]+a)
for i in range(n+1):
  SA[i]=s-SA[i]
l,r=0,2*max(A)+1
import bisect
def chk(x):
  ct=0
  for a in A:
    ct+=n-bisect.bisect_left(A,max(x-a,0))
  if ct>=m:
    return True
  else:
    return False

def count(x):
  ct=0
  for a in A:
    ct+=n-bisect.bisect_left(A,max(x-a,0))
  return ct

while r-l>1:
  mid=(l+r)//2
  if chk(mid):
    l=mid
  else:
    r=mid
ans=0
for a in A:
  aa=bisect.bisect_left(A,max(l-a,0))
  ans+=SA[aa]+a*(n-aa)
print(ans-l*(count(l)-m))

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






max_p = max(A)

X = max_p



cnt = [0 for x in range(100001)]
hpp = [0 for x in range(100001)]




for i in A:
    
    
    
    cnt[i] += 1
    hpp[i] += i


for i in range(1, 100001)[::-1]:
    cnt[i - 1] += cnt[i]
    hpp[i - 1] += hpp[i]

def countpair(left, x):
    
    if x - left <= 0:
        return N
    elif x - left > 100000:
        return 0
    else:
        return cnt[x - left]

def get_hpp(left, x):
    if x - left <= 0:
        return hpp[0]
    elif x - left > 100000:
        return 0
    else:
        return hpp[x - left]





small = 0
large = max_p * 2
threshold = 0




while (small <= large):
    
    
    X = (small + large) // 2

    nums = 0
    for i in range(N):
        left = A[i]
        
        nums += countpair(left, X)

    if nums >= M:
        threshold = X
        small = X + 1
    else:
        large = X - 1

X = threshold




Y = X + 1

hpp_sum = 0
num_sum = 0

for i in range(N):
    
    left = A[i]
    cnt_num = countpair(left, Y)
    num_sum += cnt_num
    
    
    hpp_sum += left * cnt_num
    
    hpp_sum += get_hpp(left, Y)
    


hpp_sum += X * (M - num_sum) if (M - num_sum > 0) else 0

print(hpp_sum)

Solution (Python):
import sys
from bisect import bisect_left,bisect_right
sys.setrecursionlimit(10**9)
INF=10**18
def input():
    return sys.stdin.readline().rstrip()

def main():
    N,M=map(int,input().split())
    A=sorted(list(map(int,input().split())))
    S=[0]*(N+1)
    for i in range(N):
        S[i+1]=S[i]+A[i]
    def nibutan(ok,ng):
        while abs(ok-ng) > 1:
            mid = (ok + ng) // 2
            if solve(mid):
                ok = mid
            else:
                ng = mid
        return ok
    
    def solve(mid):
        c=0
        for i in range(N):
            c+=N-bisect_left(A,mid-A[i])
        if c>=M:
            return True
        else:
            return False
        
    x=nibutan(0,10**11)
    ans=0
    count=0
    for i in range(N):
        b_l=bisect_left(A,x-A[i])
        count+=(N-b_l)
        ans+=S[N]-S[b_l]+A[i]*(N-b_l)
    if count==M:
        print(ans)
    else:
        print(ans+(M-count)*x)
    

if __name__ == '__main__':
    main()

Solution (Python):
import numpy as np
from numpy.fft import rfft, irfft
N,M = list(map(int,input().split()))
A = list(map(int,input().split()))
X = np.zeros(10**5+1)
l = 2*(10**5)+10
for i in range(N):
    X[A[i]] += 1
data = np.rint(irfft(rfft(X, l)**2 ))
num = 0
ans = 0
for i in range(len(data) -1 ,0,-1):
    if data[i] != 0:
        if num + data[i] >= M:
            ans += (M - num) * i
            break
        else:
            ans += data[i]*i
            num += data[i]
print(int(ans))

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

N,M = map(int,input().split())
A = sorted(list(map(int,input().split())))
A_r = list(reversed(A))
B = [0] + list(accumulate(A_r))

def func(x):
    count = 0
    for Ai in A:
        idx = bisect_left(A,x-Ai)
        count += N-idx
    if count >= M:
        return True
    else:
        return False

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
count = 0

for Ai in A_r:
    idx = bisect_left(A,MIN-Ai)
    ans += Ai*(N-idx) + B[N-idx]
    count += N-idx
print(ans-(count-M)*MIN)

Solution (Python):
import sys
import bisect
import itertools


def solve():
    readline = sys.stdin.buffer.readline
    mod = 10 ** 9 + 7
    n, m = list(map(int, readline().split()))
    a = list(map(int, readline().split()))
    a.sort()

    cor_v = -1
    inc_v = a[-1] * 2 + 1
    while - cor_v + inc_v > 1:
        bin_v = (cor_v + inc_v) // 2
        cost = 0
        
        p = n
        for v in a:
            p = bisect.bisect_left(a, bin_v - v, 0, p)
            cost += n - p
            if cost > m:
                break
        
        if cost >= m:
            cor_v = bin_v
        else:
            inc_v = bin_v
    acm = [0] + list(itertools.accumulate(a))
    c = 0
    t = 0
    for v in a:
        p = bisect.bisect_left(a, cor_v - v)
        c += n - p
        t += acm[-1] - acm[p] + v * (n - p)
    print(t - (c - m) * cor_v)


if __name__ == '__main__':
    solve()

Solution (Python):
import sys
import bisect
import itertools


def solve():
    readline = sys.stdin.buffer.readline
    mod = 10 ** 9 + 7
    n, m = list(map(int, readline().split()))
    a = list(map(int, readline().split()))
    a.sort()

    cor_v = -1
    inc_v = 10 ** 11
    while - cor_v + inc_v > 1:
        bin_v = (cor_v + inc_v) // 2
        cost = 0
        
        for v in a:
            cost += n - bisect.bisect_left(a, bin_v - v)
            if cost > m:
                break
        
        if cost >= m:
            cor_v = bin_v
        else:
            inc_v = bin_v
    acm = [0] + list(itertools.accumulate(a))
    c = 0
    t = 0
    for v in a:
        p = bisect.bisect_left(a, cor_v - v)
        c += n - p
        t += acm[-1] - acm[p] + v * (n - p)
    print(t - (c - m) * cor_v)


if __name__ == '__main__':
    solve()

Solution (Python):
import numpy as np

N ,M= map(int,input().split())
l=list(map(int,input().split()))
n0 = 2**int(np.ceil(np.log2(2*(max(l))-1)))

A = np.zeros(n0)

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

C = np.fft.ifft( np.fft.fft(A)*np.fft.fft(A) )
ans=0
j=2
q=[]
for ci in np.real(C[:2*(max(l))-1] + 0.5):
    q.append([int(ci),j])
    j+=1

q.reverse()

z=0
j=0
while z<M:
    ans+=q[j][0]*q[j][1]
    z+=q[j][0]
    j+=1
ans-=(z-M)*q[j-1][1]
print(ans)

Solution (Python):
import numpy as np

N ,M= map(int,input().split())
l=list(map(int,input().split()))
n0 = 2**int(np.ceil(np.log2(2*(max(l))-1)))

A = np.zeros(n0)



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

C = np.fft.ifft( np.fft.fft(A)*np.fft.fft(A) )
ans=0
j=2
q=[]
for ci in np.real(C[:2*(max(l))-1] + 0.5):
    
    q.append([int(ci),j])
    j+=1
q.reverse()

z=0
j=0
while z<M:
    ans+=q[j][0]*q[j][1]
    z+=q[j][0]
    j+=1
ans-=(z-M)*q[j-1][1]
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):
n, m=map(int, input().split())
a=list(map(int, input().split()))
a.sort()
aa=a[:]
for i in range(n-2, -1, -1):
  aa[i]+=aa[i+1]

g=[0]*(a[-1]+1)
for i in range(n):
  g[a[i]]+=1
for i in range(a[-1]-1, -1, -1):
  g[i]+=g[i+1]
  
def f(index, point):
  tmp=point-a[index]
  if tmp>a[-1]:
    return 0
  elif tmp<=0:
    return n
  else:
    return g[tmp]
  
def ff(point):
  return sum([f(i, point) for i in range(n)])

def fff():
  x=0
  y=2*a[-1]+1
  c=(x+y)//2
  while y-x>1:
    if ff(c)>m:
      x=(x+y)//2
    else:
      y=(x+y)//2
    c=(x+y)//2
  return y

def ffff(index, point):
  tmp=f(index, point)
  if tmp>0:
    return aa[-tmp]+a[index]*tmp
  else:
    return 0
  
p=fff()
q=m-ff(p)
out=0
for i in range(n):
  out+=ffff(i, p)
out+=(p-1)*q
print(out)

Solution (Python):
import bisect
N, M = list(map(int,input().split()))
data = list(map(int,input().split()))
data.sort()
def cost(x):
    re = 0
    for i in range(N):
        idx = bisect.bisect_left(data,x - data[i])
        re += N - idx
    return re

lb = 0
ub = 2*max(data) + 1
while ub - lb > 1:
    mid = (ub+lb)//2
    if cost(mid) >= M:
        lb = mid
    else:
        ub = mid
ans = 0 - (cost(lb) - M)* lb
cum_data = [0]*(N+1)
for i in range(1,N+1):
    cum_data[i] = cum_data[i -1] + data[i-1]
data_idx = [0]*N
for i in range(N):
    idx = bisect.bisect_left(data,lb - data[i])
    data_idx[i] = idx
for i in range(N):
    ans += (N - data_idx[i]) * data[i]
    ans += cum_data[N] - cum_data[data_idx[i]]
print(ans)

Solution (Python):
N,M=map(int,input().split())
A=list(map(int,input().split()))
A.sort(reverse=True)
B=[0]*(N+1)
for i in range(1,N+1):
    B[i]=B[i-1]+A[i-1]
D=[0]*(2*10**5+1)
for i in range(N):
    D[A[i]]+=1
for i in range(len(D)-1,0,-1):
    D[i-1]+=D[i]



l=-1
r=2*10**5+1
while r-l>1:
    m=(l+r)//2
    s=0
    for i in range(N):
        s+=D[max(1,m-A[i])]
    if s>=M:
        l=m
    else:
        r=m
ans=0
s=0
for i in range(N):
    v=max(0,r-A[i])
    t=max(0,min(D[v],M-s))
    ans+=B[t]+t*A[i]
    s+=t
ans+=l*(M-s)
print(ans)

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

def check(X):
    i = N-1
    n = 0
    for ax in A:
        while i>=0 and A[i]+ax<X: i-=1
        n+=i+1
    return n>=M



min_ng = 10**6
max_ok = 0
while min_ng-max_ok>1:
    tgt = (min_ng+max_ok)//2
    if check(tgt):
        max_ok = tgt
    else:
        min_ng = tgt

i = N-1
k = sum(A)
ans = 0
n = 0
for ax in A:
    while i>=0 and A[i]+ax<max_ok:
        k-=A[i]
        i-=1
    ans += k + ax*(i+1)
    n += i+1
print(ans-max_ok*(n-M))

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):
N, M = map(int, input().split())
A = list(map(int, input().split()))

A.sort(reverse=True)
U = A[0] * 2


cgs = [0] * (U + 1)
for a in A:
    cgs[a] += 1
for i in range(U, 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 = A[-1] * 2
ng = U + 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)