// #pragma GCC optimize ("03")
#include <bits/stdc++.h>
using namespace std;

#define pb push_back
#define fastio cin.tie(nullptr); ios_base::sync_with_stdio(false)
typedef long long ll;
typedef vector<int> vi;
typedef pair<int, int> ii;
typedef long double ld;

//BEGIN-LIB:
typedef complex<double> CX;
//FROM-FILE:math/polynomial/fast-fourier-transform.cpp
const double PI = acos(-1.0);

// inv == false:
//   a: coefficients of the polynomial p(x) = 
//   a[0] + a[1]*x + a[2]*x^2 ... a[a.size()-1]*x^(a.size()-1),
//   with degree a.size()-1
//   return: the discrete fourier transformation applied of p(x),
//   i.e., the evaluation of p(x) to the n-th complex roots of unity
// inv == true:
//   applies the inverse fourier transformation, getting back the
//   real coefficients from the polynomial whose evalutaion in the
//   n-th roots of unity yields the complex values in a
// complexity: O(N lg N),
// either the direct, or the inverse transformation
// NOTE: a.size() must be a power of 2
// NOTE: this function could be also plugged-in with
// STL's built-in std::complex type, by defining:
//      typedef complex<double> CX;
// however, an user-defined complex type turns out to be
// more efficient in practice
void fft (vector<CX>& a, bool inv) {
    int n = a.size();
    for (int i = 0; i < n; i++) {
        int j = 0;
        for (int bn = n, bi = i; bn >>= 1; bi >>= 1) {
            j <<= 1;
            j |= bi & 1;
        }
        if (i < j) swap(a[i], a[j]);
    }
    for (int h = 1, s = 2; h < n; h = s, s <<= 1) {
        double t = (inv ? -2.0 : 2.0)*PI/s;
        CX wn(cos(t), sin(t));
        for (int l = 0, m = h; l < n; l += s, m += s) {
            CX w = 1.0;
            for (int i = l, j = m; i < m; i++, j++) {
                CX u = a[i], v = a[j]*w;
                a[i] = u + v;
                a[j] = u - v;
                w *= wn;
            }
        }
    }
    if (inv) {
        for (CX& z : a) z /= n;
    }
}
//FROM-FILE:math/polynomial/polynomial-multiplication.cpp
// a, b: polynomials with integer coefficients
// a(x) = a[0] + a[1]*x + ... + a[a.size()-1]*x^(a.size()-1)
// b(x) = b[0] + b[1]*x + ... + b[b.size()-1]*x^(b.size()-1)
// return: coefficients of c(x) = a(x)*b(x)
// complexity: O(N lg N),
// N = degree of a + degree of b
vector<ll> multiply (const vi& a, const vi& b) {
    vector<ll> c(a.size() + b.size() - 1);
    int n = 1;
    while (n < c.size()) n <<= 1;
    vector<CX> ta(a.begin(), a.end()), tb(b.begin(), b.end()), tc(n);
    ta.resize(n);
    tb.resize(n);
    fft(ta, false);
    fft(tb, false);
    for (int i = 0; i < n; i++) tc[i] = ta[i]*tb[i];
    fft(tc, true);
    for (int i = 0; i < c.size(); i++) c[i] = llrint(tc[i].real());
    return c;
}
//FROM-FILE:math/fast-exponentiation/modular.cpp
// computes a**n mod m
// assumes 0 <= a < m
// complexity: O(lg n)
// NOTE: a and m can't be long long for overflow reasons
int modExp (int a, ll n, int m) {
    int ans = 1;
    while (n) {
        if (n & 1) ans = (ll) ans*a % m;
        a = (ll) a*a % m;
        n >>= 1;
    }
    return ans;
}
//END-LIB:

const int N = 200000, P = 200003, R = 2;

int a[N], dl[P];

int main () {
    fastio;
    int n;
    cin >> n;
    vi freq(P), sub(P - 1);
    int x = 1;
    for (int i = 0; i < P - 1; i++) {
        dl[x] = i;
        x = 2*x % P;
    }
    for (int i = 0; i < n; i++) {
        cin >> a[i];
        if (a[i] == 0) continue;
        int x = dl[a[i]];
        freq[x]++;
        sub[2*x % (P - 1)]++;
    }
    vector<ll> v = multiply(freq, freq);
    for (int i = P - 1; i < v.size(); i++) {
        v[i % (P - 1)] += v[i];
    }
    v.resize(P - 1);
    ll ans = 0;
    x = 1;
    for (int i = 0; i < P - 1; i++) {
        ll y = v[i] - sub[i] >> 1;
        ans += y*x;
        x = R*x % P;
    }
    cout << ans << '\n';
}

Measured memory_bytes: 37204