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: Runtime Error

Solution code:
#include <iostream>
#include <algorithm>
#include <vector>
#include <cmath>
#include <complex>

template <int MOD>
struct ModInt {
    using lint = long long;
    int val;

    // constructor
    ModInt(lint v = 0) : val(v % MOD) {
        if (val < 0) val += MOD;
    };

    // unary operator
    ModInt operator+() const { return ModInt(val); }
    ModInt operator-() const { return ModInt(MOD - val); }
    ModInt inv() const { return this->pow(MOD - 2); }

    // arithmetic
    ModInt operator+(const ModInt& x) const { return ModInt(*this) += x; }
    ModInt operator-(const ModInt& x) const { return ModInt(*this) -= x; }
    ModInt operator*(const ModInt& x) const { return ModInt(*this) *= x; }
    ModInt operator/(const ModInt& x) const { return ModInt(*this) /= x; }
    ModInt pow(lint n) const {
        auto x = ModInt(1);
        auto b = *this;
        while (n > 0) {
            if (n & 1) x *= b;
            n >>= 1;
            b *= b;
        }
        return x;
    }

    // compound assignment
    ModInt& operator+=(const ModInt& x) {
        if ((val += x.val) >= MOD) val -= MOD;
        return *this;
    }
    ModInt& operator-=(const ModInt& x) {
        if ((val -= x.val) < 0) val += MOD;
        return *this;
    }
    ModInt& operator*=(const ModInt& x) {
        val = lint(val) * x.val % MOD;
        return *this;
    }
    ModInt& operator/=(const ModInt& x) { return *this *= x.inv(); }

    // compare
    bool operator==(const ModInt& b) const { return val == b.val; }
    bool operator!=(const ModInt& b) const { return val != b.val; }
    bool operator<(const ModInt& b) const { return val < b.val; }
    bool operator<=(const ModInt& b) const { return val <= b.val; }
    bool operator>(const ModInt& b) const { return val > b.val; }
    bool operator>=(const ModInt& b) const { return val >= b.val; }

    // I/O
    friend std::istream& operator>>(std::istream& is, ModInt& x) noexcept {
        lint v;
        is >> v;
        x = v;
        return is;
    }
    friend std::ostream& operator<<(std::ostream& os, const ModInt& x) noexcept { return os << x.val; }
};

template <int K>
struct FastFourierTransform {
    using Real = double;

    struct cplx {
        Real real, imag;

        explicit cplx(Real real = 0, Real imag = 0)
            : real(real), imag(imag) {}

        cplx operator+(const cplx& other) const {
            return cplx(real + other.real, imag + other.imag);
        }
        cplx operator-(const cplx& other) const {
            return cplx(real - other.real, imag - other.imag);
        }
        cplx operator*(const cplx& other) const {
            return cplx(real * other.real - imag * other.imag,
                        real * other.imag + imag * other.real);
        }
        cplx conj() const { return cplx(real, -imag); }
    };
    using cplxs = std::vector<cplx>;

    static constexpr Real PI = 3.14159265358979323846L;
    const cplx unit_i{0, 1};

    cplxs zetas;

    explicit FastFourierTransform() : zetas(K) {
        for (int i = 0; i < K; ++i) {
            zetas[i] = cplx(std::cos(PI * 2 / (1 << i)),
                            std::sin(PI * 2 / (1 << i)));
        }
    }

    // ceil(log_2 n)
    static int clog2(int n) {
        int k = 0;
        while ((1 << k) < n) ++k;
        return k;
    }

    // ceil(log_4 n)
    static int clog4(int n) {
        int k = 0;
        while ((1 << (k * 2)) < n) ++k;
        return k;
    }

    // 4-radix cooley-tukey algorithm without bit reverse
    // the size of f must be a power of 4
    void fft(cplxs& f) const {
        int n = f.size();

        for (int m = n >> 2; m >= 1; m >>= 2) {
            auto zeta = zetas[clog2(m) + 2];
            cplx zeta1(1), zeta2, zeta3;

            for (int p = 0; p < m; ++p) {
                zeta2 = zeta1 * zeta1;
                zeta3 = zeta1 * zeta2;

                for (int s = 0; s < n; s += (m << 2)) {
                    auto w = f[s + p + m * 0],
                         x = f[s + p + m * 1],
                         y = f[s + p + m * 2],
                         z = f[s + p + m * 3];

                    auto wyp = w + y,
                         wym = w - y,
                         xzp = x + z,
                         xzm = (x - z) * unit_i;

                    f[s + p + m * 0] = wyp + xzp;
                    f[s + p + m * 1] = (wym + xzm) * zeta1;
                    f[s + p + m * 2] = (wyp - xzp) * zeta2;
                    f[s + p + m * 3] = (wym - xzm) * zeta3;
                }
                zeta1 = zeta1 * zeta;
            }
        }
    }

