Buggy python submission (verdict: Time Limit Exceeded):
class Convolution:
  def __init__(s, MOD):
    s.MOD = MOD
    s.PrimitiveRoot = s.primitive_root(MOD)
    s._BIT_TABLE = [0] * 64
    _BIT_HASH = 0x03F566ED27179461
    for i in range(64):
      s._BIT_TABLE[_BIT_HASH >> 58] = i
      _BIT_HASH <<= 1
      _BIT_HASH &= 0xFFFFFFFFFFFFFFFF

  def convolution(s, a, b):
    n, m = len(a), len(b)
    if not n or not m: return []
    if min(n, m) <= 60:
      if n < m:
        n, m = m, n
        a, b = b, a
      ans = [0] * (n + m - 1)
      for i in range(n):
        for j in range(m):
          ans[i + j] = (ans[i + j] + a[i] * b[j]) % s.MOD
      return ans
    z = 1 << s.ceil_pow2(n + m - 1)
    a = a + [0] * (z - n)
    a = s.butterfly(a)
    b = b + [0] * (z - m)
    b = s.butterfly(b)
    for i in range(z): a[i] = a[i] * b[i] % s.MOD
    a = s.butterfly_inv(a)
    iz = pow(z, s.MOD - 2, s.MOD)
    for i in range(n + m - 1): a[i] = a[i] * iz % s.MOD
    return a[:n + m - 1]

  def convolution_ll(s, a, b):
    n, m = len(a), len(b)
    if not n or not m: return []
    # 64bitに収める
    B = pow(2, 64)
    MOD1 = 754974721
    MOD2 = 167772161
    MOD3 = 469762049
    M2M3 = MOD2 * MOD3 % B
    M1M3 = MOD1 * MOD3 % B
    M1M2 = MOD1 * MOD2 % B
    M1M2M3 = MOD1 * MOD2 * MOD3 % B
    i1 = s.inv_gcd(MOD2 * MOD3, MOD1)[1]
    i2 = s.inv_gcd(MOD1 * MOD3, MOD2)[1]
    i3 = s.inv_gcd(MOD1 * MOD2, MOD3)[1]
    conv1 = Convolution(MOD1)
    conv2 = Convolution(MOD2)
    conv3 = Convolution(MOD3)

    c1 = conv1.convolution(a, b)
    c2 = conv2.convolution(a, b)
    c3 = conv3.convolution(a, b)

    c = [0] * (n + m - 1)

    for i in range(n + m - 1):
      x = 0
      x += (c1[i] * i1) % MOD1 * M2M3 % B
      x = (x + (c2[i] * i2) % MOD2 * M1M3) % B
      x = (x + (c3[i] * i3) % MOD3 * M1M2) % B
      diff = c1[i] - (x % MOD1)
      if diff < 0: diff += MOD1
      offset = [0, 0, M1M2M3, 2 * M1M2M3, 3 * M1M2M3]
      x -= offset[diff % 5]
      c[i] = x
    return c

  def safe_mod(s, x, m):
    x %= m
    if x < 0: x += m
    return x

  def butterfly(s, a):
    g = s.PrimitiveRoot
    n = len(a)
    h = s.ceil_pow2(n)

    first = True
    sum_e = [0] * 30
    if first:
      first = False
      es, ies = [0] * 30, [0] * 30
      cnt2 = s.bsf(s.MOD - 1)
      e = pow(g, ((s.MOD - 1) >> cnt2), s.MOD)
      ie = pow(e, s.MOD - 2, s.MOD)
      for i in range(cnt2, 1, -1):
        es[i - 2] = e
        ies[i - 2] = ie
        e = e * e % s.MOD
        ie = ie * ie % s.MOD
      now = 1
      for i in range(cnt2 - 2):
        sum_e[i] = es[i] * now % s.MOD
        now = now * ies[i] % s.MOD
    for ph in range(1, h + 1):
      w, p = 1 << (ph - 1), 1 << (h - ph)
      now = 1
      for ss in range(w):
        offset = ss << (h - ph + 1)
        for i in range(p):
          l = a[i + offset]
          r = a[i + offset + p] * now % s.MOD
          a[i + offset] = (l + r) % s.MOD
          a[i + offset + p] = (l - r) % s.MOD
        now = now * sum_e[s.bsf(~ss)] % s.MOD
    return a

  def butterfly_inv(s, a):
    g = s.PrimitiveRoot
    n = len(a)
    h = s.ceil_pow2(n)

    first = True
    sum_ie = [0] * 30
    if first:
      first = False
      es, ies = [0] * 30, [0] * 30
      cnt2 = s.bsf(s.MOD - 1)
      e = pow(g, ((s.MOD - 1) >> cnt2), s.MOD)
      ie = pow(e, s.MOD - 2, s.MOD)
      for i in range(cnt2, 1, -1):
        es[i - 2] = e
        ies[i - 2] = ie
        e = e * e % s.MOD
        ie = ie * ie % s.MOD
      now = 1
      for i in range(cnt2 - 2):
        sum_ie[i] = ies[i] * now % s.MOD
        now = now * es[i] % s.MOD
    for ph in range(h, 0, -1):
      w, p = 1 << (ph - 1), 1 << (h - ph)
      inow = 1
      for ss in range(w):
        offset = ss << (h - ph + 1)
        for i in range(p):
          l = a[i + offset]
          r = a[i + offset + p]
          a[i + offset] = (l + r) % s.MOD
          a[i + offset + p] = (s.MOD + l - r) * inow % s.MOD
        inow = inow * sum_ie[s.bsf(~ss)] % s.MOD
    return a

  def ceil_pow2(s, n):
    x = 0
    while (1 << x) < n: x += 1
    return x

  def bsf(s, n):
    if n == 0: return 64
    x = n
    y = x & -x
    i = ((y * 0x03F566ED27179461) & 0xFFFFFFFFFFFFFFFF) >> 58
    return s._BIT_TABLE[i]

  def prime(s, x):
    L, U = [], []
    i = 1
    while i * i <= x:
      if x % i == 0:
        L.append(i)
        if i != x // i:
          U.append(x // i)
      i += 1
    return L + U[::-1]

  def primitive_root(s, x):
    L = s.prime(x - 1)
    a, i = 2, 0
    k = len(L) - 1
    while i != k:
      if pow(a, L[i], x) == 1:
        a += 1
        i = 0
      else:
        i += 1
    return a

  def inv_gcd(s, a, b):
    a %= b
    if a == 0: return (b, 0)
    v, t = b, a
    m0, m1 = 0, 1
    while t:
      u = v // t
      v -= t * u
      m0 -= m1 * u
      v, t = t, v
      m0, m1 = m1, m0
    if m0 < 0: m0 += b // v
    return (v, m0)

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

r = 2
EXP = [0] * MOD
LOG = [0] * MOD
EXP[0] = 1
# r ^ i % MOD = EXP[i]
# logr(i) = LOG[i]
for i in range(1, MOD):
  EXP[i] = EXP[i - 1] * r % MOD
  LOG[EXP[i]] = i

b = [0] * MOD
for n in a:
  # log0は無限、実際には不要なので無視
  if n == 0: continue 
  b[LOG[n]] += 1

conv = Convolution(MOD)

# x * y = z => log(x * y) = log(z) => log(x) + log(y) = log(z)
# c[i] <- log(x) + log(y) が i となる個数
c = conv.convolution_ll(b, b)

ans = 0
t = 1
for n in c:
  # log(z)をzに戻しつつ個数を掛ける
  ans += n * t
  t = t * r % MOD

# i * i が含まれているので引く
for i in a:
  ans -= i * i % MOD

# i * j と j * i が両方含まれているので割る
ans //= 2

print(ans)


Fixed submission:
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
# r ^ i % MOD = EXP[i]
# logr(i) = LOG[i]
for i in range(1, MOD):
  EXP[i] = EXP[i - 1] * r % MOD
  LOG[EXP[i]] = i

b = [0] * MOD
for n in a:
  # log0は無限、実際には不要なので無視
  if n == 0: continue 
  b[LOG[n]] += 1

# x * y = z => log(x * y) = log(z) => log(x) + log(y) = log(z)
# c[i] <- log(x) + log(y) が i となる個数
f = FFT()
c = f.exe(b, b)

ans = 0
t = 1
for n in c:
  # log(z)をzに戻しつつ個数を掛ける
  ans += n * t
  t = t * r % MOD

# i * i が含まれているので引く
for i in a:
  ans -= i * i % MOD

# i * j と j * i が両方含まれているので割る
ans //= 2

print(ans)