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; 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

	# print(B)

	C = fft_np(B, B, mod)
	# print(C)
	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()

Submission status: Accepted

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

	# print(B)

	C = fft_np(B, B, mod)
	# print(C)
	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()

Submission status: Accepted

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

Submission status: Accepted

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

Submission status: Accepted

Solution code:
k = 40 # 多倍長整数内で、一整数を何bitでもつか。
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 # cnt[i] counts 2 ** i % P
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)

Submission status: Accepted

Solution code:
k = 40 # 多倍長整数内で、一整数を何bitでもつか。
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 # cnt[i] counts 2 ** i % P
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)

Submission status: Accepted

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

Submission status: Accepted

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

Submission status: Accepted

Solution code:
# coding: utf-8
import sys
#from operator import itemgetter
sysread = sys.stdin.buffer.readline
read = sys.stdin.buffer.read
#from heapq import heappop, heappush
#from collections import defaultdict
sys.setrecursionlimit(10**7)
import math
#from itertools import product, accumulate, combinations, product
#import bisect
import numpy as np
#from copy import deepcopy
#from collections import deque
#from decimal import Decimal

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

Submission status: Accepted

Solution code:
# coding: utf-8
import sys
#from operator import itemgetter
sysread = sys.stdin.buffer.readline
read = sys.stdin.buffer.read
#from heapq import heappop, heappush
#from collections import defaultdict
sys.setrecursionlimit(10**7)
import math
#from itertools import product, accumulate, combinations, product
#import bisect
import numpy as np
#from copy import deepcopy
#from collections import deque
#from decimal import Decimal

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

Submission status: Accepted

Solution code:
# coding: utf-8
import sys
#from operator import itemgetter
sysread = sys.stdin.buffer.readline
read = sys.stdin.buffer.read
#from heapq import heappop, heappush
#from collections import defaultdict
sys.setrecursionlimit(10**7)
import math
#from itertools import product, accumulate, combinations, product
#import bisect
import numpy as np
#from copy import deepcopy
#from collections import deque
#from decimal import Decimal

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

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

Submission status: Accepted

Solution code:
import numpy as np


class FastFourierTransform:
	"""FFT
	"""

	def multiply(self, F,G):
		""" multiply by FFT
		"""
		N=1
		# 1.N > deg(g) + deg(h)となる2の冪乗を選ぶ
		while N < len(F) + len(G):
			N *= 2
		# 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)
		# for n in range(N):
			# nf[n] *= ng[n]
		nf *= ng

		# 3. 離散フーリエ逆変換
		nf=np.fft.ifft(nf)

		ans = np.zeros(N,int)
		for n in range(N):
			# 係数を1/N倍 => np.fftの場合、すでに実行済
			#ans.append( round(nf[n].real / 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 of 1 to P-1 after fft
	order=[1]
	x=rootP
	while x!=1:
		order.append(x)
		x *= rootP
		x %= P

	ORDER_SIZE=len(order)	# 200_002

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

	# count non zero items
	prevCnt = np.zeros(P,int)
	resCnt = np.zeros(P,int)
	# prevCnt2 = [ 0 for _ in range(P)]
	# resCnt = [ 0 for _ in range(P)]
	for n in range(N):
		if A[n] != 0:
			prevCnt[ where[A[n]] ] += 1
			# prevCnt2[ where[A[n]] ] += 1
			# remove Ax * Ax pattern count
			resCnt[ (A[n] * A[n]) % P ] -= 1

	# fft to cnt (calc nC2 and all Ax * Ay)
	clsFFT = FastFourierTransform()
	fftRes = clsFFT.multiply(prevCnt,prevCnt)

	for r in range(len(fftRes)):
		if fftRes[r] != 0:
			# get result(Ax*Ay==resCnt.index) Cnt by summary fftRes
			resCnt[ order[ r % ORDER_SIZE ] ] += fftRes[r]

	ans = 0
	for p in range(P):
		if resCnt[p] != 0:
			# print(p,resCnt[p])
			# remove j <= i pattern count 
			# and
			# calc summary value((Ax * Ay cnt) * (Ax * Ay)))
			ans += (resCnt[p] // 2) * p

	print(ans)

Submission status: Accepted

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

	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)