    // the inverse of above function
    void ifft(cplxs& f) const {
        int n = f.size();

        for (int m = 1; m <= (n >> 2); m <<= 2) {
            auto zeta = zetas[clog2(m) + 2].conj();
            cplx zeta1 = cplx(1), zeta2, zeta3;

            for (int p = 0; p < m; ++p) {
                zeta2 = zeta1 * zeta1;
                zeta3 = zeta1 * zeta2;

                for (int s = 0; s < n; s += (m << 2)) {
                    auto w = f[s + p + m * 0],
                         x = f[s + p + m * 1] * zeta1,
                         y = f[s + p + m * 2] * zeta2,
                         z = f[s + p + m * 3] * zeta3;

                    auto wyp = w + y,
                         wym = w - y,
                         xzp = x + z,
                         xzm = (x - z) * unit_i;

                    f[s + p + m * 0] = wyp + xzp;
                    f[s + p + m * 1] = wym - xzm;
                    f[s + p + m * 2] = wyp - xzp;
                    f[s + p + m * 3] = wym + xzm;
                }
                zeta1 = zeta1 * zeta;
            }
        }

        auto ninv = cplx(1.L / n);
        for (auto& x : f) x = x * ninv;
    }

    // main routine
    using lint = long long;
    using lints = std::vector<lint>;

    lints convolute(const lints& f, const lints& g) const {
        int fsz = f.size(),
            gsz = g.size();

        // simple convolution in small cases
        if (std::min(fsz, gsz) < 8) {
            lints ret(fsz + gsz - 1, 0);
            for (int i = 0; i < fsz; ++i) {
                for (int j = 0; j < gsz; ++j) {
                    ret[i + j] += f[i] * g[j];
                }
            }
            return ret;
        }

        auto fc = li2cp(f),
             gc = li2cp(g);

        int n = 1 << (clog4(fsz + gsz - 1) * 2);
        fc.resize(n, cplx(0));
        gc.resize(n, cplx(0));

        fft(fc);
        fft(gc);

        for (int i = 0; i < n; ++i) fc[i] = fc[i] * gc[i];

        ifft(fc);

        fc.resize(fsz + gsz - 1);
        return cp2li(fc);
    }

    // lint <-> complex
    cplxs li2cp(const lints& f) const {
        cplxs ret;
        std::transform(f.begin(), f.end(), std::back_inserter(ret),
                       [](auto x) { return cplx(x); });
        return ret;
    }

    lints cp2li(const cplxs& f) const {
        lints ret;
        std::transform(f.begin(), f.end(), std::back_inserter(ret),
                       [](auto x) { return std::llround(x.real); });
        return ret;
    }
};

constexpr int MOD = 200003;
using mint = ModInt<MOD>;
const FastFourierTransform<25> FFT;

using lint = long long;

void solve() {
    // g^log[x] = x
    std::vector<int> log(MOD);
    mint g = 2;
    for (int i = 0; i < MOD - 1; ++i) {
        log[g.pow(i).val] = i;
    }

    std::vector<lint> cnt(MOD - 1, 0);

    int n;
    std::cin >> n;
    while (n--) {
        int a;
        std::cin >> a;
        if (a != 0) ++cnt[log[a]];
    }

    // FFTで畳み込む
    auto res = FFT.convolute(cnt, cnt);

    lint ans = 0;
    for (int i = 0; i < (int)res.size(); ++i) {
        lint num = llround(res[i]);
        ans += g.pow(i).val * num;
    }

    // a_i * a_iを省く
    for (int i = 0; i < (int)cnt.size(); ++i) {
        lint num = llround(cnt[i]);
        ans -= g.pow(i * 2).val * num;
    }

    // i > jを省く
    ans /= 2;

    std::cout << ans << "\n";
}

int main() {
    std::cin.tie(nullptr);
    std::ios::sync_with_stdio(false);

    solve();

    return 0;
}

