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

Solution (in CPP):

Solution:
#include <vector>
#include <complex>
#include <cmath>
using namespace std;

vector<complex<double>> FFT(vector<complex<double>> v, bool inv)
{
    double pi = acos(-1.);

    int n = int(v.size());
    int b = 0;
    while (1<<b<n)
        b++;
    
    vector<complex<double>> f(n);
    for (int i=0; i<n; i++)
    {
        int t=0;
        for (int j=0; j<b; j++)
            t = t<<1 | i>>j&1;
        f[t] = v[i];
    }

    vector<complex<double>> w(n);
    for (int i=0; i<n; i++)
    {
        double t = (inv?2:-2)*pi*i/n;
        w[i] = complex<double>(cos(t), sin(t));
    }

    for (int i=0; i<b; i++)
    for (int j=0; j<n; j+=1<<-~i)
    for (int k=0; k<1<<i; k++)
    {
        int l = k+(1<<i);
        complex<double> x = f[j+k];
        complex<double> y = f[j+l];
        f[j+k] = x + y*w[k<<(b-i-1)];
        f[j+l] = x + y*w[l<<(b-i-1)];
    }

    if (inv)
        for (int i=0; i<n; i++)
            f[i] /= n;

    return f;
}

//  ret[i] = ΣA[t]B[i-t]
vector<double> convolution(vector<double> A, vector<double> B)
{
    int n = 1;
    while (n<(int)A.size()+(int)B.size()-1)
        n *= 2;
    vector<complex<double>> a(n), b(n);
    for (int i=0; i<(int)A.size(); i++)
        a[i] = A[i];
    for (int i=0; i<(int)B.size(); i++)
        b[i] = B[i];
    vector<complex<double>> at = FFT(a, false);
    vector<complex<double>> bt = FFT(b, false);
    vector<complex<double>> ct(n);
    for (int i=0; i<n; i++)
        ct[i] = at[i]*bt[i];
    vector<complex<double>> c = FFT(ct, true);
    vector<double> C((int)A.size()+(int)B.size());
    for (int i=0; i<(int)C.size(); i++)
        C[i] = c[i].real();
    return C;
}

#include <iostream>

int main()
{
    int N;
    cin>>N;
    vector<long long> A(N);
    for (long long &a: A)
        cin>>a;

    long long P = 200003;
    
    vector<long long> pow(P);
    long long g = 1;
    for (int i=0; i<P; i++)
    {
        pow[i] = g;
        g = g*2%P;
    }
    vector<long long> log(P);
    for (int i=0; i<P; i++)
        log[pow[i]] = i;

    vector<double> C(P);
    for (long long a: A)
        if (a>0)
            C[log[a]] += 1.0;
    vector<double> C2 = convolution(C, C);

    long long ans = 0;
    for (int i=0; i<2*P-1; i++)
        ans += pow[i%(P-1)]*(long long)(C2[i]+.5);
    for (long long a: A)
        ans -= a*a%P;
    ans /= 2;
    cout<<ans<<endl;
}