题解 [AGC047C] Product Modulo

显然不能暴力算两两的乘积，而积取模而结果不取模提示我们模数肯定有用。

所有为 \(0\) 的 \(a_i\) 对答案不会产生任何贡献，可以直接删除，下文不再考虑这种情况。同时我们约定，下文的运算（除特殊说明外）均为模 \(P\) 意义下的结果。

设 \(g\) 为 \(P\) 的原根，则每个 \(a_i\) 都可以表示为 \(g^{c_i}\) 的形式，而 \(a_ia_j=g^{c_i+c_j}\)。因此设 \(a\) 中等于 \(g^{i}\) 的数有 \(A_i\) 个，则乘积为 \(g^k\) 的有序数对有 \(\sum\limits_{i+j=k}A_iA_j\) 个，是卷积形式，用 FFT 优化之。

注意根据题意，答案需要去掉所有 \(a_i^2\) 的贡献。

时间复杂度为 \(\Theta(P\log P)\)。

我的代码中 \(g=114513\)，当然你也可以用 \(2\) 或者其他原根。

// Problem: C - Product Modulo
// Contest: AtCoder - AtCoder Grand Contest 047
// URL: https://atcoder.jp/contests/agc047/tasks/agc047_c
// Memory Limit: 1024 MB
// Time Limit: 2000 ms
// 
// Powered by CP Editor (https://cpeditor.org)

//By: OIer rui_er
#include <bits/stdc++.h>
#define rep(x,y,z) for(ll x=(y);x<=(z);x++)
#define per(x,y,z) for(ll x=(y);x>=(z);x--)
#define debug(format...) fprintf(stderr, format)
#define fileIO(s) do{freopen(s".in","r",stdin);freopen(s".out","w",stdout);}while(false)
using namespace std;
typedef long long ll;
typedef complex<double> cp;
const ll N = 2e6+5, P = 200003, g = 114513;
const double pi = acos(-1);

ll n, a[N], pw[N], pos[N], cnt[N], rev[N], ans;
cp A[N];
template<typename T> void chkmin(T& x, T y) {if(x > y) x = y;}
template<typename T> void chkmax(T& x, T y) {if(x < y) x = y;}
void change(cp* a, ll n) {
	rep(i, 0, n-1) {
		rev[i] = rev[i >> 1] >> 1;
		if(i & 1) rev[i] |= n >> 1;
	}
	rep(i, 0, n-1) if(i < rev[i]) swap(a[i], a[rev[i]]);
}
void FFT(cp* a, ll n, ll sgn) {
	change(a, n);
	for(ll l = 2, o = 1; l <= n; l <<= 1, o <<= 1) {
		cp w(cos(2.0*pi/l), sgn*sin(2.0*pi/l));
		for(ll i = 0; i < n; i += l) {
			cp now = 1;
			rep(j, i, i+o-1) {
				cp u = a[j], v = a[j + o];
				a[j] = u + v * now;
				a[j + o] = u - v * now;
				now *= w;
			}
		}
	}
	if(sgn == -1) rep(i, 0, n-1) a[i] /= n;
}

int main() {
	pw[0] = 1;
	pos[1] = 0;
	rep(i, 1, P-2) {
		pw[i] = g * pw[i-1] % P;
		pos[pw[i]] = i;
	}
	scanf("%lld", &n);
	rep(i, 1, n) scanf("%lld", &a[i]);
	rep(i, 1, n) if(a[i]) ++cnt[pos[a[i]]];
	rep(i, 0, P-2) A[i] = cnt[i];
	ll m = 1;
	for(; m <= P + P; m <<= 1);
	FFT(A, m, 1);
	rep(i, 0, m-1) A[i] *= A[i];
	FFT(A, m, -1);
	rep(i, 0, 2*(P-2)) ans += (ll)(A[i].real() + .5) * pw[i%(P-1)];
	rep(i, 0, P-2) ans -= cnt[i] * pw[2*i%(P-1)];
	printf("%lld\n", ans / 2);
	return 0;
}