#include <bits/stdc++.h>
#define rep(i, a, b) for(int i = a; i < b; i++)
#define all(x) begin(x), end(x)
#define sz(x) (int)(x).size()
using namespace std;
using ll = long long;

const int P = 200003;
const int pr = 2;

int idx[P + 5];

typedef complex<double> C;
typedef vector<double> vd;
typedef vector<int> vi;

void fft(vector<C>& a) {
    int n = sz(a), L = 31 - __builtin_clz(n);
    static vector<complex<long double>> R(2, 1);
    static vector<C> rt(2, 1);  // (^ 10% faster if double)
    for (static int k = 2; k < n; k *= 2) {
        R.resize(n); rt.resize(n);
        auto x = polar(1.0L, acos(-1.0L) / k);
        rep(i,k,2*k) rt[i] = R[i] = i&1 ? R[i/2] * x : R[i/2];
    }
    vi rev(n);
    rep(i,0,n) rev[i] = (rev[i / 2] | (i & 1) << L) / 2;
    rep(i,0,n) if (i < rev[i]) swap(a[i], a[rev[i]]);
    for (int k = 1; k < n; k *= 2)
        for (int i = 0; i < n; i += 2 * k) rep(j,0,k) {
            // C z = rt[j+k] * a[i+j+k]; // (25% faster if hand-rolled)  /// include-line
            auto x = (double *)&rt[j+k], y = (double *)&a[i+j+k];        /// exclude-line
            C z(x[0]*y[0] - x[1]*y[1], x[0]*y[1] + x[1]*y[0]);           /// exclude-line
            a[i + j + k] = a[i + j] - z;
            a[i + j] += z;
        }
}

vd conv(const vd& a, const vd& b) {
    if (a.empty() || b.empty()) return {};
    vd res(sz(a) + sz(b) - 1);
    int L = 32 - __builtin_clz(sz(res)), n = 1 << L;
    vector<C> in(n), out(n);
    copy(all(a), begin(in));
    rep(i,0,sz(b)) in[i].imag(b[i]);
    fft(in);
    for (C& x : in) x *= x;
    rep(i,0,n) out[i] = in[-i & (n - 1)] - conj(in[i]);
    fft(out);
    rep(i,0,sz(res)) res[i] = imag(out[i]) / (4 * n);
    return res;
}

ll mpow(ll b, ll e) {
    ll res = 1;
    for(ll k = 1; k <= e; k *= 2) {
        if(k & e) res = (res * b) % P;
        b = (b * b) % P;
    }
    return b;
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(0);

    ll cur = 1;
    for(int i = 0; i < P - 1; i++) {
        idx[cur] = i;
        cur = (cur * pr) % P;
    }
    //for(int i = 1; i < P; i++)
    //    cout << idx[i] << " ";

    vd A(P - 1);
    int n;
    cin >> n;
    ll ans = 0;
    for(int i = 0; i < n; i++) {
        ll x;
        cin >> x;
        if(!x)
            continue;
        ans -= ((x * x) % P);
        A[idx[x]]++;
    }
    auto res = conv(A, A);
    for(int i = 0; i < P - 1; i++) {
        res[i] += res[i + P - 1];
        int cnt = round(res[i]);
        //if(cnt)
        //    cout << i << " " << cnt << '\n';
    }

    cur = 1;
    for(int i = 0; i < P - 1; i++) {
        ans += cur * round(res[i]);
        cur = (cur * pr) % P;
    }
    cout << ans / 2 << '\n';
}
### Prompt

Please formulate a CPP solution to the following problem:

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

Examples

Input

4
2019 0 2020 200002


Output

474287


Input

5
1 1 2 2 100000


Output

600013

### Response

```cpp
#include <bits/stdc++.h>
#define rep(i, a, b) for(int i = a; i < b; i++)
#define all(x) begin(x), end(x)
#define sz(x) (int)(x).size()
using namespace std;
using ll = long long;

const int P = 200003;
const int pr = 2;

int idx[P + 5];

typedef complex<double> C;
typedef vector<double> vd;
typedef vector<int> vi;

void fft(vector<C>& a) {
    int n = sz(a), L = 31 - __builtin_clz(n);
    static vector<complex<long double>> R(2, 1);
    static vector<C> rt(2, 1);  // (^ 10% faster if double)
    for (static int k = 2; k < n; k *= 2) {
        R.resize(n); rt.resize(n);
        auto x = polar(1.0L, acos(-1.0L) / k);
        rep(i,k,2*k) rt[i] = R[i] = i&1 ? R[i/2] * x : R[i/2];
    }
    vi rev(n);
    rep(i,0,n) rev[i] = (rev[i / 2] | (i & 1) << L) / 2;
    rep(i,0,n) if (i < rev[i]) swap(a[i], a[rev[i]]);
    for (int k = 1; k < n; k *= 2)
        for (int i = 0; i < n; i += 2 * k) rep(j,0,k) {
            // C z = rt[j+k] * a[i+j+k]; // (25% faster if hand-rolled)  /// include-line
            auto x = (double *)&rt[j+k], y = (double *)&a[i+j+k];        /// exclude-line
            C z(x[0]*y[0] - x[1]*y[1], x[0]*y[1] + x[1]*y[0]);           /// exclude-line
            a[i + j + k] = a[i + j] - z;
            a[i + j] += z;
        }
}

vd conv(const vd& a, const vd& b) {
    if (a.empty() || b.empty()) return {};
    vd res(sz(a) + sz(b) - 1);
    int L = 32 - __builtin_clz(sz(res)), n = 1 << L;
    vector<C> in(n), out(n);
    copy(all(a), begin(in));
    rep(i,0,sz(b)) in[i].imag(b[i]);
    fft(in);
    for (C& x : in) x *= x;
    rep(i,0,n) out[i] = in[-i & (n - 1)] - conj(in[i]);
    fft(out);
    rep(i,0,sz(res)) res[i] = imag(out[i]) / (4 * n);
    return res;
}

ll mpow(ll b, ll e) {
    ll res = 1;
    for(ll k = 1; k <= e; k *= 2) {
        if(k & e) res = (res * b) % P;
        b = (b * b) % P;
    }
    return b;
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(0);

    ll cur = 1;
    for(int i = 0; i < P - 1; i++) {
        idx[cur] = i;
        cur = (cur * pr) % P;
    }
    //for(int i = 1; i < P; i++)
    //    cout << idx[i] << " ";

    vd A(P - 1);
    int n;
    cin >> n;
    ll ans = 0;
    for(int i = 0; i < n; i++) {
        ll x;
        cin >> x;
        if(!x)
            continue;
        ans -= ((x * x) % P);
        A[idx[x]]++;
    }
    auto res = conv(A, A);
    for(int i = 0; i < P - 1; i++) {
        res[i] += res[i + P - 1];
        int cnt = round(res[i]);
        //if(cnt)
        //    cout << i << " " << cnt << '\n';
    }

    cur = 1;
    for(int i = 0; i < P - 1; i++) {
        ans += cur * round(res[i]);
        cur = (cur * pr) % P;
    }
    cout << ans / 2 << '\n';
}
```