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):
N,M=map(int,input().split())
A=list(map(int,input().split()))
A.sort()
from bisect import bisect_left
def isok(n):
  sum=0
  for i in range(N):
    sum+=N-bisect_left(A,n-A[i])
  return sum>=M
start,end=0,2*10**5+1
while end-start>1:
  mid=(start+end)//2
  if isok(mid):
    start=mid
  else:
    end=mid
b=[A[0]]
for i in range(1,N):
  b.append(b[-1]+A[i])
l=0
ans=0
bn=b[N-1]
for i in range(N):
  a=bisect_left(A,start-A[i])
  l+=N-a
  if a>0:
    ans+=A[i]*(N-a)+bn-b[a-1]
  else:
    ans+=A[i]*(N-a)+bn
ans-=start*(l-M)
print(ans)

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

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

MAX_VAL = 10 ** 5 + 1
hand = [0] * MAX_VAL
for a in A:
  hand[a] += 1
  
import numpy as np
def convolve(f, g):
    fft_len = 1
    while 2 * fft_len < len(f) + len(g) - 1:
        fft_len *= 2
    fft_len *= 2
    
    Ff = np.fft.rfft(f, fft_len)
    Fg = np.fft.rfft(g, fft_len)

    Fh = Ff * Fg

    h = np.fft.irfft(Fh, fft_len)
    h = np.rint(h).astype(np.int64)

    return h[:len(f) + len(g) - 1]

right_hand = np.array(hand, dtype = int)
left_hand = np.array(hand, dtype = int)

H = convolve(left_hand, right_hand)




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

print(ans)

Solution (Python):
import sys
from itertools import accumulate

read = sys.stdin.read
readline = sys.stdin.readline
readlines = sys.stdin.readlines
sys.setrecursionlimit(10 ** 9)
INF = 1 << 60
MOD = 1000000007


def main():
    N, M, *A = map(int, read().split())

    A.sort(reverse=True)

    counter = [0] * (A[0] + 1)
    for a in A:
        counter[a] += 1

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

    def is_ok(x):
        c = 0
        for a in A:
            if x - a <= A[0]:
                c += counter[max(x - a, 0)]
            pass

        return c <= M

    ok = 2 * A[0] + 1
    ng = 2 * A[-1] - 1
    while ok - ng > 1:
        mid = (ok + ng) // 2
        if is_ok(mid):
            ok = mid
        else:
            ng = mid

    x = ok

    ans = 0
    csum = [0]
    csum.extend(accumulate(A))
    for a in A:
        if x - a <= A[0]:
            c = counter[max(x - a, 0)]
            ans += a * c + csum[c]
            M -= c

    ans += (x - 1) * M

    print(ans)

    return


if __name__ == '__main__':
    main()

Solution (Python):
from bisect import bisect_left

def main():
    n, m = map(int, input().split())
    a = list(map(int, input().split()))
    a.sort()
    L, R = 0, 2 * 10**5 + 1
    while L+1 < R:
        P = (L+R)//2
        cnt = 0
        for v in a:
            x = P - v
            cnt += n - bisect_left(a, x)
        if cnt >= m:
            L = P
        else:
            R = P
    csum = [0]
    for v in a:
        csum.append(v)
    for i in range(n):
        csum[i+1] += csum[i]
    ans = 0
    cnt = 0
    for v in a:
        x = L - v
        idx = bisect_left(a, x)
        cnt += n-idx
        ans += csum[-1] - csum[idx]
        ans += v * (n - idx)
    ans -= (cnt - m) * L
    print(ans)

if __name__ == "__main__":
    main()

