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 bisect


def calc_pairs(p, A):
    
    
    pairs = 0
    N = len(A)
    for a in A:
        j = bisect.bisect_left(A, p - a)
        pairs += N - j
    return pairs


def main():
    N, M = list(map(int, input().split(' ')))
    A = list(map(int, input().split(' ')))
    A.sort()
    cumsum_A = [a for a in A]
    for i in range(1, N):
        cumsum_A[i] += cumsum_A[i - 1]
    
    ok, ng = 0, 2 * A[-1] + 1
    while ng - ok > 1:
        mid = (ok + ng) // 2
        if calc_pairs(mid, A) >= M:
            ok = mid
        else:
            ng = mid
    
    ans = 0
    for a in A:
        index = bisect.bisect_left(A, ok - a)
        v = (N - index) * a + cumsum_A[N - 1] - (0 if index == 0 else cumsum_A[index - 1])
        ans += v
    
    pairs = calc_pairs(ok, A)
    ans -= (pairs - M) * ok
    print(ans)


if __name__ == '__main__':
    main()

Solution (Python):
import sys
import bisect

stdin = sys.stdin


def ni(): return int(ns())


def na(): return list(map(int, stdin.readline().split()))


def naa(N): return [na() for _ in range(N)]


def ns(): return stdin.readline().rstrip()  


N, M = na()
A_array = na()
A_array.sort()

hap_min = 0
hap_max = 3 * (10 ** 5)

while(hap_max - hap_min > 1):
    hap = (hap_min + hap_max) // 2
    cnt = 0
    for a in A_array:
        idx = bisect.bisect_left(A_array, hap - a)
        cnt += N - idx
    
    if cnt >= M:
        hap_min = hap
    else:
        hap_max = hap

hap = hap_min

cnt_array = [0] * N
sum_array = [0] * N

for i, a in enumerate(A_array):
    idx = bisect.bisect_left(A_array, hap - a)
    cnt_array[i] = N - idx
    if idx != N:
        sum_array[idx] += 1



ans = 0
cusum = 0
cntsum = 0
for a, c, s in zip(A_array, cnt_array, sum_array):
    cusum += s
    cntsum += c
    ans += a * (cusum + c)

print(ans - (cntsum - M) * hap)

Solution (Python):

8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
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):
from bisect import bisect_left, bisect_right
from itertools import accumulate

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

    As.sort()

    def isOK(x):
        num = 0
        for k, A in enumerate(As, start=1):
            i = bisect_left(As, x-A)
            num += N-i
            if num >= M:
                return True
            elif num + (N-k)*N < M:
                return False
        return num >= M

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


    score = 0
    num = 0
    accAs = list(accumulate([0]+As))
    for A in As:
        i = bisect_right(As, ok-A)
        score += A*(N-i) + accAs[N] - accAs[i]
        num += N-i


    ans = score + ok*(M-num)
    print(ans)


solve()

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

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

    As.sort()

    def isOK(x):
        num = 0
        for k, A in enumerate(As, start=1):
            i = bisect_left(As, x-A)
            num += N-i
            if num >= M:
                return True
            elif num + (N-k)*N < M:
                return False
        return num >= M

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


    score = 0
    num = 0
    accAs = list(accumulate([0]+As))
    for A in As:
        i = bisect_right(As, ok-A)
        score += A*(N-i) + accAs[N] - accAs[i]
        num += N-i


    ans = score + ok*(M-num)
    print(ans)


solve()

Solution (Python):
k = 34
K = 1<<k
nu = lambda L: int("".join([bin(K+a)[-k:] for a in L[::-1]]), 2)
st = lambda n: bin(n)[2:] + "0"
li = lambda s: [int(a, 2) if len(a) else 0 for a in [s[-(i+1)*k-1:-i*k-1] for i in range(len(B)*2-1)]]

N, M = map(int, input().split())
A = [int(a) for a in input().split()]
B = [0] * 100001
for a in A:
    B[a] += 1
C = li(st(nu(B) ** 2))
ans = 0
for i in range(200001)[::-1]:
    a = min(M, C[i])
    M -= a
    ans += a * i
print(ans)