Submission status: Runtime Error

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
for i in range(p):
    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: Runtime Error

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:
  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: Runtime Error

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):
    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: Runtime Error

Solution code:
import numpy as np
fft = np.fft.fft
ifft = np.fft.ifft
#from scipy.fftpack import fft, ifft

(N, ) = [int(n) for n in input().split()]
A = [int(n) for n in input().split()]

P = 200003
PP = 524288

#P = 13
#PP = 32

log = [0 for _ in range(P)] # log2x (1 <= x <= P-1)
exp = [0 for _ in range(P)] # 2^x (1 <= x <= P-1)

base = 2

a = 1
for i in range(P):
    if i != 0:
        #print("log {} = {}, exp {} = {}".format(a, i, i, a))
        exp[i] = a
        log[a] = i
    a *= base
    a %= P

c = [0 for _ in range(P)]
for a in A:
    c[a] += 1

d = [c[exp[i]] for i in range(P)]
dtype = np.float128
d = np.array(d[1:], dtype=dtype)

d1 = np.zeros((PP,), dtype=dtype)
d2 = np.zeros((PP,), dtype=dtype)
d1[0:P-1] = d
d2[0:P-1] = d
d2[-(P-1):] = d
D1 = fft(d1)
D2 = fft(d2)
D = D1 * D2
dd = ifft(D)[0:P-1]
#print(np.abs(dd))

ans = 0
for k in range(P-1):
    ans += dd[(k-1+P-1)%(P-1)] * exp[k+1]

ans = int(round(abs(ans)))

for a in A:
    ans -= (a * a) % P

