Problem p02821 (python)

Fixed submission:
# ----------- gcd ----------- #
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

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

# --- fraction calculator --- #
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

# ------ factorization ------ #
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
# --------------fact. calculator-------------- #
mod = 0 # <-- input modulo
maxf = 0           # <-- input factional limitation

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 = 40961
#mod = 998244353, primod = 50475, invmod = 305218526
mod = [2013265921, 2281701377, 3892314113] # 998244353 = 2^23 * (odd number)
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()))
#A = [1 for i in range(N)]
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)
#pl, c = powlimit(2*N+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)]
#print(prim)
#print(invp)
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)]
#Y = NTT(G, pl, c)
#print(X)
#print(Y)
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)
#print(NTT(Z, pl, c, 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])
#print(W)
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)
#for i in range(1, max(A)+1):
#    print(Z[i]*ipl%mod)

Buggy submission (verdict: Runtime Error):
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):
        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)+1
pl, cnt = ntt.powlimit(2*m)
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)

Buggy submission (verdict: Runtime Error):
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):
        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(1, 28)
ntt.make_prime_root()
ntt.make_basis()
A = list(map(int, input().split()))
m = max(A)+1
pl, cnt = ntt.powlimit(2*m)
X = [0 for _ in range(pl)]
for i in range(N):
    X[A[i]] += 1
X = ntt.second_NTT(X, pl, 0, cnt)
Z = [X[i]*X[i]%ntt.modulo_list[0] for i in range(pl)]
Z = ntt.second_NTT(Z, pl, 0, cnt, True)
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)

Buggy submission (verdict: Runtime Error):
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):
        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(1, 27)
ntt.make_prime_root()
ntt.make_basis()
A = list(map(int, input().split()))
m = max(A)+1
pl, cnt = ntt.powlimit(2*m)
X = [0 for _ in range(pl)]
for i in range(N):
    X[A[i]] += 1
X = ntt.second_NTT(X, pl, 0, cnt)
Z = [X[i]*X[i]%ntt.modulo_list[0] for i in range(pl)]
Z = ntt.second_NTT(Z, pl, 0, cnt, True)
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)

Buggy submission (verdict: Runtime Error):
# ----------- gcd ----------- #
def gcd(x, y):
    k, l = x, y
    if y < 0:
        k, l = -k, -l
    while l:
        k, l = l, k % l
    return k

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

# --- fraction calculator --- #
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

# ------ factorization ------ #
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
# --------------fact. calculator-------------- #
mod = 0 # <-- input modulo
maxf = 0           # <-- input factional limitation

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
mod = 2013265921
def inved(a, modulo=mod):
    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(39, 40)}
for i in range(39, 40):
    j = 1
    bas = 2**i
    while len(D[i]) < 1:
        if IsPrime(j*bas+1):
            D[i].append(j*bas+1)
        j += 2
