#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
struct Complex{
    double r,i;
    Complex(double r,double i):r(r),i(i){}
    Complex(){}
	inline Complex operator + ( const Complex& rhs )  const  {
		return Complex ( r + rhs.r, i + rhs.i ) ;
	}
	inline Complex operator - ( const Complex& rhs )  const  {
		return Complex ( r - rhs.r, i - rhs.i ) ;
	}
	inline Complex operator * ( const Complex& rhs )  const  {
		return Complex ( r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r ) ;
	}
	inline void operator /= ( const double& x )   {
		r /= x, i /= x ;
	}
	inline void operator *= ( const Complex& rhs )   {
		*this = Complex ( r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r ) ;
	}
	inline void operator += ( const Complex& rhs )   {
		r += rhs.r, i += rhs.i ;
	}
};

const int maxn=4e6+6;
#define PI 3.14159265354
int pos[maxn];
void init ( const int& n )  {
    for ( int i = 0 ,j=0; i < n ; ++ i )  {
        pos[i]=j;for (int l = n >> 1; (j ^= l) < l; l >>= 1);
    }
}

void transform ( Complex *a, const int& n, bool inverse ) {
    for ( int i = 0; i < n ; ++ i ){
        if (i > pos[i]) std ::swap(a[i], a[pos[i]]);
    }

    for ( int l = 2 ; l <= n ; l <<= 1 )  {
        int m = l / 2;Complex omega={cos(2*PI/l),inverse?-sin(2*PI/l):sin(2*PI/l)};
        for ( Complex *p = a ; p != a + n ; p += l )  {
            Complex x={1,0};
            for ( int i = 0 ; i < m ; ++ i ,x *= omega)  {
                Complex t = x  * p [m + i] ;
                p [m + i] = p [i] - t ;p [i] += t ;
            }
        }
    }
}

void dft ( Complex *a, const int& n )  {
    transform ( a, n, 0 ) ;
}

void idft ( Complex *a, const int& n )  {
    transform ( a, n, 1) ;
    for ( int i = 0 ; i < n ; ++ i ) a [i] /= n ;
}
Complex A[maxn],B[maxn],C[maxn];
void FFT(int n,int m) {// len(A),len(B)
    int cnt=1;while(cnt<=(n+m)) cnt<<=1;
    init(cnt);dft(A,cnt);dft(B,cnt);
    for(int i=0;i<cnt;i++) C[i]=A[i]*B[i];
    idft(C,cnt);
    for(int i=0;i<=n+m;i++) C[i].r=ll(C[i].r+0.3);
}
const int mod=200003;
const int g=2;
std::vector<ll> pw(mod),rev(mod);
int main(){
	ll now=1;pw[0]=1;rev[1]=0;
	for(int i=1;i<mod-1;i++) {
		now=now*g%mod;pw[i]=now;
		rev[now]=i;
	}
	int n;cin>>n;
	ll ans=0;
	for(int i=1;i<=n;i++) {
		ll tem;cin>>tem;
		if(tem==0) continue;
		ans-=tem*tem%mod;
		A[rev[tem]].r+=1;
		B[rev[tem]].r+=1;
	}
	FFT(mod,mod);	
	for(int i=0;i<=2*mod;i++) {
		ll x=C[i].r;
		ans+=pw[i%(mod-1)]*x;
	}
	cout<<ans/2<<endl;
}
### Prompt

Your challenge is to write 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>
using namespace std;
typedef long long ll;
struct Complex{
    double r,i;
    Complex(double r,double i):r(r),i(i){}
    Complex(){}
	inline Complex operator + ( const Complex& rhs )  const  {
		return Complex ( r + rhs.r, i + rhs.i ) ;
	}
	inline Complex operator - ( const Complex& rhs )  const  {
		return Complex ( r - rhs.r, i - rhs.i ) ;
	}
	inline Complex operator * ( const Complex& rhs )  const  {
		return Complex ( r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r ) ;
	}
	inline void operator /= ( const double& x )   {
		r /= x, i /= x ;
	}
	inline void operator *= ( const Complex& rhs )   {
		*this = Complex ( r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r ) ;
	}
	inline void operator += ( const Complex& rhs )   {
		r += rhs.r, i += rhs.i ;
	}
};

const int maxn=4e6+6;
#define PI 3.14159265354
int pos[maxn];
void init ( const int& n )  {
    for ( int i = 0 ,j=0; i < n ; ++ i )  {
        pos[i]=j;for (int l = n >> 1; (j ^= l) < l; l >>= 1);
    }
}

void transform ( Complex *a, const int& n, bool inverse ) {
    for ( int i = 0; i < n ; ++ i ){
        if (i > pos[i]) std ::swap(a[i], a[pos[i]]);
    }

    for ( int l = 2 ; l <= n ; l <<= 1 )  {
        int m = l / 2;Complex omega={cos(2*PI/l),inverse?-sin(2*PI/l):sin(2*PI/l)};
        for ( Complex *p = a ; p != a + n ; p += l )  {
            Complex x={1,0};
            for ( int i = 0 ; i < m ; ++ i ,x *= omega)  {
                Complex t = x  * p [m + i] ;
                p [m + i] = p [i] - t ;p [i] += t ;
            }
        }
    }
}

void dft ( Complex *a, const int& n )  {
    transform ( a, n, 0 ) ;
}

void idft ( Complex *a, const int& n )  {
    transform ( a, n, 1) ;
    for ( int i = 0 ; i < n ; ++ i ) a [i] /= n ;
}
Complex A[maxn],B[maxn],C[maxn];
void FFT(int n,int m) {// len(A),len(B)
    int cnt=1;while(cnt<=(n+m)) cnt<<=1;
    init(cnt);dft(A,cnt);dft(B,cnt);
    for(int i=0;i<cnt;i++) C[i]=A[i]*B[i];
    idft(C,cnt);
    for(int i=0;i<=n+m;i++) C[i].r=ll(C[i].r+0.3);
}
const int mod=200003;
const int g=2;
std::vector<ll> pw(mod),rev(mod);
int main(){
	ll now=1;pw[0]=1;rev[1]=0;
	for(int i=1;i<mod-1;i++) {
		now=now*g%mod;pw[i]=now;
		rev[now]=i;
	}
	int n;cin>>n;
	ll ans=0;
	for(int i=1;i<=n;i++) {
		ll tem;cin>>tem;
		if(tem==0) continue;
		ans-=tem*tem%mod;
		A[rev[tem]].r+=1;
		B[rev[tem]].r+=1;
	}
	FFT(mod,mod);	
	for(int i=0;i<=2*mod;i++) {
		ll x=C[i].r;
		ans+=pw[i%(mod-1)]*x;
	}
	cout<<ans/2<<endl;
}
```