Problem:
 Score : 800 points    Problem Statement Let’s take a prime P = 200\,003 .
You are given N integers A_1, A_2, \ldots, A_N .
Find the sum of ((A_i \cdot A_j) \bmod P) over all N \cdot (N-1) / 2 unordered pairs of elements ( i < j ).  Please note that the sum isn't computed modulo P .      Constraints  2 \leq N \leq 200\,000  0 \leq A_i < P = 200\,003  All values in input are integers.         Input Input is given from Standard Input in the following format.  N  A_1  A_2  \cdots  A_N       Output Print one integer — the sum over ((A_i \cdot A_j) \bmod P) .        Sample Input 1 4
2019 0 2020 200002      Sample Output 1 474287  The non-zero products are:   2019 \cdot 2020 \bmod P = 78320  2019 \cdot 200002 \bmod P = 197984  2020 \cdot 200002 \bmod P = 197983   So the answer is 0 + 78320 + 197984 + 0 + 0 + 197983 = 474287 .       Sample Input 2 5
1 1 2 2 100000      Sample Output 2 600013

Solution (Python):
import sys
import numpy as np
import numba
from numba import njit, b1, i4, i8, f8

read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

MOD = 200_003
R = 2  

ntt_mod = 998_244_353
ntt_root = 3
max_ntt = 20
roots = np.empty(1 << max_ntt, np.int64)
iroots = np.empty(1 << max_ntt, np.int64)

ntt_data = (ntt_mod, ntt_root, max_ntt, roots, iroots)
ntt_t = numba.typeof(ntt_data)


@njit
def mpow(a, n):
    p = 1
    while n:
        if n & 1:
            p = p * a % ntt_mod
        a = a * a % ntt_mod
        n >>= 1
    return p


@njit((ntt_t, ), cache=True)
def ntt_precompute(ntt_data):
    MOD, root, max_ntt, roots, iroots = ntt_data

    roots[0] = iroots[0] = 1
    for n in range(max_ntt):
        x = mpow(root, (MOD - 1) >> n + 2)
        y = mpow(x, MOD - 2)
        for i in range(1 << n):
            roots[(1 << n) | i] = roots[i] * x % MOD
            iroots[(1 << n) | i] = iroots[i] * y % MOD
    return roots, iroots