print(D)"""
#mod = 40961
#mod = 998244353, primod = 50475, invmod = 305218526
mod = 2013265921 # 998244353 = 2^23 * (odd number)
primod = 137
"""S = factorization(mod - 1)
S = {2: S[2]}
LL = Vanilla(S)
print(LL)
LL.sort(reverse=True)
for i in range(1, 500):
    flg = (doubling(i, doubling(2, S[2]), mod) == 1)
    for j in range(1, len(LL)):
        flg *= (doubling(i, LL[j], mod) != 1)
    if flg:
        print(i)"""
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)
#pl, c = powlimit(2*N+2)
F = F[:] + [0]*(pl-mm-1)
prim = [doubling(primod, doubling(2, 27-i), mod) for i in range(27)]
invp = [inved(prim[i]) for i in range(27)]
def NTT(f, n, cn, inverse=False):
    if n == 1:
        return f
    if n == 2:
        return [(f[0]+f[1])%mod, (f[0]-f[1])%mod]
    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, inverse)
    fo = NTT(fo, dn, cn-1, inverse)
    seed = inverse * prim[cn] + (1 - inverse) * invp[cn]
    grow = 1
    for i in range(dn):
        right = fo[i] * grow % mod
        f[i+0*dn] = (fe[i] + right) % mod
        f[i+1*dn] = (fe[i] - right) % mod
        grow *= seed
        grow %= mod
    return f

X = NTT(F, pl, c)
#Y = NTT(G, pl, c)
#print(X)
#print(Y)
Z = [X[i]*X[i]%mod for i in range(pl)]
Z = NTT(Z, pl, c, 1)
#print(NTT(Z, pl, c, 1))
ipl = inved(pl)
for i in range(pl):
    Z[i] *= ipl
    Z[i] %= mod
S = 2*max(A)
hp = 0
ccnt = 0
while ccnt < M:
    if ccnt + Z[S] < M:
        ccnt += Z[S]
        hp += S * Z[S]
    else:
        hp += S * (M - ccnt)
        ccnt += M - ccnt
    S -= 1
print(hp)
#for i in range(1, max(A)+1):
#    print(Z[i]*ipl%mod)

Buggy submission (verdict: Runtime Error):
# ----------- gcd ----------- #
def gcd(x, y):
    k, l = x, y
    if y < 0:
        k, l = -k, -l
    while l:
        k, l = l, k % l
    return k

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

# --- fraction calculator --- #
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

# ------ factorization ------ #
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
# --------------fact. calculator-------------- #
mod = 0 # <-- input modulo
maxf = 0           # <-- input factional limitation

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
mod = 2013265921
def inved(a, modulo=mod):
    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(39, 40)}
for i in range(39, 40):
    j = 1
    bas = 2**i
    while len(D[i]) < 1:
        if IsPrime(j*bas+1):
            D[i].append(j*bas+1)
        j += 2
print(D)"""
#mod = 40961
#mod = 998244353, primod = 50475, invmod = 305218526
mod = 2013265921 # 998244353 = 2^23 * (odd number)
primod = 137
"""S = factorization(mod - 1)
S = {2: S[2]}
LL = Vanilla(S)
print(LL)
LL.sort(reverse=True)
for i in range(1, 500):
    flg = (doubling(i, doubling(2, S[2]), mod) == 1)
    for j in range(1, len(LL)):
        flg *= (doubling(i, LL[j], mod) != 1)
    if flg:
        print(i)"""
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()))
F = [0 for _ in range(max(A)+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*max(A)+2)
#pl, c = powlimit(2*N+2)
F = F[:] + [0]*(pl-N-1)
prim = [doubling(primod, doubling(2, 27-i), mod) for i in range(27)]
invp = [inved(prim[i]) for i in range(27)]
def NTT(f, n, cn, inverse=False):
    if n == 1:
        return f
    if n == 2:
        return [(f[0]+f[1])%mod, (f[0]-f[1])%mod]
    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, inverse)
    fo = NTT(fo, dn, cn-1, inverse)
    seed = inverse * prim[cn] + (1 - inverse) * invp[cn]
    grow = 1
    for i in range(dn):
        right = fo[i] * grow % mod
        f[i+0*dn] = (fe[i] + right) % mod
        f[i+1*dn] = (fe[i] - right) % mod
        grow *= seed
        grow %= mod
    return f

X = NTT(F, pl, c)
#Y = NTT(G, pl, c)
#print(X)
#print(Y)
Z = [X[i]*X[i]%mod for i in range(pl)]
Z = NTT(Z, pl, c, 1)
#print(NTT(Z, pl, c, 1))
ipl = inved(pl)
for i in range(pl):
    Z[i] *= ipl
    Z[i] %= mod
S = 2*max(A)
hp = 0
ccnt = 0
while ccnt < M:
    if ccnt + Z[S] < M:
        ccnt += Z[S]
        hp += S * Z[S]
    else:
        hp += S * (M - ccnt)
        ccnt += M - ccnt
    S -= 1
print(hp)
#for i in range(1, max(A)+1):
#    print(Z[i]*ipl%mod)

Buggy submission (verdict: Runtime Error):
#include <iostream> // cout, endl, cin
#include <string> // string, to_string, stoi
#include <vector> // vector
#include <algorithm> // min, max, swap, sort, reverse, lower_bound, upper_bound
#include <utility> // pair, make_pair
#include <tuple> // tuple, make_tuple
#include <cstdint> // int64_t, int*_t
#include <cstdio> // printf
#include <map> // map
#include <queue> // queue, priority_queue
#include <set> // set
#include <stack> // stack
#include <deque> // deque
#include <unordered_map> // unordered_map
#include <unordered_set> // unordered_set
#include <bitset> // bitset
#include <cctype> // isupper, islower, isdigit, toupper, tolower
#include <cmath>
#include <complex>
using namespace std;
const double pi = 3.14159265358979323846264338;
long long easy_factor(long long n){
  if (n == 1){
    return 1;
  }
  else{
    long long p = 2;
    long long f = 2;
    while (n % p != 0){
      p += 1;
      if (n % p == 0){
        f = p;
        break;
      }
      if (p * p > n){
        f = n;
        break;
      }
    }
    return f;
  }
}