Solution (Python):
class NTT:
    def __init__(self, convolution_rank, binary_level):
        self.cr = convolution_rank
        self.bl = binary_level
        self.modulo_list = [1]*convolution_rank
        self.primal_root_list = [1]*convolution_rank
        self.bin_list = [1]*(binary_level+1)
        self.primal_base_matrix = [[0 for _ in range(binary_level+1)] for __ in range(convolution_rank)]
        self.inverse_base_matrix = [[0 for _ in range(binary_level+1)] for __ in range(convolution_rank)]
        for i in range(binary_level):
            self.bin_list[i+1] = self.bin_list[i] * 2
    
    def doubling(self, n, m, modulo=0):
        y = 1
        tmp = m
        bas = n
        while tmp:
            if tmp % 2:
                y *= bas
                if modulo:
                    y %= modulo
            bas *= bas
            if modulo:
                bas %= modulo
            tmp >>= 1
        return y
    
    def powlimit(self, n):
        y = 1
        cnt = 0
        while y < n:
            y *= 2
            cnt += 1
        return y, cnt
    
    def IsPrime(self, num):
        p = 2
        if num <= 1:
            return False
        while p * p <= num:
            if num % p == 0:
                return False
            p += 1
        return True
    
    def extgcd(self, a, b, c):
        if b < 0:
            a, b, c = -a, -b, -c
        x, y, u, v, k, l = 1, 0, 0, 1, a, b
        while l:
            x, y, u, v = u, v, x - u * (k // l), y - v * (k // l)
            k, l = l, k % l
        if c % k:
            return "No Solution"
        return x * c, y * c
    
    def inved(self, a, modulo):
        x, y = self.extgcd(a, modulo, 1)
        return (x+modulo)%modulo
    
    def root_manual(self):
        for i in range(self.cr):
            r = 1
            flg = True
            while flg:
                fflg = True
                for j in range(self.bl):
                    if self.doubling(r, self.bin_list[j], self.modulo_list[i]) == 1:
                        fflg = False
                if self.doubling(r, self.bin_list[-1], self.modulo_list[i]) != 1:
                    fflg = False
                if fflg:
                    flg = False
                else:
                    r += 1
                    if r >= self.modulo_list[i]:
                        break
            self.primal_root_list[i] = r
    
    def make_prime_root(self):
        cnt = 0
        j = 1
        last = self.bin_list[-1]
        while cnt < self.cr:
            if self.IsPrime(j*last+1):
                flg = True
                r = 1
                while flg:
                    fflg = True
                    for i in range(self.bl):
                        if self.doubling(r, self.bin_list[i], j*last+1) == 1:
                            fflg = False
                    if self.doubling(r, last, j*last+1) != 1:
                        fflg = False
                    if fflg:
                        flg = False
                    else:
                        r += 1
                        if r >= j*last+1:
                            break
                if flg==False:
                    self.modulo_list[cnt] = j*last+1
                    self.primal_root_list[cnt] = r
                    cnt += 1
            j += 2
    
    def make_basis(self):
        for i in range(self.cr):
            for j in range(self.bl):
                tmp = self.doubling(2, self.bl-j)
                self.primal_base_matrix[i][j] = self.doubling(self.primal_root_list[i], tmp, self.modulo_list[i])
                self.inverse_base_matrix[i][j] = self.inved(self.primal_base_matrix[i][j], self.modulo_list[i])
    
    def NTT(self, f, n, idx, depth, inverse=False, surface=True):
        mod = self.modulo_list[idx]
        i2 = self.inved(2, mod)
        if inverse:
            seed = self.primal_base_matrix[idx][depth]
        else:
            seed = self.inverse_base_matrix[idx][depth]
        if n == 1:
            return f
        if n == 2:
            xf = [(f[0] + f[1])%mod, (f[0] - f[1])%mod]
            if inverse * surface:
                xf[0] *= i2
                xf[0] %= mod
                xf[1] *= i2
                xf[1] %= mod
            return xf
        hn = n // 2
        fe = [f[2*i+0] for i in range(hn)]
        fo = [f[2*i+1] for i in range(hn)]
        fe = self.NTT(fe, hn, idx, depth-1, inverse, False)
        fo = self.NTT(fo, hn, idx, depth-1, inverse, False)
        xf = [0 for _ in range(n)]
        grow = 1
        for i in range(hn):
            right = (fo[i] * grow) % mod
            xf[i+0*hn] = (fe[i] + right) % mod
            xf[i+1*hn] = (fe[i] - right) % mod
            grow *= seed
            grow %= mod
        if inverse * surface:
            invn = self.inved(n, mod)
            for i in range(n):
                xf[i] *= invn
                xf[i] %= mod
        return xf
    
    def second_NTT(self, f, n, idx, depth, inverse=False):
        res = [0 for _ in range(n)]
        tmp = [0 for _ in range(n)]
        MOD = self.modulo_list[idx]
        ipl = self.inved(n, MOD)
        for i in range(n):
            bas = 1
            pos = 0
            for j in range(depth, 0, -1):
                pos += bas * ((i>>(j-1)) % 2)
                bas *= 2
            res[i] = f[pos]
        for i in range(depth):
            grow = 1
            seed = inverse * self.primal_base_matrix[idx][i+1] + (1 - inverse) * self.inverse_base_matrix[idx][i+1]
            for k in range(1<<i):
                for j in range(1<<(depth-i-1)):
                    tmp[j*(1<<(i+1))+k+0*(1<<i)] = (res[j*(1<<(i+1))+k] + grow * res[j*(1<<(i+1))+k+(1<<i)]) % MOD
                    tmp[j*(1<<(i+1))+k+1*(1<<i)] = (res[j*(1<<(i+1))+k] - grow * res[j*(1<<(i+1))+k+(1<<i)]) % MOD
                grow *= seed
                grow %= MOD
            for j in range(n):
                res[j] = tmp[j]
        if inverse:
            for i in range(n):
                res[i] *= ipl
                res[i] %= MOD
        return res
    
    def CRT(self, num, a_list, m_list):
        r = a_list[0]
        bas = m_list[0]
        x, y = 0, 0
        for i in range(1, num):
            x, y = self.extgcd(bas, -m_list[i], a_list[i]-r)
            r += bas * x
            bas *= m_list[i]
        return r % bas

N, M = map(int, input().split())
ntt = NTT(2, 27)
ntt.make_prime_root()
ntt.make_basis()
A = list(map(int, input().split()))
m = max(A)
pl, cnt = ntt.powlimit(2*m+2)
X = [0 for _ in range(pl)]
for i in range(N):
    X[A[i]] += 1
X1 = ntt.second_NTT(X, pl, 0, cnt)
X2 = ntt.second_NTT(X, pl, 1, cnt)
Z1 = [X1[i]*X1[i]%ntt.modulo_list[0] for i in range(pl)]
Z2 = [X2[i]*X2[i]%ntt.modulo_list[1] for i in range(pl)]
Z1 = ntt.second_NTT(Z1, pl, 0, cnt, True)
Z2 = ntt.second_NTT(Z2, pl, 1, cnt, True)
Z = [ntt.CRT(2, [Z1[i], Z2[i]], ntt.modulo_list) for i in range(pl)]
pnt = 2*m
hap = 0
num = 0
while num < M:
    if Z[pnt] + num < M:
        hap += Z[pnt] * pnt
        num += Z[pnt]
    else:
        hap += pnt * (M - num)
        num = M
    pnt -= 1
print(hap)

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

def I(): return int(input())
def MI(): return map(int, input().split())
def LI(): return list(map(int, input().split()))


"""
「1回の握手でX以上幸福度が上がるならば採用」としたときに，何個採用することになるか？
此の値が初めてM以上になったところのXを利用すれば答えが出る．
例えばM+2開催消してしまったら，2*X引けば良い

「X以上なら採用」は左手ごとに，右手を二分探索

"""
def main():
    import bisect
    
    mod=10**9+7
    N,M=MI()
    A=LI()
    A.sort()
    
    def ch(X):
        temp=0
        
        for i in range(N):
            l=A[i]
            t=X-l
            num=bisect.bisect_left(A,t)
            temp+=N-num
            
        return temp>=M
    
    
    ok=0
    ng=2*(10**5) + 1
    
    while ng - ok > 1:
        med=(ok+ng)//2
        if ch(med):
            ok=med
        else:
            ng=med
            
    
    S=[0]*(N+1)
    for i in range(N):
        S[i+1]=S[i]+A[i]
        
    ans=0
    temp=0
    X=ok
    for i in range(N):
        l=A[i]
        t=X-l
        num=bisect.bisect_left(A,t)
        temp+=N-num
        ans+=S[-1] - S[num] + (N-num)*l
        

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

main()

Solution (Python):
from collections import Counter,defaultdict,deque
from heapq import heappop,heappush
from bisect import bisect_left,bisect_right 
import sys,math,itertools,fractions,pprint
sys.setrecursionlimit(10**8)
mod = 10**9+7
INF = float('inf')
def inp(): return int(sys.stdin.readline())
def inpl(): return list(map(int, sys.stdin.readline().split()))

n,m = inpl()
a = inpl(); a.sort()
ok = a[0] + a[0] - 1
ng = a[-1] + a[-1] + 1
def sol(x):
    res = 0
    for i in range(n):
        y = x-a[i]
        res += n-bisect_left(a,y)
    return True if res >= m else False
while abs(ok-ng) > 1:
    mid = (ok+ng)//2
    if sol(mid):
        ok = mid
    else:
        ng = mid
b = [0]
for i in range(n)[::-1]:
    b.append(b[-1]+a[i])
res = 0
cnt = 0

for x in a:
    ind = n-bisect_left(a,ok-x)
    cnt += ind
    res += b[ind] + x*ind
print(res-(cnt-m)*ok)

Solution (Python):
import sys
import math
import cmath
sys.setrecursionlimit(1000000)

def dft_divide(f, l):
    if l == 1:
        return
    ll = l // 2
    f0, f1 = [0] * ll, [0] * ll
    for i in range(ll):
        f0[i] = f[2 * i]
        f1[i] = f[2 * i + 1]
    dft_divide(f0, ll)
    dft_divide(f1, ll)
    dft_solve(f, f0, f1, l)
    return

def dft_solve(f, f0, f1, l):
    ll = l // 2
    zz = cmath.exp(2j * cmath.pi / l)
    z = 1
    for i in range(l):
        f[i] = f0[i % ll] + z * f1[i % ll]
        z *= zz
    return

def idft_divide(f, l):
    if l == 1:
        return
    ll = l // 2
    f0, f1 = [0] * ll, [0] * ll
    for i in range(ll):
        f0[i] = f[2 * i]
        f1[i] = f[2 * i + 1]
    idft_divide(f0, ll)
    idft_divide(f1, ll)
    idft_solve(f, f0, f1, l)
    return

def idft_solve(f, f0, f1, l):
    ll = l // 2
    zz = cmath.exp(-2j * cmath.pi / l)
    z = 1
    for i in range(l):
        f[i] = f0[i % ll] + z * f1[i % ll]
        z *= zz
    return

def multiply(f, g, ma):
    lf, lg = len(f), len(g)
    m = 2 ** (int(math.log2(lf + lg)) + 1)
    for i in range(m - lf):
        f.append(0)
    for i in range(m - lg):
        g.append(0)
    dft_divide(f, m)
    dft_divide(g, m)
    fg = [0] * m
    for i in range(m):
        fg[i] = f[i] * g[i]
    idft_divide(fg, m)
    c = [0] * (2 * ma + 1)
    for i in range(2 * ma + 1):
        c[i] = round(fg[i].real / m)
    return c

n, m = map(int, input().split())
a = list(map(int, input().split()))
ma = max(a)
happiness1 = [0] * (ma + 1)
happiness2 = [0] * (ma + 1)
for i in a:
    happiness1[i] += 1
    happiness2[i] += 1
c = multiply(happiness1, happiness2, ma)
cnt = 0
ans = 0
for i in range(2 * ma, -1, -1):
    if c[i] + cnt >= m:
        ans += (m - cnt) * i
        cnt += (m - cnt)
    else:
        ans += c[i] * i
        cnt += c[i]
    if cnt == m:
        break
print(ans)

Solution (Python):
import sys
import math
import cmath
sys.setrecursionlimit(1000000)

def dft_divide(f, l):
    if l == 1:
        return
    ll = l // 2
    f0, f1 = [0] * ll, [0] * ll
    for i in range(ll):
        f0[i] = f[2 * i]
        f1[i] = f[2 * i + 1]
    dft_divide(f0, ll)
    dft_divide(f1, ll)
    dft_solve(f, f0, f1, l)
    return

def dft_solve(f, f0, f1, l):
    ll = l // 2
    for i in range(l):
        z = cmath.exp(2j * i * cmath.pi / l)
        f[i] = f0[i % ll] + z * f1[i % ll]
    return

def idft_divide(f, l):
    if l == 1:
        return
    ll = l // 2
    f0, f1 = [0] * ll, [0] * ll
    for i in range(ll):
        f0[i] = f[2 * i]
        f1[i] = f[2 * i + 1]
    idft_divide(f0, ll)
    idft_divide(f1, ll)
    idft_solve(f, f0, f1, l)
    return

def idft_solve(f, f0, f1, l):
    ll = l // 2
    for i in range(l):
        z = cmath.exp(-2j * i * cmath.pi / l)
        f[i] = f0[i % ll] + z * f1[i % ll]
    return

def multiply(f, g, ma):
    lf, lg = len(f), len(g)
    m = 2 ** (int(math.log2(lf + lg)) + 1)
    for i in range(m - lf):
        f.append(0)
    for i in range(m - lg):
        g.append(0)
    dft_divide(f, m)
    dft_divide(g, m)
    fg = [0] * m
    for i in range(m):
        fg[i] = f[i] * g[i]
    idft_divide(fg, m)
    
    c = [0] * (2 * ma + 1)
    for i in range(2 * ma + 1):
        c[i] = round(fg[i].real / m)
    return c

n, m = map(int, input().split())
a = list(map(int, input().split()))
ma = max(a)
happiness1 = [0] * (ma + 1)
happiness2 = [0] * (ma + 1)
for i in a:
    happiness1[i] += 1
    happiness2[i] += 1
c = multiply(happiness1, happiness2, ma)
cnt = 0
ans = 0
for i in range(2 * ma, -1, -1):
    if c[i] + cnt >= m:
        ans += (m - cnt) * i
        cnt += (m - cnt)
    else:
        ans += c[i] * i
        cnt += c[i]
    if cnt == m:
        break
print(ans)

Solution (Python):
import sys
import math
import cmath
sys.setrecursionlimit(1000000)

def dft_divide(f, l):
    if l == 1:
        return
    ll = l // 2
    f0 = [f[2 * i] for i in range(ll)]
    f1 = [f[2 * i + 1] for i in range(ll)]
    dft_divide(f0, ll)
    dft_divide(f1, ll)
    dft_solve(f, f0, f1, l)
    return

def dft_solve(f, f0, f1, l):
    ll = l // 2
    for i in range(l):
        z = cmath.exp(2j * i * cmath.pi / l)
        f[i] = f0[i % ll] + z * f1[i % ll]
    return

def idft_divide(f, l):
    if l == 1:
        return
    ll = l // 2
    f0 = [f[2 * i] for i in range(ll)]
    f1 = [f[2 * i + 1] for i in range(ll)]
    idft_divide(f0, ll)
    idft_divide(f1, ll)
    idft_solve(f, f0, f1, l)
    return

def idft_solve(f, f0, f1, l):
    ll = l // 2
    for i in range(l):
        z = cmath.exp(-2j * i * cmath.pi / l)
        f[i] = f0[i % ll] + z * f1[i % ll]
    return

def multiply(f, g):
    lf, lg = len(f), len(g)
    m = 2 ** (int(math.log2(lf + lg)) + 1)
    for i in range(m - lf):
        f.append(0)
    for i in range(m - lg):
        g.append(0)
    dft_divide(f, m)
    dft_divide(g, m)
    fg = [f[i] * g[i] for i in range(m)]
    idft_divide(fg, m)
    c = [round(fg[i].real / m) for i in range(m)]
    return c, m

n, m = map(int, input().split())
a = list(map(int, input().split()))
ma = max(a)
happiness1 = [0] * (ma + 1)
happiness2 = [0] * (ma + 1)
for i in a:
    happiness1[i] += 1
    happiness2[i] += 1
c, lc = multiply(happiness1, happiness2)
cnt = 0
ans = 0
for i in range(min(200000, lc - 1), -1, -1):
    if c[i] + cnt >= m:
        ans += (m - cnt) * i
        cnt += (m - cnt)
    else:
        ans += c[i] * i
        cnt += c[i]
    if cnt == m:
        break
print(ans)

Solution (Python):
from bisect import bisect_left, bisect_right

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

ok, ng = 0, 10**18+1
while abs(ok-ng) > 1:
    mid = (ok+ng)//2
    ok_flag = True
    cnt = 0
    
    
    for i,a in enumerate(al):
        rem = mid-a
        cnt_a = n-bisect_left(al,rem)
        
        cnt += cnt_a
        
    if cnt >= m:
        ok = mid
    else:
    	ng = mid

min_sum = ok
cum_sums = [0]
csum = 0
for a in al:
    csum += a
    cum_sums.append(csum)
    


ans = 0
cnt = 0
for i,a in enumerate(al):
    rem = min_sum - a
    ind = bisect_left(al, rem)
    
    
    csum = cum_sums[n] - cum_sums[ind]
    
    ans += csum
    ans += a*(n-ind)
    cnt += (n-ind)

ans -= (cnt-m)*min_sum

print(ans)

Solution (Python):
from bisect import bisect_left
n, m = map(int, input().split())
import math
import numpy as np
import decimal
import collections
import itertools
import sys
def find(x):
    if par[x] < 0:
        return x
    else:
        par[x] = find(par[x])
        return par[x]
def unite(x, y):
    p = find(x)
    q = find(y)
    if p == q:
        return None
    if p > q:
        p,q = q,p
    par[p] += par[q]
    par[q] = p
def same(x, y):
    return find(x) == find(y)
def size(x):
    return -par[find(x)]
par = [-1 for i in range(n)]
def prime_numbers(x):
    if x < 2:
        return []
    prime_numbers = [i for i in range(x)]
    prime_numbers[1] = 0
    for prime_number in prime_numbers:
        if prime_number > math.sqrt(x):
            break
        if prime_number == 0:
            continue
        for composite_number in range(2 * prime_number, x, prime_number):
            prime_numbers[composite_number] = 0
    return [prime_number for prime_number in prime_numbers if prime_number != 0]
def is_prime(x):
    if x < 2:
        return False
    if x == 2 or x == 3 or x == 5:
        return True
    if x % 2 == 0 or x % 3 == 0 or x % 5 == 0:
        return False
    prime_number = 7
    difference = 4
    while prime_number <= math.sqrt(x):
        if x % prime_number == 0:
            return False
        prime_number += difference
        difference = 6 - difference
    return True
BIT = [0] * (n + 1)
def add(i, x):
    while i <= n:
        BIT[i] += x
        i += i & -i
def query(i):
    s = 0
    while i > 0:
        s += BIT[i]
        i -= i & -i
    return s
a = list(map(int, input().split()))
a.sort(reverse=True)
l, r = 0, 10 ** 6
while r - l > 1:
    m2 = (l + r) // 2
    cnt = 0
    ind = n - 1
    for i in a:
        while ind >= 0 and a[ind] + i < m2:
            ind -= 1
        cnt += ind + 1
    if cnt >= m:
        l = m2
    else:
        r = m2
cum = [0]
for i in a:
    cum.append(cum[-1] + i)
ans = 0
cnt2 = 0
ind2 = n - 1
b = l
for i in a:
    while ind2 >= 0 and a[ind2] + i < b:
        ind2 -= 1
    cnt2 += ind2 + 1
    ans += (ind2 + 1) * i + cum[ind2 + 1]
print(ans - (cnt2 - m) * b)

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

def count_bigger_x(x):
    cnt = 0
    piv = N-1
    for a in A:
        while piv >= 0 and a + A[piv] >= x:
            piv -= 1
        cnt += N - piv - 1
    return cnt


def is_ok(x):
    cnt = count_bigger_x(x)
    return cnt >= M


def bisect(ng, ok):
    while (abs(ok - ng) > 1):
        mid = (ok + ng) // 2
        if is_ok(mid):
            ok = mid
        else:
            ng = mid
    return ok

x = bisect(2 * A[-1] + 1, 2 * A[0] - 1)
ans = 0
piv = N-1
for a in A:
    while piv >= 0 and a + A[piv] >= x:
        piv -= 1
    ans += (N - piv - 1) * a
ans *= 2
ans -= 1 * x * (count_bigger_x(x) - M)
print(ans)

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

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

A.sort(reverse=True)
l = 0
r = 3*10**5
while r - l > 1:
    m = (r + l)//2 
    ind = N-1
    cnt = 0
    for a in A:
        while ind >= 0 and A[ind] + a < m:
            ind -= 1
        cnt += ind+1
    if cnt >= M:
        l = m
    else:
        r = m

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

ans = 0
border = l
ind = N-1
cnt = 0
for a in A:
    while ind >= 0 and A[ind] + a < border:
        ind -= 1
    cnt += (ind+1)
    ans += a*(ind+1) + B[ind+1]

ans -= (cnt-M)*border

print(ans)

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

sys.setrecursionlimit(10 ** 7)
input = sys.stdin.readline
f_inf = float('inf')
mod = 10 ** 9 + 7


def resolve():
    def meguru_bisect(ok, ng):
        while abs(ok - ng) > 1:
            mid = (ok + ng) // 2
            if is_ok(mid):
                ok = mid
            else:
                ng = mid
        return ok

    def is_ok(x):
        cnt = 0
        for a in A:
            t = x - a
            idx = bisect_right(A, t)
            cnt += n - idx
        return cnt < m

    n, m = map(int, input().split())
    A = list(map(int, input().split()))
    A.sort()
    ng, ok = 0, 10 ** 15 + 1
    mth = meguru_bisect(ok, ng)
    R = [0] + list(accumulate(A))

    res = 0
    cnt = 0
    for a in A:
        s = mth - a
        left = bisect_left(A, s)
        res += (n - left) * a + R[-1] - R[left]
        cnt += n - left
    print(res - mth * (cnt - m))


if __name__ == '__main__':
    resolve()

Solution (Python):
import sys,queue,math,copy,itertools,bisect,collections,heapq

def main():
    LI = lambda : [int(x) for x in sys.stdin.readline().split()]

    N,M = LI()
    A = LI()
    a = sorted(A,reverse=True)

    high = a[0] * 2 + 1
    low = 0
    while high > low + 1:
        mid = (high + low) // 2

        cnt = 0
        j = N - 1
        for i in range(N):
            while a[i] + a[j] < mid:
                j -= 1
                if j < 0: break
            else:
                cnt += (j+1)
                continue
            break
        if cnt < M:
            high = mid
        else:
            low = mid

    s = [a[0]]
    for i in range(1,N):
        s.append(s[-1]+a[i])

    ans = 0
    cnt = 0
    j = N - 1
    for i in range(N):
        while a[i] + a[j] < low:
            j -= 1
            if j < 0: break
        else:
            ans += a[i] * (j+1) + s[j]
            cnt += j + 1
            continue
        break

    print(ans - (cnt-M) * low)


if __name__ == '__main__':
    main()

Solution (Python):
import sys,queue,math,copy,itertools,bisect,collections,heapq

def main():
    LI = lambda : [int(x) for x in sys.stdin.readline().split()]

    N,M = LI()
    A = LI()
    a = sorted(A,reverse=True)

    high = a[0] * 2 + 1
    low = 0
    while high > low + 1:
        mid = (high + low) // 2

        cnt = 0
        j = N - 1
        for i in range(N):
            while a[i] + a[j] < mid:
                j -= 1
                if j < 0: break
            else:
                cnt += (j+1)
                continue
            break
        if cnt <= M:
            high = mid
        else:
            low = mid

    s = [a[0]]
    for i in range(1,N):
        s.append(s[-1]+a[i])

    ans = 0
    cnt = 0
    m = a[0] * 2
    j = N - 1
    for i in range(N):
        while a[i] + a[j] < low:
            j -= 1
            if j < 0: break
        else:
            ans += a[i] * (j+1) + s[j]
            cnt += j + 1
            m = min(m,a[i] + a[j])
            continue
        break

    print(ans - (cnt-M) * low)


if __name__ == '__main__':
    main()

Solution (Python):
import sys,queue,math,copy,itertools,bisect,collections,heapq

def main():
    LI = lambda : [int(x) for x in sys.stdin.readline().split()]

    N,M = LI()
    A = LI()
    a = sorted(A,reverse=True)

    high = a[0] * 2
    low = 0
    while high > low + 1:
        mid = (high + low) // 2

        cnt = 0
        j = N - 1
        for i in range(N):
            while a[i] + a[j] < mid:
                j -= 1
                if j < 0: break
            else:
                cnt += (j+1)
                continue
            break
        if cnt <= M:
            high = mid
        else:
            low = mid

    s = [a[0]]
    for i in range(1,N):
        s.append(s[-1]+a[i])

    ans = 0
    cnt = 0
    m = a[0] * 2
    j = N - 1
    for i in range(N):
        while a[i] + a[j] < low:
            j -= 1
            if j < 0: break
        else:
            ans += a[i] * (j+1) + s[j]
            cnt += j + 1
            m = min(m,a[i] + a[j])
            continue
        break

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


if __name__ == '__main__':
    main()

Solution (Python):
import numpy as np
from numpy.fft import rfft,irfft
import sys
sys.setrecursionlimit(10**9)

def mi(): return map(int,input().split())
def ii(): return int(input())
def isp(): return input().split()
def deb(text): print("-------\n{}\n-------".format(text))

INF=10**20
class Polynomial:
    def __init__(self,coefficient=None,dim=0,const=1):
        if coefficient == None:
            self.coefficient = np.zeros(dim+1, np.int64)
            self.coefficient[dim] = const
        else:
            self.coefficient = coefficient

    def __add__(self,other): 
        f,g = self.coefficient, other.coefficient
        if len(f) > len(g): f,g = g,f

        h = Polynomial(dim=len(g)-1,const=0)
        h.coefficient[len(f):] += g[len(f):]
        h.coefficient[:len(f)] += f + g[:len(f)]
        return h
    
    def __iadd__(self,other): 
        h = self.__add__(other)
        self.coefficient = h.coefficient
        return self

    def __mul__(self,other): 
        f = self.coefficient
        g = other.coefficient
        h = Polynomial()
        h.coefficient = self.fft(f,g)[:len(f)+len(g)-1]
        return h
    
    def __len__(self):
        return len(self.coefficient)

    def __getitem__(self,key):
        return self.coefficient[key]

    def get_coefficient(self,x):
        return self.coefficient[x]


    def fft(self,A,B,fft_len=1<<18):
        x = irfft(rfft(A,fft_len) * rfft(B,fft_len))
        return np.rint(x).astype(np.int64)



def main():
    N,M=mi()
    A=list(mi())
    
    coefficient = [0] * (2*10**5+10)
    for a in A:
        coefficient[a] += 1
    
    f = Polynomial(coefficient=coefficient)
            
    g = f * f
    ans = 0
    handshake = M
    for x in range(len(g)-1,1,-1):
        count = g[x]
        if count <= 0: continue

        if count > handshake:
            ans += handshake * x
        else:
            ans += count * x
 
        handshake = max(handshake-count,0)
        if handshake == 0:
            break

    print(ans)

if __name__ == "__main__":
    main()