As you see, the problem is about “Discrete Logarithm”. All the numbers mentioned below are unnegetive intergers.

First of all, define \(\delta_p(x)\) as the minimum positive interger such as \(x^{\delta_p(x)} \equiv 1 \pmod p\),if \(p\) is a prime number. It’s obvious that if \(0 \le i < j < \delta_p(x)\), we can claim that \(x^i \not \equiv x^j \pmod p\). And a famous Femart’s Theorem tells us if \(p \nmid x\), we know \(x^{p-1} \equiv 1 \pmod p\). If \(\delta_p(x) = p-1\), we call \(x\) a Primitive Root of \(p\). There is also a significant property which I won’t prove: if \(x^k \equiv 1 \pmod p\), we know \(\delta_p(x) \mid k\).

If \(t\) is a Primitive Root of \(p\), for all \(1 \le x \le p-1\), there exists a \(k\) such that \(t^k \equiv x \pmod p\), and the \(k\) less than \(p-1\) is unique.

If we know \(t\) and \(x\), to find the \(k\) such that \(t^k = x\), we only need to let \(x = \log_t x\).

In the mod calculations,we can have similar operations. Define \(\text{ind}_t x\) as the minumum positive value such as \(t^{\text{ind}_t x} \equiv x \pmod p\), if \(t\) is a Primitive Root of \(p\).

Note that in most places, \(\text{ind}_t 1 =0\). But here I make it be \(p-1\) to make the problem easy to solve.

So the problem is equal to find pairs \((i,j)\) such as:

• \(1 \le i , j \le n\).
• There is a positive \(k\) such that \(k \times \text{ind}_t a_i \equiv \text{ind}_t a_j \pmod {p-1}\).

The proof is for exercise.

Let \(b_i\) be \(\gcd(\text{ind}_t a_i,p-1)\), the problem is also equal to:

• \(1 \le i , j \le n\).
• \(b_i | b_j\).

And most importantly, \(\dfrac{p-1}{b_i} = \delta_p (a_i)\).

Proof is for exercise.

So we can calculate \(\delta_p(a_i)\) ( I won’t explain. You can search the internet ), then we know \(b_i\).

Notice that \(b_i | p-1\), so there are at most \(d(p-1) \le 1.1 \times 10^4\) different \(b\).

We can solve that in a time complexity of \(O(\sqrt p + n \log^2 p + d^2(p-1))\).

Code :

#include<bits/stdc++.h>
#define int long long
#define ffor(i,a,b) for(int i=(a);i<=(b);i++)
#define roff(i,a,b) for(int i=(a);i>=(b);i--)
using namespace std;
const int MAXN=2e5+10;
int n,p,a[MAXN];
int qpow(__int128 base,int p,int mod) {
	__int128 ans=1;
	while(p) {
		if(p&1) ans=ans*base%mod;
		base=base*base%mod,p>>=1;	
	}
	return ans;
}
vector<int> frac;
map<int,int> mp;
int calc_b(int v) {
	int ans=p-1;
	for(auto id:frac) if(qpow(v,ans/id,p)==1) ans/=id;
	return (p-1)/ans;
}
signed main() {
	ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
	cin>>n>>p; ffor(i,1,n) cin>>a[i];
	int v=p-1;
	ffor(i,2,v/i) while(v%i==0) frac.push_back(i),v/=i;
	if(v) frac.push_back(v);
	sort(frac.begin(),frac.end());
	ffor(i,1,n) mp[calc_b(a[i])]++;
	int ans=0;
	for(auto pr1:mp) for(auto pr2:mp) if(pr1.first%pr2.first==0) ans+=pr1.second*pr2.second;
	cout<<ans;
	return 0;
}

Refference : OI-Wiki

Chinese version:Link