Solution (Python):
import numpy as np
def convolve(A, B):
    
    
    dtype = np.int64
    fft, ifft = np.fft.rfft, np.fft.irfft
    a, b = len(A), len(B)
    if a == b == 1:
        return np.array([A[0]*B[0]])
    n = a+b-1  
    k = 1 << (n-1).bit_length()
    AB = np.zeros((2, k), dtype=dtype)
    AB[0, :a] = A
    AB[1, :b] = B
    return np.rint(ifft(fft(AB[0]) * fft(AB[1]))).astype(np.int64)[:n]

import sys
input = sys.stdin.readline
n,m = map(int, input().split())
a = list(map(int, input().split()))
 
cnt = np.zeros(100001)
for i in a:
    cnt[i] += 1 

c = convolve(cnt,cnt)
ans = 0
for i in range(len(c))[::-1]:
    if c[i] > 0:
      p = min(m,c[i])
      m -= p
      ans += i*p
      if m == 0:
        break
print(ans)

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

mm = 10**10
k = mm.bit_length()
K = 1<<k
nu = lambda L: int("".join([bin(K+a)[-k:] for a in L[::-1]]), 2)
st = lambda n: bin(n)[2:] + "0"
li = lambda s: [int(a, 2) if len(a) else 0 for a in [s[-(i+1)*k-1:-i*k-1] for i in range(200001)]]

n,m = map(int, input().split())
a = [int(i) for i in input().split()]
b = [0]*100001

for i in a:
    b[i] += 1
    
c = li(st(nu(b)*nu(b)))
ans = 0
for i in range(200001)[::-1]:
    if c[i] > 0:
      p = min(m,c[i])
      m -= p
      ans += i*p
      if m == 0:
        break
print(ans)

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:]

import sys
input = sys.stdin.readline
n,m = map(int, input().split())
a = list(map(int, input().split()))
 
cnt = [0] * (100001)
for i in a:
    cnt[i] += 1 

c = convolve(cnt[0:], cnt[0:])
ans = 0
for i in range(len(c))[::-1]:
    if c[i] > 0:
      p = min(m,c[i])
      m -= p
      ans += i*p+p
      if m == 0:
        break
print(ans)

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

mm = 10**10
k = mm.bit_length()
K = 1<<k
nu = lambda L: int("".join([bin(K+a)[-k:] for a in L[::-1]]), 2)
st = lambda n: bin(n)[2:] + "0"
li = lambda s: [int(a, 2) if len(a) else 0 for a in [s[-(i+1)*k-1:-i*k-1] for i in range(200001)]]

n,m = map(int, input().split())
a = [int(i) for i in input().split()]
b = [0]*100001

for i in a:
    b[i] += 1
    
c = li(st(nu(b)*nu(b)))
ans = 0
for i in range(200001)[::-1]:
    if c[i] > 0:
      p = min(m,c[i])
      m -= p
      ans += i*p
      if m == 0:
        break
print(ans)

Solution (Python):
from bisect import bisect_right

N, M = map(int, input().split())
*A, = map(int, input().split())
INF = 10**12
A.sort()
ac = [0]*(N+1)
for i in range(1, N+1):
    ac[i] = ac[i-1]+A[i-1]

def f(x):
    v = 0
    for i in A:
        v += N-bisect_right(A, x-i)
    return v

def g(l):
    ans = INF
    for i in A:
        res = bisect_right(A, l-i)
        if 0<=res<N:
            ans = min(ans, A[res]+i)
    return ans

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

ans = 0
for i in A:
    v1 = N-bisect_right(A, r-i)
    ans += ac[N]-ac[N-v1]+v1*i
ans += (M-f(r))*g(l)

print(ans)

Solution (Python):
from bisect import bisect_left, bisect_right

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







def f(x):
    v = 0
    for i in A:
        v += N-bisect_right(A, x-i)
    return v

def get_nearest(l):
    ans = 10**12
    for i in A:
        res = bisect_right(A, l-i)
        if 0<=res<N:
            ans = min(ans, A[res]+i)
    return ans

l, r = 0, 10**12+1
while r-l>1:
    m = (l+r)//2
    if f(m)>=M:
        l = m
    else:
        r = m

ans = 0
for i in A:
    v1 = N-bisect_right(A, r-i)
    ans += ac[N]-ac[N-v1]+v1*i
ans += (M-f(r))*get_nearest(l)

print(ans)

Solution (Python):
from bisect import bisect_left


