#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef complex<double> CD;
typedef long long LL;
const double PI = acos(-1);

struct FFT {
    int N;
    vector<int> perm;

    void precalculate() {
        perm.resize(N);
        perm[0] = 0;

        for (int k=1; k<N; k<<=1) {
            for (int i=0; i<k; i++) {
                perm[i] <<= 1;
                perm[i+k] = 1 + perm[i];
            }
        }
    }

    void fft(vector<CD> &v, bool invert = false) {
        if (v.size() != perm.size()) {
            N = v.size();
            assert(N && (N&(N-1)) == 0);
            precalculate();
        }

        for (int i=0; i<N; i++)
            if (i < perm[i])
                swap(v[i], v[perm[i]]);

        for (int len = 2; len <= N; len <<= 1) {
            double angle = 2 * PI / len;
            if (invert) angle = -angle;
            CD factor = polar(1.0, angle);

            for (int i=0; i<N; i+=len) {
                CD w(1);
                for (int j=0; j<len/2; j++) {
                    CD x = v[i+j], y = w * v[i+j+len/2];
                    v[i+j] = x+y;
                    v[i+j+len/2] = x-y;
                    w *= factor;
                }
            }
        }
        if (invert)
            for (CD &x : v) x/=N;
    }

    vector<LL> multiply(const vector<LL> &a, const vector<LL> &b) {
        vector<CD> fa(a.begin(), a.end()), fb(b.begin(), b.end());

        int n = 1;
        while (n < a.size()+ b.size())  n<<=1;
        fa.resize(n);
        fb.resize(n);

        fft(fa);
        fft(fb);
        for (int i=0; i<n; i++) fa[i] *= fb[i];
        fft(fa, true);

        vector<LL> ans(n);
        for (int i=0; i<n; i++)
            ans[i] = round(fa[i].real());
        return ans;
    }
};

const LL p = 200003;
int found[p];
LL ghat[p];

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

//    LL x = generator(p);
    LL x = 2;
//    cout << x << endl;

    memset(found, -1, sizeof found);
    LL y = 1;
    for (int k = 0; k < p-1; k++) {
        ghat[k] = y;
        assert(found[y]==-1);
        found[y] = k;
        y = (y*x)%p;
    }

    LL singles = 0;

    int n;
    cin >> n;
    vector< LL >v(p-1, 0);
    for (int i = 0; i < n; i++) {
        LL x;
        cin >> x;
        if (x==0) continue;
        singles += (x*x)%p;
        v[found[x]]++;
    }

    FFT fft;
    vector< LL >w = fft.multiply(v, v);

    LL ans = 0;
    for (int i = 0; i < w.size(); i++) {
        ans += ghat[i%(p-1)]*w[i];
    }

    ans = (ans-singles)/2;

    cout << ans << "\n";

    return 0;
}

### Prompt

Construct a cpp code solution to the problem outlined:

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>
using namespace std;
typedef long long LL;
typedef complex<double> CD;
typedef long long LL;
const double PI = acos(-1);

struct FFT {
    int N;
    vector<int> perm;

    void precalculate() {
        perm.resize(N);
        perm[0] = 0;

        for (int k=1; k<N; k<<=1) {
            for (int i=0; i<k; i++) {
                perm[i] <<= 1;
                perm[i+k] = 1 + perm[i];
            }
        }
    }

    void fft(vector<CD> &v, bool invert = false) {
        if (v.size() != perm.size()) {
            N = v.size();
            assert(N && (N&(N-1)) == 0);
            precalculate();
        }

        for (int i=0; i<N; i++)
            if (i < perm[i])
                swap(v[i], v[perm[i]]);

        for (int len = 2; len <= N; len <<= 1) {
            double angle = 2 * PI / len;
            if (invert) angle = -angle;
            CD factor = polar(1.0, angle);

            for (int i=0; i<N; i+=len) {
                CD w(1);
                for (int j=0; j<len/2; j++) {
                    CD x = v[i+j], y = w * v[i+j+len/2];
                    v[i+j] = x+y;
                    v[i+j+len/2] = x-y;
                    w *= factor;
                }
            }
        }
        if (invert)
            for (CD &x : v) x/=N;
    }

    vector<LL> multiply(const vector<LL> &a, const vector<LL> &b) {
        vector<CD> fa(a.begin(), a.end()), fb(b.begin(), b.end());

        int n = 1;
        while (n < a.size()+ b.size())  n<<=1;
        fa.resize(n);
        fb.resize(n);

        fft(fa);
        fft(fb);
        for (int i=0; i<n; i++) fa[i] *= fb[i];
        fft(fa, true);

        vector<LL> ans(n);
        for (int i=0; i<n; i++)
            ans[i] = round(fa[i].real());
        return ans;
    }
};

const LL p = 200003;
int found[p];
LL ghat[p];

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

//    LL x = generator(p);
    LL x = 2;
//    cout << x << endl;

    memset(found, -1, sizeof found);
    LL y = 1;
    for (int k = 0; k < p-1; k++) {
        ghat[k] = y;
        assert(found[y]==-1);
        found[y] = k;
        y = (y*x)%p;
    }

    LL singles = 0;

    int n;
    cin >> n;
    vector< LL >v(p-1, 0);
    for (int i = 0; i < n; i++) {
        LL x;
        cin >> x;
        if (x==0) continue;
        singles += (x*x)%p;
        v[found[x]]++;
    }

    FFT fft;
    vector< LL >w = fft.multiply(v, v);

    LL ans = 0;
    for (int i = 0; i < w.size(); i++) {
        ans += ghat[i%(p-1)]*w[i];
    }

    ans = (ans-singles)/2;

    cout << ans << "\n";

    return 0;
}

```