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



Sample input:
4
2019 0 2020 200002


Sample output:
474287

Submission status: Accepted

Solution code:
import sys
INF = 1 << 60
MOD = 10**9 + 7 # 998244353
sys.setrecursionlimit(2147483647)
input = lambda:sys.stdin.readline().rstrip()

prime = 998244353
root = 3
def _fmt(A, inverse = False):
    N = len(A)
    logN = (N - 1).bit_length()
    base = pow(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(f, g):
    N = 1 << (len(f) + len(g) - 2).bit_length()
    Ff, Fg = _fmt(f + [0] * (N - len(f))), _fmt(g + [0] * (N - len(g)))
    N_inv = pow(N, prime - 2, prime)
    fg = _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

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
    Bh = [b >> 9 for b in B]
    Bl = [b & 511 for b in B]
    for i, b in enumerate(convolution(Bh, Bh)):
        C[exp[i % (P - 1)]] += b << 18
    for i, b in enumerate(convolution(Bh, Bl)):
        C[exp[i % (P - 1)]] += b << 10
    for i, b in enumerate(convolution(Bl, Bl)):
        C[exp[i % (P - 1)]] += 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()

Submission status: Accepted

Solution code:
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  # primitive root

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
    # r^k の個数を B[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))

Submission status: Accepted

Solution code:
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  # primitive root

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
    # r^k の個数を B[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))

Submission status: Accepted

Solution code:
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)

Submission status: Accepted

Solution code:
import sys
INF = 1 << 60
MOD = 10**9 + 7 # 998244353
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()

Submission status: Accepted

Solution code:
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

# 2はPの原始根

C = []  # C[i] = (2**i mod P)
d = {}  # d[n] = (2**i == n たるi)
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):  # fを離散フーリエ変換
    if N == 1:
        return f
    return dfs(f,0,1,N)

def iDFT(F,N):  # Fを離散フーリエ逆変換
    if N == 1:
        return F
    f = dfs(F,0,1,N,inverse=True)
    for i in range(N):
        f[i] /= N
    return f

# deg(f)+deg(g)<Nとなる2冪Nを選び、f,gの長さがNになるように0を埋めておく

def product(f,g,N):  # 多項式f,gの積
    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)

Submission status: Accepted

Solution code:
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)

Submission status: Accepted

Solution code:
# フーリエ変換は、多項式の積動詞の積において指数の和が一定となるようなものをまとめる感じ
# Ci=Aj*Bk となるときはi=j+kのものをまとめている
# 今回はmod=200003でi+jが一致するものを集めたいので
# 係数→log→フーリエ変換→指数で戻す、という流れ？
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

# x**(log(Ai)) * x**(log(Aj)) とすれば
# convoluteすればx**(log(Ai)+log(Aj))となり、
# x**log(log(Ai * Aj))となってくれる
for i in range(n):
    if A[i]:
        # なのでlog(A[i])のところに1を立てる
        pos = log2[A[i]]
        F[pos] += 1


# maspyさんのyukicoder解答から
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]  # why mod-1????????
for i in range(n):
    ans1 -= A[i] * A[i] % mod
print(ans1 // 2)

# ans2 = 0
# for i in range(n):
#     for j in range(i + 1, n):
#         ans2 += A[i] * A[j] % mod
#         print((i, j), A[i] * A[j] % mod)
#         print(log2[A[i] * A[j] % mod])
# print(ans2)

Submission status: Accepted

Solution code:
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)

Submission status: Accepted

Solution code:
#https://qiita.com/sage-git/items/90c342fd2127ec39b874
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)

Submission status: Accepted

Solution code:
def convolve(A, B):
    # 畳み込み  # 要素は整数
    # 3 つ以上の場合は一度にやった方がいい
    import numpy as np
    dtype = np.int64  # np.float128 は windows では動かない？
    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()  # n 以上の最小の 2 冪
    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
#
# print(Perm[:9])
# print(Perm_inv[:9])

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)
# for i, c in enumerate(C_):
#     if c:
#         print(i%P)
Ans = C_ * np.arange(P)
an = Ans.sum()
diag = sum(a*a%P for a in A)

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

Submission status: Accepted

Solution code:
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)

Submission status: Accepted

Solution code:
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)

Submission status: Accepted

Solution code:
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  # primitive 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
    # r^k の個数を B[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))

Submission status: Runtime Error

Solution code:
from sys import stdin,stdout
def INPUT():return list(int(i) for i in stdin.readline().split())
def inp():return stdin.readline()
def out(x):return stdout.write(x)
import math
import random
J=int(1e9)+7
P=200003
###############################################################################\n=17
n=int(input())
A=INPUT()
sum=0
for i in ragne(n):
    for j in range(i+1,n):
        sum+=(A[i]*A[j])%P
print(sum)

Submission status: Runtime Error

Solution code:
import sys
sysread = sys.stdin.buffer.readline
read = sys.stdin.buffer.read
 
sys.setrecursionlimit(10**7)
import math
import numpy as np
 
INF = 1 << 50
 
 
	def run():
    N, *A = map(int, read().split())
    M = 200003
    P = [0] * (M-1)
    A = [a for a in A if a != 0]
    val2arr = {}
    arr2val = {}
    val = 1
    for p in range(M-1):
        val2arr[val] = p
        arr2val[p] = val
        val *= 2
        val %= M
 
 
    for a in A:
        k = val2arr[a]
        P[k] += 1
 
    L = 1 << (2 * M - 3).bit_length()
 
    fft_f = np.fft.ifft(np.fft.fft(P, L) ** 2)
 
    sum = 0
    for i in range(L):
        k = i % (M-1)
        sum += fft_f[i].real * arr2val[k]
 
    sub = 0
    for a in A:
        sub += (a * a) % M
    print(math.ceil((sum - sub) / 2))
 
 
 
	if __name__ == "__main__":
    run()

Submission status: Runtime Error

Solution code:
import sys
sysread = sys.stdin.buffer.readline
read = sys.stdin.buffer.read

sys.setrecursionlimit(10**7)
import math
import numpy as np
 
INF = 1 << 50
 
 
 def run():
    N, *A = map(int, read().split())
    M = 200003
    P = [0] * (M-1)
    A = [a for a in A if a != 0]
    val2arr = {}
    arr2val = {}
    val = 1
    for p in range(M-1):
        val2arr[val] = p
        arr2val[p] = val
        val *= 2
        val %= M
 
 
    for a in A:
        k = val2arr[a]
        P[k] += 1
 
    L = 1 << (2 * M - 3).bit_length()
 
    fft_f = np.fft.ifft(np.fft.fft(P, L) ** 2)
 
    sum = 0
    for i in range(L):
        k = i % (M-1)
        sum += fft_f[i].real * arr2val[k]
 
    sub = 0
    for a in A:
        sub += (a * a) % M
    print(math.ceil((sum - sub) / 2))
 
 
 
if __name__ == "__main__":
    run()