def solve(m):
    r = N - 1
    l = 0
    cnt = 0
    while 0 <= r and l < N:
        if A[r] + A[l] < m:
            l += 1
        else:
            cnt += N - l
            r -= 1
    return cnt


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




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

X = ok
amari = float('inf')
cnt = 0
score = 0
for i, Ai in enumerate(A):
    Aj = X - Ai
    j = bisect_left(A, Aj)
    score += S[N] - S[j] + Ai * (N - j)
    cnt += N - j
    if j < N:
        amari = min(amari, A[i] + A[j])

print(score - (cnt - M) * amari)

Solution (Python):
from bisect import bisect_left


def solve(m):
    r = N - 1
    l = 0
    cnt = 0
    while 0 <= r and l < N:
        if A[r] + A[l] < m:
            l += 1
        else:
            cnt += N - l
            r -= 1
    return cnt


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




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

X = ok
amari = float('inf')
cnt = 0
score = 0
for i, Ai in enumerate(A):
    Aj = X - Ai
    j = bisect_left(A, Aj)
    score += S[N] - S[j] + Ai * (N - j)
    cnt += N - j
    if j < N:
        amari = min(amari, A[i] + A[j])

print(score - (cnt - M) * amari)

Solution (Python):
import math

import queue
from collections import deque,defaultdict
import heapq as hpq
from sys import stdin,setrecursionlimit


ipt = stdin.readline
setrecursionlimit(10**7)

def main():
    n,m = map(int,ipt().split())
    a = [int(i) for i in ipt().split()]
    a.sort()
    maa = max(a)

    pn = [0]*(2*maa+1)
    pn[0] = n
    pa = 0
    for i in a:
        for j in range(pa+1,i+1):
            pn[j] = pn[pa]
        pn[i+1] = pn[pa]-1
        pa = i+1

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

    pt = [0]*(2*maa+1)
    pp = maa+1
    for i in a[::-1]:
        for j in range(pp-1,i,-1):
            pt[j] = pt[pp]
        pt[i] = pt[pp]+i
        pp = i
    for j in range(pp-1,-1,-1):
        pt[j] = pt[pp]

    ans = 0
    sm = 0
    for i in a:
        sm += pn[max(0,l-i)]
        ans += pt[max(0,l-i)]+pn[max(0,l-i)]*i
    print(ans-l*(sm-m))
    return

if __name__ == '__main__':
    main()

Solution (Python):

from bisect import bisect_left
from itertools import accumulate

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

def f(x):
    count=0
    for i in A:
        idx=bisect_left(A,x-i)
        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 f(MID):
        MIN=MID
    else:
        MAX=MID

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

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

def expD(t):
  nexp, pexp = {}, {}
  while t != 0:
    pexp[t] = cmath.exp(2j*pi/t)
    nexp[t] = cmath.exp(-2j*pi/t)
    t //= 2
  return nexp, pexp
def powerlimit(x):
  a = 1
  while a < x:
    a *= 2
  return a

N, M = map(int, input().split())
A = list(map(int, input().split()))
m = max(A)
X = [0 for i in range(m+1)]
for i in A:
  X[i] += 1
limit = powerlimit(2*(m+1))
X = X[:] + [0 for _ in range(limit-m-1)]

nexp, pexp = expD(limit)

def fft(f, n):
  hn = n//2
  if n == 1:
    return f
  elif n == 2:
    return [f[0] + f[1], f[0] - f[1]]
  else:
    f0 = [f[2*_] for _ in range(hn)]
    f1 = [f[2*_+1] for _ in range(hn)]
    f0 = fft(f0, hn)
    f1 = fft(f1, hn)
    grow = 1
    fund = nexp[n]
    for i in range(hn):
      right = grow * f1[i]
      f[i] = f0[i] + right
      f[i+hn] = f0[i] - right
      grow *= fund
    return f
def ifft(f, n):
  hn = n//2
  if n == 1:
    return f
  elif n == 2:
    return [f[0] + f[1], f[0] - f[1]]
  else:
    f0 = [f[2*_] for _ in range(hn)]
    f1 = [f[2*_+1] for _ in range(hn)]
    f0 = ifft(f0, hn)
    f1 = ifft(f1, hn)
    grow = 1
    fund = pexp[n]
    for i in range(hn):
      right = grow * f1[i]
      f[i] = f0[i] + right
      f[i+hn] = f0[i] - right
      grow *= fund
    return f
  