assert(ans % 2 == 0)
print(ans // 2)

Submission status: Runtime Error

Solution code:
import numpy as np
#fft = np.fft.fft
#ifft = np.fft.ifft
from scipy.fftpack import fft, ifft

(N, ) = [int(n) for n in input().split()]
A = [int(n) for n in input().split()]

P = 200003
PP = 524288

#P = 13
#PP = 32

log = [0 for _ in range(P)] # log2x (1 <= x <= P-1)
exp = [0 for _ in range(P)] # 2^x (1 <= x <= P-1)

base = 2

a = 1
for i in range(P):
    if i != 0:
        #print("log {} = {}, exp {} = {}".format(a, i, i, a))
        exp[i] = a
        log[a] = i
    a *= base
    a %= P

c = [0 for _ in range(P)]
for a in A:
    c[a] += 1

d = [c[exp[i]] for i in range(P)]
d = np.array(d[1:], dtype=np.float64)

d1 = np.zeros((PP,), dtype=np.float64)
d2 = np.zeros((PP,), dtype=np.float64)
d1[0:P-1] = d
d2[0:P-1] = d
d2[-(P-1):] = d
D1 = fft(d1)
D2 = fft(d2)
D = D1 * D2
dd = ifft(D)[0:P-1]
#print(np.abs(dd))

ans = 0
for k in range(P-1):
    ans += dd[(k-1+P-1)%(P-1)] * exp[k+1]

ans = int(round(abs(ans)))

for a in A:
    ans -= (a * a) % P

assert(ans % 2 == 0)
print(ans // 2)

Submission status: Runtime Error

Solution code:
import numpy as np
#fft = np.fft.fft
#ifft = np.fft.ifft
from scipy.fftpack import fft, ifft

(N, ) = [int(n) for n in input().split()]
A = [int(n) for n in input().split()]

P = 200003
#P = 13

log = [0 for _ in range(P)] # log2x (1 <= x <= P-1)
exp = [0 for _ in range(P)] # 2^x (1 <= x <= P-1)

base = 2

a = 1
for i in range(P):
    if i != 0:
        #print("log {} = {}, exp {} = {}".format(a, i, i, a))
        exp[i] = a
        log[a] = i
    a *= base
    a %= P

c = [0 for _ in range(P)]
for a in A:
    c[a] += 1

d = [c[exp[i]] for i in range(P)]

d = np.array(d[1:], dtype=np.float64)
D = fft(d)
D = D ** 2
dd = ifft(D)

ans = 0
for k in range(P-1):
    ans += dd[(k-1+P-1)%(P-1)] * exp[k+1]

ans = int(round(abs(ans)))

for a in A:
    ans -= (a * a) % P

assert(ans % 2 == 0)
print(ans // 2)

Submission status: Runtime Error

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: Runtime Error

Solution code:
N = int(input())
a = list(map(int,input().split()))
sum = 0
P =200003
d = 0
count = N-1
temp = 0
a.sort()
for num in a:
	d += num

for i in range(N):
	b = P //a[i]
	for j in range(count+1):
		if a[count-j] <= b:
			break
		temp = count-j
		count = temp
	temp = d
	for q in range(count,N):
		sum += a[i]*a[count-q]%P
		temp -= a[q]
	sum += a[i]*temp

print(sum)

Submission status: Runtime Error

Solution code:
N = int(input())
a = list(map(int,input().split()))
sum = 0
P =200003
d = 0
count = N-1
temp = 0
a.sort()
for num in a:
	d += num
for num in a:
	sum += num*(d-num)

for i in range(N):
	b = P //a[i]
	for j in range(count):
		if a[count-j] <= b:
			break
		sum -= a[i]*a[count-j]
		sum += a[i]*a[count-j]%P
		temp = count-j
	count = temp
	
	

print(sum)


#9872..8332

Submission status: Runtime Error

Solution code:
import sys;input=sys.stdin.readline
N, = map(int, input().split())
S = []
for _ in range(N):
    s = input().strip()
    S.append(s)
S.sort(key=lambda x:len(x))
D = dict()
dd = set()
R = 0
for s in S:
    l = len(s)
    d = [0]*26
    for c in dd:
        d[c] += 1
    mm = 0
    tmp=0
    for i in range(l-1, -1, -1):
        cc = ord(s[i]) - 97
        tmp=26*tmp+cc
        R += d[cc]
        d[cc] = 0
        if tmp in D:
            for c in D[tmp]:
                d[c] += 1
        if i == 0:
            f = cc
        if i == 1:
            nn = tmp
        mm += 1
    if l != 1:
        if nn not in D:
            D[nn] = set()
        D[nn].add(f)
    else:
        dd.add(f)
print(R)

Submission status: Runtime Error

Solution code:
import pandas as pd
P = 200003
N = int(input())
s = pd.Series(map(int,input().split(" ")))

ans = 0
for i in range(N):
    ans += sum(s[i]*s[i+1:]%P)

print(ans)

Submission status: Runtime Error

Solution code:
N = int(input())
a = list(map(int,input().split()))
sum = 0
P =200003
d = 0
for num in a:
	d += num
for num in a:
	sum += num*(d-num)

for i in range(N):
	b = P // a[i]
	list =[a[j] for j in range(i+1,N) if a[j]>b]
	for x in list:
		sum -= (a[i]*x//P)*P
	

print(sum)

Submission status: Runtime Error

Solution code:
import numpy as np
def main():
  n=int(input())
  a=np.array(list(map(int,(input()).split())))
  p=200003
  aouter = np.outer(a, a)
  adot=a*a
  aouter=aouter%p
  adot=adot%p
  print((np.sum(aouter)-np.sum(adot))//2)
main()

Submission status: Runtime Error

Solution code:
import numpy as np
def main():
  n=int(input())
  a=np.array(list(map(int,(input().replace(" 0 "," ")).split())))
  p=200003
  aouter = np.outer(a, a)
  adot=a*a
  aouter=aouter%p
  adot=adot%p
  print((np.sum(aouter)-np.sum(adot))//2)
main()

Submission status: Runtime Error

Solution code:
P = 200003
N = int(input())
s = list(map(int,input().split(" ")))
if 0 in s:
	s.remove(0)
 
ans = 0
while len(s) != 0:
	ans += sum(list(map(lambda x: (x * s[0])%P ,s[1:])))
    s.pop(0)
print(ans)

Submission status: Runtime Error

Solution code:
import itertools
n = int(input())
arr = list(map(int, input().split()))
arr.remove(0)
s = 0
P = 200003

for pair in itertools.combinations(arr, r=2):
    s += ((pair[0]*pair[1])%P)

print(s)

Submission status: Runtime Error

Solution code:
P = 200003
N = int(input())
s = list(map(int,input().split(" ")))
s = s.remove(0)

ans = 0
while len(s) != 0:
    ans += sum(list(map(lambda x: (x * s[0])%P ,s[1:])))
    s.pop(0)
print(ans)

Submission status: Runtime Error

Solution code:
// Dmitry _kun_ Sayutin (2020)

#include <bits/stdc++.h>

using std::cin;
using std::cout;
using std::cerr;

using std::vector;
using std::map;
using std::array;
using std::set;
using std::string;
using std::queue;

using std::pair;
using std::make_pair;

using std::tuple;
using std::make_tuple;
using std::get;

using std::min;
using std::abs;
using std::max;
using std::swap;

using std::unique;
using std::sort;
using std::generate;
using std::reverse;
using std::min_element;
using std::max_element;

#ifdef LOCAL
#define LASSERT(X) assert(X)
#else
#define LASSERT(X) {}
#endif

template <typename T>
T input() {
    T res;
    cin >> res;
    LASSERT(cin);
    return res;
}

template <typename IT>
void input_seq(IT b, IT e) {
    std::generate(b, e, input<typename std::remove_reference<decltype(*b)>::type>);
}

#define SZ(vec)         int((vec).size())
#define ALL(data)       data.begin(),data.end()
#define RALL(data)      data.rbegin(),data.rend()
#define TYPEMAX(type)   std::numeric_limits<type>::max()
#define TYPEMIN(type)   std::numeric_limits<type>::min()

const int LOG = 20;
const int N = (1 << LOG);

typedef std::complex<double> cd;

int rev[N];
cd W[N];

void precalc() {
    const double pi = std::acos(-1);
    for (int i = 0; i != N; ++i)
        W[i] = cd(std::cos(2 * pi * i / N), std::sin(2 * pi * i / N));

    int last = 0;
    for (int i = 1; i != N; ++i) {
        if (i == (2 << last))
            ++last;

        rev[i] = rev[i ^ (1 << last)] | (1 << (LOG - 1 - last));
    }
}

void fft(vector<cd>& a) {
    for (int i = 0; i != N; ++i)
        if (i < rev[i])
            std::swap(a[i], a[rev[i]]);

    for (int lvl = 0; lvl != LOG; ++lvl)
        for (int start = 0; start != N; start += (2 << lvl))
            for (int pos = 0; pos != (1 << lvl); ++pos) {
                cd x = a[start + pos];
                cd y = a[start + pos + (1 << lvl)];

                y *= W[pos << (LOG - 1 - lvl)];
                
                a[start + pos] = x + y;
                a[start + pos + (1 << lvl)] = x - y;
            }
}

void inv_fft(vector<cd>& a) {
    fft(a);
    std::reverse(a.begin() + 1, a.end());
    
    for (cd& elem: a)
        elem /= N;
}

#define int64_t __int128_t

int main() {
    std::iostream::sync_with_stdio(false);
    cin.tie(nullptr);
    cout.tie(nullptr);

    // code here
    int root = 2;
    int mod = 200 * 1000 + 3;

    vector<int> logg(mod, -100);
    vector<int> poww(mod, -100);

    int64_t sum = 0;
    
    if (1) {
        int pw = 1;

        for (int i = 0; i <= mod - 2; ++i) {
            logg[pw] = i;
            poww[i] = pw;
            
            pw = (pw * 1ll * root) % mod;
        }
    }
    
    int n = input<int>();

    vector<cd> arr(1 << LOG);

    vector<int> in(n);
    for (int i = 0; i < n; ++i) {
        int val = input<int>();
        in[i] = val;
        
        if (val)
            arr[logg[val]] += 1;
    }

    precalc();
    fft(arr);
    for (int i = 0; i < SZ(arr); ++i)
        arr[i] *= arr[i];

    inv_fft(arr);

    for (int i = 0; i < SZ(arr); ++i) {
        sum += poww[i % (mod - 1)] * int64_t(arr[i].real());
    }

    for (int i = 0; i < n; ++i)
        sum -= (in[i] * 1ll * in[i]) % mod;

    sum /= 2;
    
    cout << (long long)sum << "\n";
    
    return 0;
}

Submission status: Runtime Error

Solution code:
n = int(input())
arr = list(map(int, input().split()))
ans = 0
while True:
    if 0 in arr:
        arr.remove(0)
    else:
        break
for i in range(n):
    for j in range(i + 1, n):
        ans += (arr[i] * arr[j]) % 200003
print(ans)

Submission status: Runtime Error

Solution code:
n = int(input())
arr = list(map(int, input().split()))
for i in range(n):
  for j in range(i + 1, n):
    ans += (arr[i] * arr[j]) % 200003
print(ans)