@njit((ntt_t, i8[:], b1), cache=True)
def ntt(ntt_data, A, inverse):
    MOD, _, _, roots, iroots = ntt_data

    N = len(A)

    if not inverse:
        m = N >> 1
        while m:
            for k in range(N // m >> 1):
                s = 2 * m * k
                for i in range(s, s + m):
                    x, y = A[i], A[i + m] * roots[k]
                    A[i], A[i + m] = x + y, x - y
            A %= MOD
            m >>= 1
    else:
        m = 1
        while m < N:
            for k in range(N // m >> 1):
                s = 2 * m * k
                for i in range(s, s + m):
                    x, y = A[i], A[i + m]
                    A[i], A[i + m] = x + y, (x - y) * iroots[k]
            A %= MOD
            m <<= 1
        for i in range(N):
            A[i] = A[i] * mpow(N, MOD - 2) % MOD
    return A


@njit((ntt_t, i8[:], i8[:]), cache=True)
def ntt_convolve(ntt_data, A, B):
    L = 1
    LA, LB = len(A), len(B)
    while 2 * L < LA + LB - 1:
        L *= 2
    L *= 2
    A = np.append(A, np.zeros(L - LA, np.int64))
    B = np.append(B, np.zeros(L - LB, np.int64))
    A = ntt(ntt_data, A, False)
    B = ntt(ntt_data, B, False)
    for i in range(len(A)):
        A[i] = A[i] * B[i] % ntt_data[0]
    return ntt(ntt_data, A, True)[:LA + LB - 1]


ntt_precompute(ntt_data)

@njit((i8[:], ), cache=True)
def precompute(A):
    exp = np.zeros(MOD - 1, np.int64)
    log = np.zeros(MOD, np.int64)
    exp[0] = 1
    for k in range(1, MOD - 1):
        exp[k] = exp[k - 1] * R % MOD
        log[exp[k]] = k
    
    B = np.zeros(MOD - 1, np.int64)
    for x in A:
        if x != 0:
            B[log[x]] += 1
    return exp, log, B

@njit((i8[:], i8[:], i8[:], i8[:]), cache=True)
def get_ans(A, C, exp, log):
    
    ans = 0
    for i in range(len(C)):
        x = exp[i % (MOD - 1)]
        ans += x * C[i]
    for a in A:
        x = a * a % MOD
        ans -= x
    return ans // 2

A = np.array(read().split(), np.int64)[1:]

exp, log, B = precompute(A)

Bh, Bl = B >> 9, B & 511
x = ntt_convolve(ntt_data, Bh, Bh) << 18
x += ntt_convolve(ntt_data, Bh, Bl) << 10
x += ntt_convolve(ntt_data, Bl, Bl)

print(get_ans(A, x, exp, log))

Solution (Python):
import sys
import numpy as np
import numba
from numba import njit, b1, i4, i8, f8

read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

MOD = 200_003
R = 2  

NTTMOD = 998_244_353
NTT_root = 3
max_ntt = 20


@njit((i8, i8, i8), cache=True)
def mpow(a, n, MOD):
    p = 1
    while n:
        if n & 1:
            p = p * a % MOD
        a = a * a % MOD
        n >>= 1
    return p


@njit((i8, i8), cache=True)
def ntt_precompute(N, primitive_root):
    roots = np.empty(1 << N, np.int64)
    iroots = np.empty(1 << N, np.int64)
    roots[0] = iroots[0] = 1
    for n in range(N):
        x = mpow(primitive_root, (NTTMOD - 1) >> n + 2, NTTMOD)
        y = mpow(x, NTTMOD - 2, NTTMOD)
        for i in range(1 << n):
            roots[(1 << n) | i] = roots[i] * x % NTTMOD
            iroots[(1 << n) | i] = iroots[i] * y % NTTMOD
    return roots, iroots


@njit((i8[:], i8[:], i8[:], numba.b1), cache=True)
def ntt(roots, iroots, A, inverse):
    N = len(A)
    if not inverse:
        m = N >> 1
        while m:
            for k in range(N // m >> 1):
                s = 2 * m * k
                for i in range(s, s + m):
                    x, y = A[i], A[i + m] * roots[k]
                    A[i], A[i + m] = x + y, x - y
            A %= NTTMOD
            m >>= 1
    else:
        m = 1
        while m < N:
            for k in range(N // m >> 1):
                s = 2 * m * k
                for i in range(s, s + m):
                    x, y = A[i], A[i + m]
                    A[i], A[i + m] = x + y, (x - y) * iroots[k]
            A %= NTTMOD
            m <<= 1
        for i in range(N):
            A[i] = A[i] * mpow(N, NTTMOD - 2, NTTMOD) % NTTMOD
    return A


@njit((i8[:], i8[:], i8[:], i8[:]), cache=True)
def ntt_convolve(roots, iroots, A, B):
    A = ntt(roots, iroots, A.copy(), False)
    B = ntt(roots, iroots, B.copy(), False)
    for i in range(len(A)):
        A[i] = A[i] * B[i] % NTTMOD
    ntt(roots, iroots, A, True)
    return A


roots, iroots = ntt_precompute(max_ntt, NTT_root)

@njit((i8[:], ), cache=True)
def precompute(A):
    exp = np.zeros(MOD - 1, np.int64)
    log = np.zeros(MOD, np.int64)
    exp[0] = 1
    for k in range(1, MOD - 1):
        exp[k] = exp[k - 1] * R % MOD
        log[exp[k]] = k
    
    B = np.zeros(MOD - 1, np.int64)
    for x in A:
        if x != 0:
            B[log[x]] += 1
    return exp, log, B

@njit((i8[:], i8[:], i8[:], i8[:]), cache=True)
def get_ans(A, C, exp, log):
    
    ans = 0
    for i in range(len(C)):
        x = exp[i % (MOD - 1)]
        ans += x * C[i]
    for a in A:
        x = a * a % MOD
        ans -= x
    return ans // 2

A = np.array(read().split(), np.int64)[1:]

exp, log, B = precompute(A)

B = np.append(B, [0] * ((1 << 19) - len(B)))
Bh, Bl = B >> 9, B & 511
x = ntt_convolve(roots, iroots, Bh, Bh) << 18
x += ntt_convolve(roots, iroots, Bh, Bl) << 10
x += ntt_convolve(roots, iroots, Bl, Bl)

print(get_ans(A, x, exp, log))

Solution (Python):
import numpy as np
import scipy.fftpack as fft

mod = 200003
n = int(input())
a = []
for i in map(int, input().split()):
    if i % mod == 0:
        n -= 1
    else:
        a.append(i % mod)

tab = [0] * mod
j = 1
for i in range(mod-1):
    tab[j] = i
    j = (j*2) % mod

b = [0] * (mod-1)
for i in a:
    b[tab[i]] += 1

b = np.array(b, dtype=np.float64)
b = fft.fft(b)
b = b*b
b = fft.ifft(b)

ans = 0
for i, j in enumerate(b):
    ans += int(round(j.real))*pow(2, i, mod)

for i in a:
    ans -= i*i % mod
print(ans//2)

Solution (Python):
import sys
INF = 1 << 60
MOD = 10**9 + 7 
sys.setrecursionlimit(2147483647)
input = lambda:sys.stdin.readline().rstrip()
from cmath import pi, rect

def _fft(A, inverse = False):
    step = N
    for k in range(logN):
        step >>= 1
        w = rect(1, pi / (1 << k) * (1 - 2 * inverse))
        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, A[ps] - A[pt] * wj
            wj *= w
        A = nA
    return A

def convolution(f, g):
    Ff, Fg = _fft(f + [0] * (N - len(f))), _fft(g + [0] * (N - len(g)))
    fg = _fft([a * b / N for a, b in zip(Ff, Fg)], inverse = True)
    del fg[len(f) + len(g) - 1:]
    return fg

P = 200_003
logN = (2 * P - 1).bit_length()
N = 1 << logN
r = 2
exp = [1] * (P - 1)
log = [0] * P

e = 1
for i in range(1, P - 1):
    e *= r
    e %= P
    exp[i] = e
    log[e] = i

def resolve():
    n = int(input())
    A = list(map(int, input().split()))
    B = [0] * P
    for a in A:
        if a == 0:
            continue
        B[log[a]] += 1
    C = [0] * P
    for i, b in enumerate(convolution(B, B)):
        C[exp[i % (P - 1)]] += round(abs(b))
    for a in A:
        C[(a * a) % P] -= 1
    ans = 0
    for i in range(1, P):
        ans += i * C[i] // 2
    print(ans)
resolve()

Solution (Python):
import sys
def I(): return int(sys.stdin.readline().rstrip())
def LI(): return list(map(int,sys.stdin.readline().rstrip().split()))  

N = I()
A = LI()
B = []
for i in range(N):
    if A[i] != 0:
        B.append(A[i])
N = len(B)
P = 200003



C = []  
d = {}  
a = 1
for i in range(P-1):
    C.append(a)
    d[a] = i
    a *= 2
    a %= P

D = [0]*(P-1)
for i in range(N):
    D[d[B[i]]] += 1

n = 2**((2*P-3).bit_length())
D += [0]*(n-(P-1))


from cmath import pi,exp

def dfs(f,start,step,N,inverse=False):
    if N == 2:
        a,b = f[start],f[start+step]
        return [a+b,a-b]
    if inverse is False:
        f0 = dfs(f,start,2*step,N//2,inverse=False)
        f1 = dfs(f,start+step,2*step,N//2,inverse=False)
        z = exp(2j*pi/N)
    else:
        f0 = dfs(f,start,2*step,N//2,inverse=True)
        f1 = dfs(f,start+step,2*step,N//2,inverse=True)
        z = exp(-2j*pi/N)
    z0 = 1.0
    for i in range(N//2):
        a,b = f0[i],f1[i]*z0
        f0[i],f1[i] = a+b,a-b
        z0 *= z
    return f0+f1

def DFT(f,N):  
    if N == 1:
        return f
    return dfs(f,0,1,N)

def iDFT(F,N):  
    if N == 1:
        return F
    f = dfs(F,0,1,N,inverse=True)
    for i in range(N):
        f[i] /= N
    return f



def product(f,g,N):  
    F = DFT(f,N)
    G = DFT(g,N)
    FG = [F[i]*G[i] for i in range(N)]
    fg = iDFT(FG,N)
    for i in range(N):
        fg[i] = int(fg[i].real+0.5)
    return fg


E = product(D,D,n)
ans = 0
for i in range(P-1):
    ans += (C[i]*(E[i]+E[i+(P-1)]))

for i in range(N):
    ans -= (B[i]**2) % P

print(ans//2)

Solution (Python):
import numpy as np
from scipy import signal
import pprint

n = input()
a = np.array(list(map(int, input().split())), dtype=np.int64)

p = 200003
g = 2

c = np.zeros(p, dtype=np.int64)
d = np.zeros(p, dtype=np.int64)

k = 1
for i in range(p-1):
  c[k] = i
  d[i] = k
  k = k*g%p

b = np.zeros(p, dtype=np.float64)

for x in a:
  if x > 0:
    w = int(c[x])
    b[w] += 1

res = signal.fftconvolve(b, b, mode='full')


ans = 0

for i in range(len(res)):
  if res[i] > 0.5:
    w = int(round(res[i]))
    ans += w * d[i % (p-1)]

for x in a:
  ans -= x*x%p

ans = ans//2

print(ans)

Solution (Python):




import numpy as np
import sys

read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines
sys.setrecursionlimit(10 ** 7)

mod = 200003
n = int(readline())
A = list(map(int, input().split()))
F = np.zeros(mod, dtype=int)
prime = 2
now = 1
pow2 = [0] * (mod)
log2 = [-1] * (mod)

for i in range(mod):
    pow2[i] = now
    log2[pow2[i]] = i
    now *= prime
    now %= mod




for i in range(n):
    if A[i]:
        
        pos = log2[A[i]]
        F[pos] += 1



def fft_convolve(f, g):
    """
    数列 (多項式) f, g の畳み込みの計算．上下 15 bitずつ分けて計算することで，
    30 bit以下の整数，長さ 250000 程度の数列での計算が正確に行える．
    """
    fft = np.fft.rfft
    ifft = np.fft.irfft
    Lf = len(f)
    Lg = len(g)
    L = Lf + Lg - 1
    fft_len = 1 << L.bit_length()
    fl = f & (1 << 15) - 1
    fh = f >> 15
    gl = g & (1 << 15) - 1
    gh = g >> 15
    conv = lambda f, g: ifft(fft(f, fft_len) * fft(g, fft_len))[:L]
    x = conv(fl, gl)
    y = conv(fl + fh, gl + gh)
    z = conv(fh, gh)
    a, b, c = map(lambda x: (x + 0.5).astype(np.int64), [x, y, z])
    return a + ((b - a - c) << 15) + (c << 30)


x = fft_convolve(F, F)
ans1 = 0
L = len(x)
for i in range(L):
    ans1 += pow2[i % (mod - 1)] * x[i]  
for i in range(n):
    ans1 -= A[i] * A[i] % mod
print(ans1 // 2)

Solution (Python):
import sys
import numpy as np

def fft_convolve(f, g):
    """
    数列 (多項式) f, g の畳み込みの計算．上下 15 bitずつ分けて計算することで，
    30 bit以下の整数，長さ 250000 程度の数列での計算が正確に行える．
    """
    fft = np.fft.rfft
    ifft = np.fft.irfft
    Lf = len(f)
    Lg = len(g)
    L = Lf + Lg - 1
    fft_len = 1 << L.bit_length()
    fl = f & (1 << 15) - 1
    fh = f >> 15
    gl = g & (1 << 15) - 1
    gh = g >> 15

    def conv(f, g):
        return ifft(fft(f, fft_len) * fft(g, fft_len))[:L]

    x = conv(fl, gl)
    y = conv(fl + fh, gl + gh)
    z = conv(fh, gh)
    a, b, c = map(lambda x: (x + .5).astype(np.int64), [x, y, z])
    return (a + ((b - a - c) << 15) + (c << 30))


ROOT=2
P=200003
pw=[1]*(P-1)
lg=[-1]*P

for i in range(1,len(pw)):
    pw[i]=ROOT*pw[i-1]%P
for i in range(len(pw)):
    lg[pw[i]]=i

N=int(input())
A=np.array([int(x) for x in input().split()],dtype=np.int64)
B=[0]*P
for a in A:
    if a==0:
        continue
    B[lg[a]]+=1
B=np.array(B,dtype=np.int64)
C=fft_convolve(B,B)
ans=0
for a in A:
    ans-=a*a%P
for i in range(len(C)):
    ans+=C[i]*pw[i%(P-1)]
ans//=2
print(ans)

Solution (Python):

import numpy as np
M=int(input())
P=200003
a=[int(i) for i in input().split()]
b=[0 for i in range(P)]
for i in range(M):
    b[a[i]]+=1
c=[0 for i in range(P-1)]
tmp=1
for i in range(P-1):
    c[i]=b[tmp]
    tmp=tmp*2
    tmp%=P
N = P-1
n0 = 2**int(np.ceil(np.log2(2*N - 1)))

A = np.zeros(n0)
B = np.zeros(n0)
for i in range(N):
    A[i],B[i]=c[i],c[i]

C = np.fft.ifft( np.fft.fft(A)*np.fft.fft(B) )

ret = [int(ci) for ci in np.real(C[:2*N - 1] + 0.5)]
ans=0
for k in range(2*N-1):
    ans+=ret[k]*pow(2,k,P)
for i in range(M):
    ans-=(a[i]*a[i])%P
print(ans//2)

Solution (Python):
def convolve(A, B):
    
    
    import numpy as np
    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
import numpy as np
input = sys.stdin.readline

P = 200003
N = int(input())
A = sorted(map(int, input().split()))
A = [a for a in A if a]
Bucket = [0] * P
for a in A:
    Bucket[a] += 1

Perm = []
Perm_inv = [0] * (P)
p = 1
for i in range(P):
    Perm.append(p)
    p = p*2 % P
for i, p in enumerate(Perm):
    Perm_inv[p] = i




def perm(arr, P):
    res = [0] * len(arr)
    for a, p in zip(arr, P):
        res[p] = a
    return res

B = perm(Bucket, Perm_inv)
C = convolve(B, B)
C_ = C[:P].copy()
C_[:len(C)-P+1] += C[P-1:]
C_ = perm(C_, Perm)



Ans = C_ * np.arange(P)
an = Ans.sum()
diag = sum(a*a%P for a in A)

ans = (an-diag) // 2
print(ans)

Solution (Python):
import sys
import numpy as np

def fft_convolve(f, g):
    """
    数列 (多項式) f, g の畳み込みの計算．上下 15 bitずつ分けて計算することで，
    30 bit以下の整数，長さ 250000 程度の数列での計算が正確に行える．
    """
    fft = np.fft.rfft
    ifft = np.fft.irfft
    Lf = len(f)
    Lg = len(g)
    L = Lf + Lg - 1
    fft_len = 1 << L.bit_length()
    fl = f & (1 << 15) - 1
    fh = f >> 15
    gl = g & (1 << 15) - 1
    gh = g >> 15

    def conv(f, g):
        return ifft(fft(f, fft_len) * fft(g, fft_len))[:L]

    x = conv(fl, gl)
    y = conv(fl + fh, gl + gh)
    z = conv(fh, gh)
    a, b, c = map(lambda x: (x + .5).astype(np.int64), [x, y, z])
    return (a + ((b - a - c) << 15) + (c << 30))


ROOT=2
P=200003
pw=[1]*(P-1)
lg=[-1]*P

for i in range(1,len(pw)):
    pw[i]=ROOT*pw[i-1]%P
for i in range(len(pw)):
    lg[pw[i]]=i

N=int(input())
A=np.array([int(x) for x in input().split()],dtype=np.int64)
B=[0]*P
for a in A:
    if a==0:
        continue
    B[lg[a]]+=1
B=np.array(B,dtype=np.int64)
C=fft_convolve(B,B)
for i in range(P-1,len(C)):
    C[i%(P-1)]+=C[i]
C=C[:P-1]

ans=0
for i in range(len(B)):
    C[2*i%(P-1)]-=B[i]*B[i]
    ans+=B[i]*(B[i]-1)//2*pw[2*i%(P-1)]
for i in range(len(C)):
    C[i]//=2

for i in range(len(C)):
    ans+=pw[i]*C[i]

print(ans)

Solution (Python):
k = 36
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(X)*2-1)]]

P = 200003
g = 2
D = {}
E = []
s = 1
for i in range(P - 1):
    E.append(s)
    D[s] = i
    s = s * g % P

N = int(input())
A = [int(a) for a in input().split()]
X = [0] * (P - 1)
ans = 0
for a in A:
    if a:
        ans -= a ** 2 % P
        X[D[a]] += 1

XX = li(st(nu(X) ** 2))
for i in range(P-2):
    XX[i] += XX[i+P-1]

for i, x in enumerate(XX[:P-1]):
    ans += E[i] * x
print(ans // 2)

Solution (Python):
import sys
import numpy as np
import numba
from numba import njit, b1, i4, i8, f8

read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

MOD = 200_003
R = 2  

@njit((i8[:], ), cache=True)
def precompute(A):
    exp = np.zeros(MOD - 1, np.int64)
    log = np.zeros(MOD, np.int64)
    exp[0] = 1
    for k in range(1, MOD - 1):
        exp[k] = exp[k - 1] * R % MOD
        log[exp[k]] = k
    
    B = np.zeros(MOD - 1, np.int64)
    for x in A:
        if x != 0:
            B[log[x]] += 1
    return exp, log, B

@njit((i8[:], i8[:], i8[:], i8[:]), cache=True)
def get_ans(A, C, exp, log):
    
    ans = 0
    for i in range(len(C)):
        x = exp[i % (MOD - 1)]
        ans += x * C[i]
    for a in A:
        x = a * a % MOD
        ans -= x
    return ans // 2

A = np.array(read().split(), np.int64)[1:]

exp, log, B = precompute(A)
fft, ifft = np.fft.fft, np.fft.ifft
fft_len = 1 << 20
FB = fft(B, fft_len)
C = np.rint(ifft(FB * FB, fft_len)).astype(np.int64)

print(get_ans(A, C, exp, log))