Submission status: Accepted

Solution code:
import numpy as np


class FastFourierTransform:
	"""FFT
	"""

	def multiplySameArr(self,A):
		""" multiply by FFT(A*A)
		"""
		N=1
		# 1. N > deg(g) + deg(h)となる2の冪乗を選ぶ
		while N < 2*len(A): N *= 2

		# 2. 離散フーリエ変換の計算
		nf = np.zeros(N,int)
		nf[:len(A)]=A
		# nf=np.fft.fft(nf)
		nf=np.fft.rfft(nf)
		nf = nf * nf

		# 3. 離散フーリエ逆変換
		# nf=np.fft.ifft(nf)
		nf=np.fft.irfft(nf)

		# コピー
		ans = np.rint( nf )
		# ans = np.zeros(N,int)
		# for n in range(N):
		# 	# 係数を1/N倍 => np.fftの場合、すでに実行済
		# 	#ans.append( round(nf[n].real / 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 of 1 to P-1 after fft
	order=[1]
	x=rootP
	while x!=1:
		order.append(x)
		x *= rootP
		x %= P

	ORDER_SIZE=len(order)	# 200_002

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

	# count non zero items
	prevCnt = np.zeros(P,int)
	resCnt = np.zeros(P,int)
	for n in range(N):
		if A[n] != 0:
			prevCnt[ where[A[n]] ] += 1
			# remove Ax * Ax pattern count
			resCnt[ (A[n] * A[n]) % P ] -= 1

	# fft to cnt (calc nC2 and all Ax * Ay)
	clsFFT = FastFourierTransform()
	fftRes = clsFFT.multiplySameArr(prevCnt)

	for r in range(len(fftRes)):
		if fftRes[r] != 0:
			# get result(Ax*Ay==resCnt.index) Cnt by summary fftRes
			resCnt[ order[ r % ORDER_SIZE ] ] += fftRes[r]

	ans = 0
	for p in range(P):
		if resCnt[p] != 0:
			# remove j <= i pattern count 
			# and
			# calc summary value((Ax * Ay cnt) * (Ax * Ay)))
			ans += (resCnt[p] // 2) * p

	print(ans)

Submission status: Accepted

Solution code:
import numpy as np


class FastFourierTransform:
	"""FFT
	"""

	def multiplySameArr(self,A):
		""" multiply by FFT(A*A)
		"""
		N=1
		# 1. N > deg(g) + deg(h)となる2の冪乗を選ぶ
		while N < 2*len(A): N *= 2

		# 2. 離散フーリエ変換の計算
		nf = np.zeros(N,int)
		nf[:len(A)]=A
		nf=np.fft.fft(nf)
		nf = nf * nf

		# 3. 離散フーリエ逆変換
		nf=np.fft.ifft(nf)

		# コピー
		ans = np.zeros(N,int)
		for n in range(N):
			# 係数を1/N倍 => np.fftの場合、すでに実行済
			#ans.append( round(nf[n].real / 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 of 1 to P-1 after fft
	order=[1]
	x=rootP
	while x!=1:
		order.append(x)
		x *= rootP
		x %= P

	ORDER_SIZE=len(order)	# 200_002

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

	# count non zero items
	prevCnt = np.zeros(P,int)
	resCnt = np.zeros(P,int)
	for n in range(N):
		if A[n] != 0:
			prevCnt[ where[A[n]] ] += 1
			# remove Ax * Ax pattern count
			resCnt[ (A[n] * A[n]) % P ] -= 1

	# fft to cnt (calc nC2 and all Ax * Ay)
	clsFFT = FastFourierTransform()
	fftRes = clsFFT.multiplySameArr(prevCnt)

	for r in range(len(fftRes)):
		if fftRes[r] != 0:
			# get result(Ax*Ay==resCnt.index) Cnt by summary fftRes
			resCnt[ order[ r % ORDER_SIZE ] ] += fftRes[r]

	ans = 0
	for p in range(P):
		if resCnt[p] != 0:
			# remove j <= i pattern count 
			# and
			# calc summary value((Ax * Ay cnt) * (Ax * Ay)))
			ans += (resCnt[p] // 2) * p

	print(ans)

Submission status: Accepted

Solution code:
import sys
sys.setrecursionlimit(10**7) #再帰関数の上限,10**5以上の場合python
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#2分探索
#bisect_left(l,x), bisect(l,x)#aはソート済みである必要あり。aの中からx未満の要素数を返す。rightだと以下
from collections import deque, defaultdict
#deque(l), pop(), append(x), popleft(), appendleft(x)
#q.rotate(n)で → にn回ローテート
from collections import Counter#文字列を個数カウント辞書に、
#S=Counter(l),S.most_common(x),S.keys(),S.values(),S.items()
from itertools import accumulate,combinations,permutations,product#累積和
#list(accumulate(l))
from heapq import heapify,heappop,heappush
#heapify(q),heappush(q,a),heappop(q) #q=heapify(q)としないこと、返り値はNone
from functools import reduce,lru_cache#pypyでもうごく
#@lru_cache(maxsize = None)#maxsizeは保存するデータ数の最大値、2**nが最も高効率
from decimal import Decimal
# from numpy.fft import rfft, irfft
import numpy as np
def input(): 
    x=sys.stdin.readline()
    return x[:-1] if x[-1]=="\n" else x
def printe(*x):print("## ",*x,file=sys.stderr)
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は64bit整数
    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

    #w.sort(key=itemgetter(1),reverse=True)  #二個目の要素で降順並び替え

    N = int(input())
    #N, K = map(int, input().split())
    A = list(map(int, input().split())) #1行ベクトル
    #L = tuple(int(input()) for i in range(N)) #改行ベクトル
    #S = tuple(tuple(map(int, input().split())) for i in range(N)) #改行行列
    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)
    #printe(X,sum(X))
    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()

Submission status: Accepted

Solution code:
import sys
sys.setrecursionlimit(10**7) #再帰関数の上限,10**5以上の場合python
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#2分探索
#bisect_left(l,x), bisect(l,x)#aはソート済みである必要あり。aの中からx未満の要素数を返す。rightだと以下
from collections import deque, defaultdict
#deque(l), pop(), append(x), popleft(), appendleft(x)
#q.rotate(n)で → にn回ローテート
from collections import Counter#文字列を個数カウント辞書に、
#S=Counter(l),S.most_common(x),S.keys(),S.values(),S.items()
from itertools import accumulate,combinations,permutations,product#累積和
#list(accumulate(l))
from heapq import heapify,heappop,heappush
#heapify(q),heappush(q,a),heappop(q) #q=heapify(q)としないこと、返り値はNone
from functools import reduce,lru_cache#pypyでもうごく
#@lru_cache(maxsize = None)#maxsizeは保存するデータ数の最大値、2**nが最も高効率
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("## ",*x,file=sys.stderr)
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は64bit整数
    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

    #w.sort(key=itemgetter(1),reverse=True)  #二個目の要素で降順並び替え

    N = int(input())
    #N, K = map(int, input().split())
    A = list(map(int, input().split())) #1行ベクトル
    #L = tuple(int(input()) for i in range(N)) #改行ベクトル
    #S = tuple(tuple(map(int, input().split())) for i in range(N)) #改行行列
    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)
    #printe(X,sum(X))
    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()

Submission status: Accepted

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

Submission status: Accepted

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

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

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 << 19
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: Accepted

Solution code:
'''
自宅用PCでの解答
'''
import math
#import numpy as np
import itertools
import queue
import bisect
from collections import deque,defaultdict
import heapq as hpq
from sys import stdin,setrecursionlimit
#from scipy.sparse.csgraph import dijkstra
#from scipy.sparse import csr_matrix
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()

Submission status: Accepted

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

Submission status: Accepted

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

Submission status: Accepted

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