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 numpy as np

def convolve(a, b):
	n = 1
	n_a, n_b = len(a), len(b)
	while n < n_a + n_b - 1:
		n <<= 1
	A = np.fft.rfft(a, n)
	B = np.fft.rfft(b, n)

	C = A * B
	c = np.fft.irfft(C, n)
	c = np.rint(c).astype(np.int64)

	return c[:n_a + n_b - 1]

p = 200003
n = int(input())
b = list(map(int, input().split()))
a = [0 for _ in range(p-1)]
ref = [0 for _ in range(p)]
po = [0 for _ in range(p-1)]
y = 1
for i in range(p-1):
	ref[y] = i
	po[i] = y
	y = (y * 2) % p
ans = 0
for x in b:
	if x == 0:
		continue
	a[ref[x]] += 1
	ans -= (x*x) % p
conv = convolve(a, a)
for i in range(2*p-3):
	ans += po[i % (p-1)] * conv[i]
print(ans // 2)

Solution (Python):
import numpy as np
from numpy.fft import fft
from numpy.fft import ifft
class FFT:
  def exe(s, A, B):
    N, M = 1, len(A) + len(B) - 1
    while N < M: N <<= 1
    A, B = s.arrN(N, A), s.arrN(N, B)
    A = fft(A) * fft(B)
    A = ifft(A).real[:M] + 0.5
    return list(map(int, A))
  def arrN(s, N, L):
    return np.zeros(N) + (L + [0] * (N - len(L)))

N = int(input())
a = list(map(int, input().split()))
MOD = 200003

r = 2
EXP = [0] * MOD
LOG = [0] * MOD
EXP[0] = 1


for i in range(1, MOD):
  EXP[i] = EXP[i - 1] * r % MOD
  LOG[EXP[i]] = i

b = [0] * MOD
for n in a:
  
  if n == 0: continue 
  b[LOG[n]] += 1



f = FFT()
c = f.exe(b, b)

ans = 0
t = 1
for n in c:
  
  ans += n * t
  t = t * r % MOD


for i in a:
  ans -= i * i % MOD


ans //= 2

print(ans)

Solution (Python):
import sys; input = sys.stdin.buffer.readline
sys.setrecursionlimit(10**7)
from collections import defaultdict
mod = 200003; INF = float("inf")

import numpy as np

def fft_np(a, b, n):
	N = 2 ** ((n - 1).bit_length() + 1) - n
	a = np.append(a, [0] * N)
	b = np.append(b, [0] * N)
	A = np.fft.rfft(a)
	B = np.fft.rfft(b)
	C = A * B
	c = np.fft.irfft(C)
	c = np.rint(c).astype(int)

	return c

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


def main():
	N = int(input())
	A = getlist()
	B = np.array([0] * mod)
	conv = defaultdict(int)
	num = 2
	for i in range(1, mod):
		conv[num] = i
		num = (num * 2) % mod

	for i in range(N):
		if A[i] == 0:
			continue
		x = conv[A[i]]
		B[x] += 1

	

	C = fft_np(B, B, mod)
	
	M = len(C)
	ans = 0
	for i in range(M):
		ans += C[i] * pow(2, i, mod)

	for i in range(mod):
		ans -= B[i] * pow(2, 2 * i, mod)

	print(ans // 2)



if __name__ == '__main__':
	main()

Solution (Python):
import sys; input = sys.stdin.buffer.readline
sys.setrecursionlimit(10**7)
from collections import defaultdict
mod = 200003; INF = float("inf")

import numpy as np

def fft_np(a, b, n):
	N = 2 ** ((n - 1).bit_length() + 1) - n
	a = np.append(a, [0] * N)
	b = np.append(b, [0] * N)
	A = np.fft.fft(a)
	B = np.fft.fft(b)
	C = A * B
	c = np.fft.ifft(C)
	c = np.rint(c).astype(int)

	return c

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


def main():
	N = int(input())
	A = getlist()
	B = np.array([0] * mod)
	conv = defaultdict(int)
	num = 2
	for i in range(1, mod):
		conv[num] = i
		num = (num * 2) % mod

	for i in range(N):
		if A[i] == 0:
			continue
		x = conv[A[i]]
		B[x] += 1

	

	C = fft_np(B, B, mod)
	
	M = len(C)
	ans = 0
	for i in range(M):
		ans += C[i] * pow(2, i, mod)

	for i in range(mod):
		ans -= B[i] * pow(2, 2 * i, mod)

	print(ans // 2)



if __name__ == '__main__':
	main()

Solution (Python):
import scipy.fftpack as f;input();m=200003;t=[0]*m;a=[0]*2*m;s=0
for i in range(m):t[pow(2,i,m)]=i
for i in map(int,input().split()):
	if i:a[t[i]]+=1;s-=i*i%m
for i,j in enumerate(f.ifft(f.fft(a)**2)):s+=int(round(j.real))*pow(2,i,m)
print(s//2)

Solution (Python):
import scipy.fftpack as f;input();m=200003;t=[0]*m;a=[0]*2*m;s=0
for i in range(m):t[pow(2,i,m)]=i
for i in map(int,input().split()):
	if i:a[t[i]]+=1;s-=i*i%m
for i,j in enumerate(f.ifft(f.fft(a)**2)):s+=int(round(j.real))*pow(2,i,m)
print(s//2)

Solution (Python):
import scipy.fftpack as f;input();m=200003;t=[0]*m;a=[0.0]*2*m;s=0
for i in range(m):t[pow(2,i,m)]=i
for i in map(int,input().split()):
	if i:a[t[i]]+=1;s-=i*i%m
for i,j in enumerate(f.ifft(f.fft(a)**2)):s+=int(round(j.real))*pow(2,i,m)
print(s//2)

Solution (Python):
k = 40 
MOD = 200003
list_size = 2 * MOD - 1


def n2l(num, list_size = list_size):
    "number to list"
    s = '{0:0{1}b}'.format(num, k * list_size)
    return [int(s[i-k:i], 2) for i in range(k * list_size, 0, -k)]

def l2n(ls):
    "list to number"
    return int(''.join(('{0:0{1}b}'.format(n, k) for n in reversed(ls))), 2)


def rn2l(num, list_size = list_size):
    "reversed n2l"
    s = '{0:0{1}b}'.format(num, k * list_size)
    return [int(s[i:i+k], 2) for i in range(0, k * list_size, k)]

def rl2n(ls):
    "reversed l2n"
    return int(''.join(('{0:0{1}b}'.format(n, k) for n in ls)), 2)


N = int(input())
As = sorted(map(int, input().split()))

two_pow = [1] * MOD
for i in range(1, MOD):
    two_pow[i] = two_pow[i-1] * 2 % MOD

two_pow_inv = [0] * MOD
for i, e in enumerate(two_pow[:-1]):
    two_pow_inv[e] = i

cnt = [0] * MOD 
for A in As:
    if A != 0:
        cnt[two_pow_inv[A]] += 1

cnt_conv = rn2l(rl2n(cnt) ** 2)
ans = 0
for i, cnt in enumerate(cnt_conv):
    if i < MOD:
        ans += two_pow[i] * cnt
    else:
        ans += two_pow[i - MOD + 1] * cnt

ans -= sum((A**2 % MOD for A in As))
ans //= 2
print(ans)

Solution (Python):
k = 40 
MOD = 200003
list_size = 2 * MOD - 1


def n2l(num, list_size = list_size):
    "number to list"
    s = '{0:0{1}b}'.format(num, k * list_size)
    return [int(s[i-k:i], 2) for i in range(k * list_size, 0, -k)]

def l2n(ls):
    "list to number"
    return int(''.join(('{0:0{1}b}'.format(n, k) for n in reversed(ls))), 2)


def rn2l(num, list_size = list_size):
    "reversed n2l"
    s = '{0:0{1}b}'.format(num, k * list_size)
    return [int(s[i:i+k], 2) for i in range(0, k * list_size, k)]

def rl2n(ls):
    "reversed l2n"
    return int(''.join(('{0:0{1}b}'.format(n, k) for n in ls)), 2)


N = int(input())
As = sorted(map(int, input().split()))

two_pow = [1] * MOD
for i in range(1, MOD):
    two_pow[i] = two_pow[i-1] * 2 % MOD

two_pow_inv = [0] * MOD
for i, e in enumerate(two_pow[:-1]):
    two_pow_inv[e] = i

cnt = [0] * MOD 
for A in As:
    if A != 0:
        cnt[two_pow_inv[A]] += 1

cnt_conv = n2l(l2n(cnt) ** 2)
ans = 0
for i, cnt in enumerate(cnt_conv):
    if i < MOD:
        ans += two_pow[i] * cnt
    else:
        ans += two_pow[i - MOD + 1] * cnt

ans -= sum((A**2 % MOD for A in As))
ans //= 2
print(ans)

Solution (Python):
import numpy as np
N=int(input())
P=200003
a=[int(i) for i in input().split()]
b=[0 for i in range(P)]
for i in range(N):
    b[a[i]]+=1
M=1<<19
c=np.zeros(M)
for i in range(P-1):
    c[i]=b[pow(2,i,P)]
X=[int(i+0.1) for i in np.real(np.fft.ifft(np.fft.fft(c)*np.fft.fft(c)))]
ans=0
for k in range(M):
    ans+=X[k]*pow(2,k,P)
for i in range(N):
    ans-=(a[i]*a[i])%P
print(ans//2)

Solution (Python):
import numpy as np
from numpy.fft import rfft,irfft

p = 200003
N = int(input())
A = list(map(int,input().split()))
a = np.zeros(p-1)

g = 1
d = {}
p2 = [1]*(p-1)
for i in range(p-1):
    d[g] = i
    p2[i] = g
    g *= 2
    g %= p

for i in range(N):
    if A[i]==0:
        continue
    a[d[A[i]]] += 1

b = rfft(a,2*p-2)
b = irfft(b*b,2*p-2)
b = np.rint(b).astype(np.int64)

ans = 0
for i in range(p-1):
    ans += (b[i]+b[i+p-1])*p2[i]

for i in range(N):
    ans -= A[i]**2%p
ans //= 2
print(ans)

Solution (Python):

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)
    fft_f = [int(np.rint(val)) for val in fft_f]
    sum = 0
    for i in range(L):
        k = i % (M-1)
        sum += fft_f[i] * arr2val[k]

    sub = 0
    for a in A:
        sub += (a * a) % M
    print(int((sum - sub) // 2))



if __name__ == "__main__":
    run()

Solution (Python):

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)
    fft_f = [round(val.real) for val in fft_f]
    sum = 0
    for i in range(L):
        k = i % (M-1)
        sum += fft_f[i] * arr2val[k]

    sub = 0
    for a in A:
        sub += (a * a) % M
    print(int((sum - sub) // 2))



if __name__ == "__main__":
    run()

Solution (Python):

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()

Solution (Python):
import sys
import numpy as np

sys.setrecursionlimit(10 ** 8)
ini = lambda: int(sys.stdin.readline())
inl = lambda: [int(x) for x in sys.stdin.readline().split()]
ins = lambda: sys.stdin.readline().rstrip()
debug = lambda *a, **kw: print("\033[33m", *a, "\033[0m", **dict(file=sys.stderr, **kw))

P = 200003
PROOT = 2


def calc_tables():
    prev = [None] * P
    val = 1
    for i in range(P):
        prev[val] = i
        val = (val * PROOT) % P
    return prev


def convolve(f, g):
    fft_len = 1
    alen = len(f) + len(g) - 1
    while fft_len < alen:
        fft_len *= 2

    fa = np.fft.rfft(f, fft_len)
    fb = np.fft.rfft(g, fft_len)

    h = np.fft.irfft(fa * fb)
    h = np.rint(h).astype(np.int64)
    return h


def solve():
    prev = calc_tables()

    _ = ini()
    a = [x for x in inl() if x != 0]
    b = [prev[x] for x in a]
    a = np.array(a, dtype=np.int64)
    m = len(a)

    b_counts = np.zeros(P, np.int64)
    for x in b:
        b_counts[x] += 1

    root_pow = np.ones(2 * P, dtype=np.int64)
    for i in range(1, 2 * P):
        root_pow[i] = root_pow[i - 1] * PROOT % P

    c = convolve(b_counts, b_counts)
    ans = 0
    for k in range(2 * P - 1):
        ans += c[k] * root_pow[k]
    for k in range(P):
        ans -= b_counts[k] * root_pow[2 * k]
    return ans // 2


print(solve())

Solution (Python):
import numpy as np


class FastFourierTransform:
	"""FFT
	"""

	def multiply(self, F,G):
		""" multiply by FFT
		"""
		N=1
		
		while N < len(F) + len(G):
			N *= 2
		
		nf = np.zeros(N,int)
		nf[:len(F)]=F
		ng = np.zeros(N,int)
		ng[:len(G)]=G
		nf=np.fft.fft(nf)
		ng=np.fft.fft(ng)
		
			
		nf *= ng

		
		nf=np.fft.ifft(nf)

		ans = np.zeros(N,int)
		for n in range(N):
			
			
			ans[n] = int(nf[n].real + 0.5)

		return ans

if __name__ == "__main__":
	N=int(input())
	A=list(map(int,input().split()))

	P=200_003
	rootP=2

	
	order=[1]
	x=rootP
	while x!=1:
		order.append(x)
		x *= rootP
		x %= P

	ORDER_SIZE=len(order)	

	
	where=np.zeros(P,int)
	for i in range(ORDER_SIZE):
		where[ order[i] ] = i

	
	prevCnt = np.zeros(P,int)
	resCnt = np.zeros(P,int)
	
	
	for n in range(N):
		if A[n] != 0:
			prevCnt[ where[A[n]] ] += 1
			
			
			resCnt[ (A[n] * A[n]) % P ] -= 1

	
	clsFFT = FastFourierTransform()
	fftRes = clsFFT.multiply(prevCnt,prevCnt)

	for r in range(len(fftRes)):
		if fftRes[r] != 0:
			
			resCnt[ order[ r % ORDER_SIZE ] ] += fftRes[r]

	ans = 0
	for p in range(P):
		if resCnt[p] != 0:
			
			
			
			
			ans += (resCnt[p] // 2) * p

	print(ans)

Solution (Python):
import numpy as np

if __name__ == "__main__":
	N=int(input())
	A=list(map(int,input().split()))

	P=200_003
	rootP=2

	order=[1]
	x=rootP
	while x!=1:
		order.append(x)
		x *= rootP
		x %= P

	ORDER_SIZE=len(order)	

	where=np.zeros(P,int)
	for i in range(ORDER_SIZE):
		where[ order[i] ] = i

	prevCnt = np.zeros(P,int)
	resCnt = np.zeros(P,int)
	for n in range(N):
		if A[n] != 0:
			prevCnt[ where[A[n]] ] += 1
			resCnt[ (A[n] * A[n]) % P ] -= 1

	N_=1
	while N_ < 2*len(prevCnt): N_ *= 2
	nf = np.zeros(N_,int)
	nf[:len(prevCnt)]=prevCnt
	nf=np.fft.rfft(nf)
	fftRes = np.rint( np.fft.irfft(nf*nf) )

	for r in range(len(fftRes)):
		if fftRes[r] != 0:
			resCnt[ order[ r % ORDER_SIZE ] ] += fftRes[r]

	ans = 0
	for p in range(P):
		if resCnt[p] != 0:
			ans += (resCnt[p] // 2) * p

	print(ans)

Solution (Python):
import numpy as np


class FastFourierTransform:
	"""FFT
	"""

	def multiplySameArr(self,A):
		""" multiply by FFT(A*A)
		"""
		N=1
		
		while N < 2*len(A): N *= 2

		
		nf = np.zeros(N,int)
		nf[:len(A)]=A
		
		nf=np.fft.rfft(nf)
		nf = nf * nf

		
		
		nf=np.fft.irfft(nf)

		
		ans = np.rint( nf )
		
		
		
		
		
		

		return ans

if __name__ == "__main__":
	N=int(input())
	A=list(map(int,input().split()))

	P=200_003
	rootP=2

	
	order=[1]
	x=rootP
	while x!=1:
		order.append(x)
		x *= rootP
		x %= P

	ORDER_SIZE=len(order)	

	
	where=np.zeros(P,int)
	for i in range(ORDER_SIZE):
		where[ order[i] ] = i

	
	prevCnt = np.zeros(P,int)
	resCnt = np.zeros(P,int)
	for n in range(N):
		if A[n] != 0:
			prevCnt[ where[A[n]] ] += 1
			
			resCnt[ (A[n] * A[n]) % P ] -= 1

	
	clsFFT = FastFourierTransform()
	fftRes = clsFFT.multiplySameArr(prevCnt)

	for r in range(len(fftRes)):
		if fftRes[r] != 0:
			
			resCnt[ order[ r % ORDER_SIZE ] ] += fftRes[r]

	ans = 0
	for p in range(P):
		if resCnt[p] != 0:
			
			
			
			ans += (resCnt[p] // 2) * p

	print(ans)

Solution (Python):
import numpy as np


class FastFourierTransform:
	"""FFT
	"""

	def multiplySameArr(self,A):
		""" multiply by FFT(A*A)
		"""
		N=1
		
		while N < 2*len(A): N *= 2

		
		nf = np.zeros(N,int)
		nf[:len(A)]=A
		nf=np.fft.fft(nf)
		nf = nf * nf

		
		nf=np.fft.ifft(nf)

		
		ans = np.zeros(N,int)
		for n in range(N):
			
			
			if nf[n].real != 0:
				ans[n] = int(nf[n].real + 0.5)

		return ans

if __name__ == "__main__":
	N=int(input())
	A=list(map(int,input().split()))

	P=200_003
	rootP=2

	
	order=[1]
	x=rootP
	while x!=1:
		order.append(x)
		x *= rootP
		x %= P

	ORDER_SIZE=len(order)	

	
	where=np.zeros(P,int)
	for i in range(ORDER_SIZE):
		where[ order[i] ] = i

	
	prevCnt = np.zeros(P,int)
	resCnt = np.zeros(P,int)
	for n in range(N):
		if A[n] != 0:
			prevCnt[ where[A[n]] ] += 1
			
			resCnt[ (A[n] * A[n]) % P ] -= 1

	
	clsFFT = FastFourierTransform()
	fftRes = clsFFT.multiplySameArr(prevCnt)

	for r in range(len(fftRes)):
		if fftRes[r] != 0:
			
			resCnt[ order[ r % ORDER_SIZE ] ] += fftRes[r]

	ans = 0
	for p in range(P):
		if resCnt[p] != 0:
			
			
			
			ans += (resCnt[p] // 2) * p

	print(ans)

Solution (Python):
import sys
sys.setrecursionlimit(10**7) 
import math
from copy import copy, deepcopy
from copy import deepcopy as dcp
from operator import itemgetter
from bisect import bisect_left, bisect, bisect_right

from collections import deque, defaultdict


from collections import Counter

from itertools import accumulate,combinations,permutations,product

from heapq import heapify,heappop,heappush

from functools import reduce,lru_cache

from decimal import Decimal

import numpy as np
def input(): 
    x=sys.stdin.readline()
    return x[:-1] if x[-1]=="\n" else x
def printe(*x):print("
def printl(li): _=print(*li, sep="\n") if li else None
def argsort(s, return_sorted=False): 
    inds=sorted(range(len(s)), key=lambda k: s[k])
    if return_sorted: return inds, [s[i] for i in inds]
    return inds
def alp2num(c,cap=False): return ord(c)-97 if not cap else ord(c)-65
def num2alp(i,cap=False): return chr(i+97) if not cap else chr(i+65)
def matmat(A,B):
    K,N,M=len(B),len(A),len(B[0])
    return [[sum([(A[i][k]*B[k][j]) for k in range(K)]) for j in range(M)] for i in range(N)]
def matvec(M,v):
    N,size=len(v),len(M)
    return [sum([M[i][j]*v[j] for j in range(N)]) for i in range(size)]
def T(M):
    n,m=len(M),len(M[0])
    return [[M[j][i] for j in range(n)] for i in range(m)]
def binr(x): return bin(x)[2:]
def bitcount(x): 
    x= x - ((x >> 1) & 0x5555555555555555)
    x= (x & 0x3333333333333333) + ((x >> 2) & 0x3333333333333333)
    x= (x + (x >> 4)) & 0x0f0f0f0f0f0f0f0f 
    x+= (x >> 8); x+= (x >> 16); x+= (x >> 32) 
    return x & 0x7f

def main():
    mod = 200003
    P=mod

    

    N = int(input())
    
    A = list(map(int, input().split())) 
    
    
    cur=1
    base=2
    atoi=[0]*P
    itoa=[0]*(P-1)
    cnt=[0]*(P-1)
    for i in range(0,P-1):
        atoi[cur]=i
        itoa[i]=cur
        cur=cur*base%mod

    ta2=0
    for a in A:
        if a==0:
            continue
        i=atoi[a]
        cnt[i]+=1
        ta2+=a*a%P
    rf=np.fft.rfft(cnt, 2*(P-1))
    X = np.rint(np.fft.irfft(rf * rf)).astype(np.int64)
    
    ans=0
    for i in range(P-1):
        ans+=itoa[i]*(X[i]+X[P-1+i])

    print(int((ans-ta2)//2))

if __name__ == "__main__":
    main()

Solution (Python):
import sys
sys.setrecursionlimit(10**7) 
import math
from copy import copy, deepcopy
from copy import deepcopy as dcp
from operator import itemgetter
from bisect import bisect_left, bisect, bisect_right

from collections import deque, defaultdict


from collections import Counter

from itertools import accumulate,combinations,permutations,product

from heapq import heapify,heappop,heappush

from functools import reduce,lru_cache

from decimal import Decimal
from numpy.fft import rfft, irfft

def input(): 
    x=sys.stdin.readline()
    return x[:-1] if x[-1]=="\n" else x
def printe(*x):print("
def printl(li): _=print(*li, sep="\n") if li else None
def argsort(s, return_sorted=False): 
    inds=sorted(range(len(s)), key=lambda k: s[k])
    if return_sorted: return inds, [s[i] for i in inds]
    return inds
def alp2num(c,cap=False): return ord(c)-97 if not cap else ord(c)-65
def num2alp(i,cap=False): return chr(i+97) if not cap else chr(i+65)
def matmat(A,B):
    K,N,M=len(B),len(A),len(B[0])
    return [[sum([(A[i][k]*B[k][j]) for k in range(K)]) for j in range(M)] for i in range(N)]
def matvec(M,v):
    N,size=len(v),len(M)
    return [sum([M[i][j]*v[j] for j in range(N)]) for i in range(size)]
def T(M):
    n,m=len(M),len(M[0])
    return [[M[j][i] for j in range(n)] for i in range(m)]
def binr(x): return bin(x)[2:]
def bitcount(x): 
    x= x - ((x >> 1) & 0x5555555555555555)
    x= (x & 0x3333333333333333) + ((x >> 2) & 0x3333333333333333)
    x= (x + (x >> 4)) & 0x0f0f0f0f0f0f0f0f 
    x+= (x >> 8); x+= (x >> 16); x+= (x >> 32) 
    return x & 0x7f

def main():
    mod = 200003
    P=mod

    

    N = int(input())
    
    A = list(map(int, input().split())) 
    
    
    cur=1
    base=2
    atoi=[0]*P
    itoa=[0]*(P-1)
    cnt=[0]*(P-1)
    for i in range(0,P-1):
        atoi[cur]=i
        itoa[i]=cur
        cur=cur*base%mod

    ta2=0
    for a in A:
        if a==0:
            continue
        i=atoi[a]
        cnt[i]+=1
        ta2+=a*a%P
    rf=rfft(cnt, 2*(P-1))
    X = irfft(rf * rf)
    
    ans=0
    for i in range(P-1):
        ans+=itoa[i]*round(X[i]+X[P-1+i])

    print(int((ans-ta2)//2))

if __name__ == "__main__":
    main()

Solution (Python):
import numpy 
 
def convolve(g, h):
  s=len(g)+len(h)+1
  n=1
  while n<s:
    n*=2
  gg=numpy.fft.fft(g+[0]*(n-len(g)))
  hh=numpy.fft.fft(h+[0]*(n-len(h)))
  ff=[gg[i]*hh[i] for i in range(n)]
  ff=numpy.fft.ifft(ff)
  ans=[int(ff[i].real+0.5) for i in range(n)]
  return ans

p=200003
n=int(input())
arr=list(map(int,input().split()))
pos={}
tmp=1
for i in range(p):
  pos[tmp]=i
  tmp*=2
  tmp%=p
cnt1=[0]*p
cnt2=[0]*p
for val in arr:
  if val!=0:
    cnt1[pos[val]]+=1
    cnt2[pos[val]]+=1
ans=convolve(cnt1,cnt2)
for val in arr:
  if val!=0:
    ans[pos[val]*2]-=1
total=0
for i in range(len(ans)):
  if ans[i]!=0:
    total+=(ans[i]//2)*pow(2,i,p)
print(total)

Solution (Python):
import numpy 
 
def convolve(g, h):
  s=len(g)+len(h)+1
  n=1
  while n<s:
    n*=2
  gg=numpy.fft.fft(g+[0]*(n-len(g)))
  hh=numpy.fft.fft(h+[0]*(n-len(h)))
  ff=[gg[i]*hh[i] for i in range(n)]
  ff=numpy.fft.ifft(ff)
  ans=[int(ff[i].real+0.5) for i in range(n)]
  return ans

p=200003
n=int(input())
arr=list(map(int,input().split()))
g=2
order=[1]
while g!=1:
  order.append(g)
  g*=2
  g%=p
pos=[0]*p
for i in range(len(order)):
  pos[order[i]]=i
cnt1=[0]*p
cnt2=[0]*p
coef=[0]*p
for val in arr:
  if val!=0:
    cnt1[pos[val]]+=1
    cnt2[pos[val]]+=1
    coef[(val*val)%p]-=1
ans=convolve(cnt1,cnt2)
for i in range(len(ans)):
  if ans[i]!=0:
    prod=order[i%(p-1)]
    coef[prod]+=ans[i]
total=0
for i in range(p):
  total+=(coef[i]//2)*i
print(total)

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 = 420_000
FB = fft(B, fft_len)
C = np.rint(ifft(FB * FB, fft_len)).astype(np.int64)

print(get_ans(A, C, 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  

@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 << 19
FB = fft(B, fft_len)
C = np.rint(ifft(FB * FB, fft_len)).astype(np.int64)

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

Solution (Python):
'''
自宅用PCでの解答
'''
import math

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


ipt = stdin.readline
setrecursionlimit(10**7)
mod = 200003
dir = [(-1,0),(0,-1),(1,0),(0,1)]
alp = "abcdefghijklmnopqrstuvwxyz"

def dft(f,n):
    if n == 1:
        return f
    f0 = dft(f[0::2],n//2)
    f1 = dft(f[1::2],n//2)
    zeta = complex(math.cos(2*math.pi/n),math.sin(2*math.pi/n))
    pow_zeta = 1
    for i in range(n):
        f[i] = f0[i%(n//2)] + pow_zeta*f1[i%(n//2)]
        pow_zeta *= zeta
    return f

def reverse_dft(f,n):
    if n == 1:
        return f
    f0 = reverse_dft(f[0::2],n//2)
    f1 = reverse_dft(f[1::2],n//2)
    zeta = complex(math.cos(2*math.pi/n),-math.sin(2*math.pi/n))
    pow_zeta = 1
    for i in range(n):
        f[i] = f0[i%(n//2)] + pow_zeta*f1[i%(n//2)]
        pow_zeta *= zeta
    return f

def multiply(g,h):
    n = 2**len(format((len(g)+len(h)+1),'b'))
    g += [0]*(n-len(g))
    h += [0]*(n-len(h))
    gg = dft(g,n)
    hh = dft(h,n)
    ff = [gi*hi for gi,hi in zip(gg,hh)]
    return [int(ffi.real/n+0.1) for ffi in reverse_dft(ff,n)]

def main():
    n = int(ipt())
    g = 2
    order = [1]
    for i in range(mod-2):
        order.append(order[-1]*g%mod)
    where = [0]*mod
    for i,oi in enumerate(order):
        where[oi] = i

    nms = [0]*(mod-1)
    for i in map(int,ipt().split()):
        if i == 0:
            continue
        nms[where[i]] += 1

    ans = 0
    for i,ni in enumerate(nms):
        ans -= ni*order[(i*2)%(mod-1)]

    nn = multiply(nms*1,nms*1)

    for i in range(2*(mod-1)-1):
        ans += nn[i]*order[i%(mod-1)]

    print(ans//2)

    return None

if __name__ == '__main__':
    main()

Solution (Python):
import scipy.fftpack as F
from numpy import *
m=200003
n=int(input())
a=[0]*m
b=[0]*(m-1)
ans=0
for i in range(m-1):
    a[pow(2,i,m)]=i
for i in map(int,input().split()):
    if i:
        b[a[i]]+=1
        ans-=i*i%m
b=array(b,dtype=float64)
b=F.fft(b)
b=b*b
b=F.ifft(b)
for i,j in enumerate(b):
    ans+=int(round(j.real))*pow(2,i,m)
print(ans//2)

Solution (Python):
import sys
import numpy as np

n, *aaa = map(int, sys.stdin.buffer.read().split())
aaa = [a for a in aaa if a != 0]

fft_len = 1 << 19

P = 200003
PR = 2
convert = [0] * P
rev_con = [0] * P
rev_con[0] = 1
tmp = 1
for i in range(1, P):
    tmp = tmp * PR % P
    convert[tmp] = i
    rev_con[i] = tmp

self_pair_products = 0
f = np.zeros(P)
for a in aaa:
    f[convert[a]] += 1
    self_pair_products += a * a % P

Ff = np.fft.rfft(f, fft_len)
Fff = Ff * Ff
g = np.rint(np.fft.irfft(Fff, fft_len))

g[:P - 1] += g[P - 1:2 * (P - 1)]
g[:fft_len - 2 * (P - 1)] += g[2 * (P - 1):]
g = g[:P - 1].astype(np.int64)
rev_con = np.array(rev_con[:-1])

ans = (g * rev_con).sum() - self_pair_products
print(ans // 2)

Solution (Python):
import sys
import numpy as np

n, *aaa = map(int, sys.stdin.buffer.read().split())
aaa = [a for a in aaa if a != 0]
n = len(aaa)

fft_len = 1 << 19

P = 200003
PR = 2
convert = [0] * P
rev_con = [0] * P
rev_con[0] = 1
tmp = 1
for i in range(1, P):
    tmp = tmp * PR % P
    convert[tmp] = i
    rev_con[i] = tmp

ans = 0
f = np.zeros(P)
for a in aaa:
    f[convert[a]] += 1
    ans -= a * a % P

Ff = np.fft.rfft(f, fft_len)
Fff = Ff * Ff
g = np.rint(np.fft.irfft(Fff, fft_len))
for i, c in enumerate(g):
    if c != 0:
        ans += rev_con[i % (P - 1)] * int(c + 1e-6)

print(ans // 2)

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