template <typename T>
T vecout(T A){
  for (int i=0;i<A.size();i++){
    cout << A.at(i) << endl;
  }
}

vector<complex<long double>> simply_fft(vector<complex<long double>> f, long long n){
  if (n == 1){
    return f;
  }
  else{
    long long dn = n / 2;
    vector<complex<long double>> fe(dn);
    vector<complex<long double>> fo(dn);
    for (long long j=0;j<dn;j++){
      fe.at(j) = f.at(2*j+0);
      fo.at(j) = f.at(2*j+1);
    }
    vector<complex<long double>> tmp0 = simply_fft(fe, dn);
    vector<complex<long double>> tmp1 = simply_fft(fo, dn);
    complex<long double> grow = 1.0 + 0.0i;
    complex<long double> seed = cos((2.0*pi)/n) - 1i*sin((2.0*pi)/n);
    complex<long double> right = 0;
    for (long long i=0;i<dn;i++){
      right = tmp1.at(i) * grow;
      f.at(i) = tmp0.at(i) + right;
      f.at(i+dn) = tmp0.at(i) - right;
      grow *= seed;
    }
    return f;
  }
}

vector<complex<long double>> simply_ifft(vector<complex<long double>> f, long long n){
  if (n == 1){
    return f;
  }
  else{
    long long dn = n / 2;
    vector<complex<long double>> fe(dn);
    vector<complex<long double>> fo(dn);
    for (long long j=0;j<dn;j++){
      fe.at(j) = f.at(2*j+0);
      fo.at(j) = f.at(2*j+1);
    }
    vector<complex<long double>> tmp0 = simply_ifft(fe, dn);
    vector<complex<long double>> tmp1 = simply_ifft(fo, dn);
    complex<long double> grow = 1.0 + 0.0i;
    complex<long double> seed = cos((2.0*pi)/n) + 1i*sin((2.0*pi)/n);
    complex<long double> right = 0;
    for (long long i=0;i<dn;i++){
      right = tmp1.at(i) * grow;
      f.at(i) = tmp0.at(i) + right;
      f.at(i+dn) = tmp0.at(i) - right;
      grow *= seed;
    }
    return f;
  }
}
vector<long long> vecmachine(){
  vector<long long> A(100000);
  for (int i=0;i<100000;i++){
    A.at(i) = 100000;
  }
  return A;
}
long long powlimit(long long x){
  long long y = 1;
  while (y < x){
    y *= 2;
  }
  return y;
}

// ------------------------------------------

int main(){
  long long N, M;
  cin >> N >> M;
  long long m = 0;
  vector<long long> A = vecmachine();
  for (int i=0;i<N;i++){
    //cin >> A.at(i);
    if (m < A.at(i)){
      m = A.at(i);
    }
  }
  long long pl = powlimit(2*m+2);
  vector<complex<long double>> X(pl, 0);
  for (int i=0;i<N;i++){
    X.at(A.at(i)) += 1.0;
  }
  vector<complex<long double>> F = simply_fft(X, pl);
  vector<complex<long double>> Y(pl, 0);
  for (int i=0;i<pl;i++){
    Y.at(i) = F.at(i) * F.at(i);
  }
  vector<complex<long double>> Z = simply_ifft(Y, pl);
  vector<long long> sol(pl, 0);
  for (long long i=0;i<pl;i++){
    sol.at(i) = floor(0.5+real(Z.at(i))/pl);
  }
  long long S = 0, i = pl - 1, pat = 0;
  while (pat < M){
    if (pat + sol.at(i) > M){
      S += (M - pat) * i;
      pat = M;
    }
    else{
      S += sol.at(i) * i;
      pat += sol.at(i);
    }
    i -= 1;
    if (i < 0){
      break;
    }
  }
  cout << S << endl;
}