F = fft(X, limit)
V = [F[i]*F[i] for i in range(limit)]
E = ifft(V, limit)
for i in range(limit):
  E[i] = int(0.5+E[i].real/limit)
cnt = 0
S = 0
point = limit - 1

while cnt < M:
  Ep = E[point]
  if Ep != 0:
    if cnt + Ep > M:
      S += (M - cnt) * point
      cnt = M
    else:
      S += point * Ep
      cnt += Ep
  point -= 1
print(S)

Solution (Python):

def gcd(x, y):
    k, l = x, y
    if y < 0:
        k, l = -k, -l
    while l:
        k, l = l, k % l
    return k
def extgcd(a, b, c):
    if b == 0:
        return c // a, 0
    if b < 0:
        a, b, c = -a, -b, -c
    x, y, u, v, k, l = 1, 0, 0, 1, a, b
    while l != 0:
        x, y, u, v = u, v, x - u * (k // l), y - v * (k // l)
        k, l = l, k % l
    return c * x, c * y


def IsPrime(n):
    if n == 1:
        return False
    p = 2
    while p * p <= n:
        if n % p == 0:
            return False
        p += 1
    return True


def fracsum(f1, f2, sign = 1):
    fs = [f1[0]*f2[1]+sign*f1[1]*f2[0], f1[1]*f2[1]]
    g = gcd(fs[0], fs[1])
    fs[0] //= g
    fs[1] //= g
    return fs
def fracprod(f1, f2):
    fs = [f1[0]*f2[0], f1[1]*f2[1]]
    g = gcd(fs[0], fs[1])
    fs[0] //= g
    fs[1] //= g
    return fs


def factorization(n):
    D = {}
    a = n
    p = 2
    while a != 1:
        cnt = 0
        while a % p == 0:
            cnt += 1
            a //= p
        if cnt != 0:
            D[p] = cnt
        p += 1
        if p * p > n and a != 1:
            D[a] = 1
            break
    return D

mod = 0 
maxf = 0           

def doubling(n, m, modulo=0):
    y = 1
    base = n
    tmp = m
    while tmp:
        if tmp % 2 == 1:
            y *= base
            if modulo > 0:
                y %= modulo
        base *= base
        if modulo > 0:
            base %= modulo
        tmp //= 2
    return y
def Vanilla(D):
    F = [1]
    for i in D:
        mL = F[:]
        dp = i
        for j in range(D[i]):
            for k in mL:
                F.append(dp*k)
            dp *= i
    return F
def bic(n):
    cnt = 0
    while n % 2 == 0:
        cnt += 1
        n //= 2
    return cnt

def inved(a, modulo=5):
    x, y, u, v, k, l = 1, 0, 0, 1, a, modulo
    while l != 0:
        x, y, u, v = u, v, x - u * (k // l), y - v * (k // l)
        k, l = l, k % l
    return x % modulo
"""D = {i: [] for i in range(27, 33)}
for i in range(27, 33):
    j = 1
    bas = 2**i
    while len(D[i]) < 5:
        if IsPrime(j*bas+1):
            D[i].append(j*bas+1)
        j += 2
print(D)"""


mod = [2013265921, 2281701377, 3892314113] 
primod = [137, 551, 499]
"""for k in range(3):
    S = factorization(mod[k] - 1)
    S = {2: S[2]}
    LL = Vanilla(S)
    print(LL)
    LL.sort(reverse=True)
    for i in range(350, 700):
        flg = (doubling(i, doubling(2, S[2]), mod[k]) == 1)
        for j in range(1, len(LL)):
            flg *= (doubling(i, LL[j], mod[k]) != 1)
        if flg:
            print(mod[k], i)
    print()"""
def powlimit(n):
    y = 1
    cnt = 0
    while y < n:
        y *= 2
        cnt += 1
    return y, cnt

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

mm = max(A)
F = [0 for _ in range(mm+1)]
for i in range(N):
    F[A[i]] += 1
"""N = int(input())
AB = [list(map(int, input().split())) for _ in range(N)]
F = [0] + [AB[i][0] for i in range(N)]
G = [0] + [AB[i][1] for i in range(N)]"""
pl, c = powlimit(2*mm+2)

F = F[:] + [0]*(pl-mm-1)
prim = [[doubling(primod[j], doubling(2, 27-i), mod[j]) for i in range(27)] for j in range(2)]
invp = [[inved(prim[j][i], mod[j]) for i in range(27)] for j in range(2)]


def NTT(f, n, cn, num, inverse=False):
    if n == 1:
        return f
    if n == 2:
        return [(f[0]+f[1])%mod[num], (f[0]-f[1])%mod[num]]
    dn = n // 2
    fe = [f[2*i+0] for i in range(dn)]
    fo = [f[2*i+1] for i in range(dn)]
    fe = NTT(fe, dn, cn-1, num, inverse)
    fo = NTT(fo, dn, cn-1, num, inverse)
    seed = inverse * prim[num][cn] + (1 - inverse) * invp[num][cn]
    grow = 1
    xf = [0 for _ in range(n)]
    for i in range(dn):
        right = fo[i] * grow % mod[num]
        xf[i+0*dn] = (fe[i] + right) % mod[num]
        xf[i+1*dn] = (fe[i] - right) % mod[num]
        grow *= seed
        grow %= mod[num]
    return xf

X = [NTT(F, pl, c, j) for j in range(2)]



Z = [[X[j][i]*X[j][i]%mod[j] for i in range(pl)] for j in range(2)]
for j in range(2):
    Z[j] = NTT(Z[j], pl, c, j, 1)

W = [0 for _ in range(2*mm+1)]
for j in range(2):
    ipl = inved(pl, mod[j])
    for i in range(pl):
        Z[j][i] *= ipl
        Z[j][i] %= mod[j]
for i in range(2*mm+1):
    k1, k2 = extgcd(mod[0], -mod[1], Z[1][i] - Z[0][i])
    W[i] = (mod[0]*k1+Z[0][i]) % (mod[0]*mod[1])

S = 2*max(A)
hp = 0
ccnt = 0
while ccnt < M:
    if ccnt + W[S] < M:
        ccnt += W[S]
        hp += S * W[S]
    else:
        hp += S * (M - ccnt)
        ccnt += M - ccnt
    S -= 1
print(hp)

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

a_cum = [0 for _ in range(n+1)]
for i in range(n):
    a_cum[i+1] = a_cum[i] + a[i]

l, r = 0, 10 ** 6
while r - l > 1:
    x = (l + r) // 2
    cnt = 0

    for i in range(n):
        idx = bisect.bisect_left(a, x-a[i])
        cnt += n - idx

    if cnt >= m:
        l = x

    else:
        r = x

ans = 0
cnt = 0
for i in range(n):
    idx = bisect.bisect_left(a, l-a[i])
    cnt += n - idx
    ans += a_cum[n] - a_cum[idx]
    ans += a[i] * (n - idx)

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

print(ans)

Solution (Python):
import numpy as np

n, m = map(int, input().split())
A = np.array(input().split(), np.int64)
A = np.sort(A)

def is_ok(x):
    cnt = n**2-np.searchsorted(A, x-A, side='left').sum()
    if cnt >= m:
        return True
    else:
        return False

l = 0
r = 2*10**5+1
while l+1 < r:
    c = (l+r)//2
    if is_ok(c):
        l = c
    else:
        r = c

from itertools import accumulate
B = [0]+list(A)
B = list(accumulate(B))
ans = 0
cnt = 0
for a in A:
    j = np.searchsorted(A, l-a, side='left').sum()
    cnt += n-j
    ans += B[-1]-B[j]+a*(n-j)
ans -= (cnt-m)*l
print(ans)

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

import bisect
def is_ok(x):
    cnt = 0
    for i, a in enumerate(A):
        j = bisect.bisect_left(A, x-a)
        cnt += n-j
    if cnt >= m:
        return True
    else:
        return False

l = 0
r = 2*10**5+1
while l+1 < r:
    c = (l+r)//2
    if is_ok(c):
        l = c
    else:
        r = c

from itertools import accumulate
B = [0]+A
B = list(accumulate(B))
ans = 0
cnt = 0
for a in A:
    j = bisect.bisect_left(A, l-a)
    cnt += (n-j)
    ans += B[-1]-B[j]+a*(n-j)
ans -= (cnt-m)*l
print(ans)

Solution (Python):
from itertools import accumulate

n, m = map(int, input().split())
a = list(map(int, input().split()))
MAX = 10 ** 5

scores = [0] * (2 * MAX + 1)
for e in a:
    scores[e] += 1

acc = list(accumulate(scores))


def f(x):
    cnt = 0
    for e in a:
        cnt += n - acc[max(0, x - e - 1)]

    return cnt


l = 0
r = 2 * MAX + 1
while r - l > 1:
    mid = (l + r) // 2
    if f(mid) >= m:
        l = mid
    else:
        r = mid

scores_sum = [i * e for i, e in enumerate(scores)]
acc_sum = list(accumulate(scores_sum))
sm = sum(a)
ans = 0
cnt_sum = 0
for e in a:
    i = max(0, l - e - 1)
    cnt = n - acc[i]
    ans += e * cnt + sm - acc_sum[i]
    cnt_sum += cnt

ans -= (cnt_sum - m) * l

print(ans)

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

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

def binary_search():
    l, r = 0, 10**10+100
    
    while l<=r:
        m = (l+r)//2
        
        if judge(m):
            l = m+1
        else:
            r = m-1
    
    return r

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

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

b = binary_search()
ans = 0
cnt = 0

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

ans += (M-cnt)*b
print(ans)

Solution (Python):
import sys
INF = 1 << 60
MOD = 10**10 + 1
sys.setrecursionlimit(2147483647)
input = lambda:sys.stdin.readline().rstrip()

class NTT(object):
    def __init__(self, prime, root):
        self._prime = prime
        self._root = root

    def _fmt(self, A, inverse = False):
        N = len(A)
        logN = (N - 1).bit_length()
        prime = self._prime
        base = pow(self._root, (prime - 1) // N * (1 - 2 * inverse) % (prime - 1), prime)
        step = N
        for k in range(logN):
            step >>= 1
            w = pow(base, step, prime)
            wj = 1
            nA = [0] * N
            for j in range(1 << k):
                for i in range(1 << logN - k - 1):
                    s, t = i + j * step, i + j * step + (N >> 1)
                    ps, pt = i + j * step * 2, i + j * step * 2 + step
                    nA[s], nA[t] = (A[ps] + A[pt] * wj) % prime, (A[ps] - A[pt] * wj) % prime
                wj = (wj * w) % prime
            A = nA
        return A

    def convolution(self, f, g):
        N = 1 << (len(f) + len(g) - 2).bit_length()
        prime = self._prime
        Ff, Fg = self._fmt(f + [0] * (N - len(f))), self._fmt(g + [0] * (N - len(g)))
        N_inv = pow(N, prime - 2, prime)
        fg = self._fmt([a * b % prime * N_inv % prime for a, b in zip(Ff, Fg)], inverse = True)
        del fg[len(f) + len(g) - 1:]
        return fg

def Garner(R, M):
    n = len(R)
    M.append(MOD)
    coeffs, consts = [1] * (n + 1), [0] * (n + 1)
    for i in range(n):
        a, b, u, v = coeffs[i], M[i], 1, 0
        while b:
            u, v = v, u - a // b * v
            a, b = b, a - a // b * b
        t = (R[i] - consts[i]) * (u % M[i]) % M[i]
        for j in range(i + 1, n + 1):
            consts[j] = (consts[j] + t * coeffs[j]) % M[j]
            coeffs[j] = coeffs[j] * M[i] % M[j]
    M.pop()
    return consts[-1]

class arbitrary_mod_NTT(object):
    def __init__(self, primes, roots):
        self._ntts = [NTT(prime, root) for prime, root in zip(primes, roots)]
        self._primes = primes

    def convolution(self, f, g):
        res = [0] * (len(f) + len(g) - 1)
        fgs = [ntt.convolution(f, g) for ntt in self._ntts]
        for i in range(len(res)):
            res[i] = Garner([fg[i] for fg in fgs], self._primes)
        return res

def resolve():
    n, m = map(int, input().split())
    f = [0] * (100001)
    for a in map(int, input().split()):
        f[a] += 1
    ff = arbitrary_mod_NTT([167772161, 469762049, 1224736769], [3, 3, 3]).convolution(f, f)

    ans = 0
    n = len(ff)
    for i in range(n - 1, -1, -1):
        add = min(m, ff[i])
        m -= add
        ans += i * add

    print(ans)
resolve()

Solution (Python):
import sys; input = sys.stdin.buffer.readline
sys.setrecursionlimit(10**7)
from collections import defaultdict
import bisect

con = 10 ** 9 + 7; INF = float("inf")

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

def Binary_Search(A, N, M):
	
	left = 0
	right = 10 ** 7
	ans = 0
	
	Asum = [0]
	for i in range(N):
		Asum.append(Asum[i] + A[-1 - i])

	leftj = [INF, INF]
	rightj = [0, 0]
	
	while left <= right:
		mid = (left + right) // 2
		var = 0
		happiness = 0
		for i in range(N):
			ind = bisect.bisect_left(A, mid - A[i])
			ind = N - ind
			var += ind
			happiness += ind * A[i] + Asum[ind]

		
		if var == M:
			return happiness
		elif var > M:
			leftj = min(leftj, [var, -mid])
			left = mid + 1
		else:
			ans = max(ans, happiness)
			rightj = max(rightj, [var, -mid])
			right = mid - 1

	
	
	

	ans = ans + (M - rightj[0]) * (-leftj[1])
			
	return ans



def main():
	N, M = getlist()
	A = getlist()
	A.sort()
	ans = Binary_Search(A, N, M)

	print(ans)

if __name__ == '__main__':
	main()

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

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

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

def is_ok(arg):
    c = 0
    f = False
    for i in a:
      if c >= m:
        break
      c += (n-bisect_left(a,arg-i))
    if c >= m:
      f = True
    return f

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

x = bisect_ok(2*a[-1]+1,1)
ans = 0
cnt = 0
for i in a:
  ind = n-bisect_left(a,x-i)
  cnt += ind
  ans += i*ind+(b[-1]-b[n-ind])
if cnt > m:
  ans -= x*(cnt-m)
print(ans)

Solution (Python):
import bisect


def is_ok(a, target, m):
    num = 0
    for val in a:
        index = bisect.bisect_left(a, target - val)
        num += len(a) - index
    if num <= m:
        return True
    else:
        return False


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

a.sort()

ok, ng = 10 ** 10, -1
while ok - ng > 1:
    mid = (ok + ng) // 2
    if is_ok(a, mid, m):
        ok = mid
    else:
        ng = mid

rui = [0] * (n + 1)

for i in range(n):
    rui[i+1] = rui[i] + a[i]

cnt = 0
ret = 0

for val in a:
    index = bisect.bisect_left(a, ok - val)
    num = len(a) - index
    cnt += num
    ret += (num * val) + (rui[-1] - rui[index])

ret += (m - cnt) * ng

print(ret)

Solution (Python):
import sys
from collections import defaultdict
readline = sys.stdin.buffer.readline



def geta(fn=lambda s: s.decode()):
    return map(fn, readline().split())


def gete(fn=lambda s: s.decode()):
    return fn(readline().rstrip())


import numpy as np


def main():
    N, M = geta(int)
    A = np.array(list(geta(int)))
    A.sort()

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

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

    while left + 1 < right:
        middle = (left + right) // 2
        if c(middle) >= M:
            left = middle
        else:
            right = middle

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

    ans = Acum[X].sum() + (A * (N - X)).sum()

    ans += left * (M - (N * N - X.sum()))

    print(ans)


if __name__ == "__main__":
    main()

Solution (Python):
import sys
from collections import defaultdict
import bisect
import itertools
readline = sys.stdin.buffer.readline



def geta(fn=lambda s: s.decode()):
    return map(fn, readline().split())


def gete(fn=lambda s: s.decode()):
    return fn(readline().rstrip())


def main():
    N, M = geta(int)
    A = list(geta(int))
    A.sort()

    def c(x):
        """
        @return  
        """
        ret = 0
        for a in A:
            ret += bisect.bisect_left(A, x - a)
        return N * N - ret

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

    while left + 1 < right:
        middle = (left + right) // 2
        if c(middle) >= M:
            left = middle
        else:
            right = middle

    ans = 0
    Acum = list(itertools.accumulate(A[::-1]))[::-1]
    for a in A:
        idx = bisect.bisect_left(A, right - a)
        if idx < N:
            ans += Acum[idx] + a * (N - idx)
    ans += left * (M - c(right))

    print(ans)


if __name__ == "__main__":
    main()