First, for each \(A_i\), find its multiplicative order modulo \(P\). The multiplicative order is the smallest positive integer \(k_i\) such that \(A_i ^ {k_i} \equiv 1 \pmod P\).

To do this, we use the fact that by Fermat’s little theorem, \(A_i ^ {P-1} \equiv 1 \pmod P\), since \(P\) is prime. This means that the powers of \(A_i\) modulo \(P\) forms a cycle, the length of which is one of the divisors of \(P-1\), and which is our desired multiplicative order. To find it out, we could first factorize \(P-1\). Now start with the tentative multiplicative order \(x {\ :=\ } P-1\). Then for each prime factor \(p\) of \(P-1\) (including duplicate factors), we test whether \(A_i ^ {x / p} \equiv 1 \pmod P\); if it is, the cycle length is a divisor of \(x/p\), so we can assign \(x {\ :=\ } x / p\) and continue. The final value of \(x\) is our desired multiplicative order. This step should take \(O(\sqrt P + N \log^2 P)\) time.

Now, an \((i, j)\) pair as in the problem statement is satisfactory iff \(k_j\) divides \(k_i\) (see note * below for a proof sketch of this). To count these pairs, first we can group our found \(k_i\) into (value, count) pairs, then we can look at pairs of (value, count) pairs to count the possible \((i, j)\) pairs.

There might be up to \(\tau(P-1)\) distinct values for multiplicative order, where \(\tau\) is the divisor count function. According to this list, \(\tau(P-1) \leq 10752\) (in fact it’s a bit less as you’re limited to primes for \(P\)), so using \(O(\tau(P-1)^2)\) time to count the satisfactory pairs is fine, given the fairly generous time limit.

It’s also possible to count these pairs in \(O(\tau(P-1) \cdot \omega(P-1)) \approx O(\tau(P-1) \log P / \log \log P)\), where \(\omega(P-1)\) is the number of distinct prime factors of \(P-1\). Specifically, you treat each distinct prime factor as a “dimension” that you can do a kind of fast zeta transform on, to obtain the count of the number of \(k_i\) divisible by each divisor of \(P-1\), via the following pseudocode:

• Let \(cnt\) be a map from each divisor of \(P-1\) to its multiplicity among \(k_i\).
• For each distinct prime factor \(p\) of \(P-1\):
  • For each divisor \(d\) of \(P-1\), from largest to smallest:
    • If \(p\) divides \(d\), add \(cnt[d]\) to \(cnt[d/p]\).
  • (the dimension representing the powers of \(p\) in the lattice of divisors has now been converted from “point value” to “suffix sum value”)

\(cnt[d]\) now counts the number of \(k_i\) that is divisible by \(d\).

* Proof sketch and derivation: take an arbitrary element with multiplicative order exactly \(P-1\); this is known as a primitive root, and call it \(g\). This always exists for prime moduli. The first \(P-1\) powers of \(g\) modulo \(P\) will cover all nonzero elements \([1, P)\), and thus each \(A_i\) will have an index \(0 \le \text{ind}_g(A_i) < P-1\) such that \(g^{\text{ind}_g(A_i)} \equiv A_i\).

Now consider a similar problem, but where all \(A_i\) has been replaced with \(a_i = \text{ind}_g(A_i)\). The problem becomes counting the pairs \((i, j)\) where there is a positive integer \(k\) such that \(k \cdot a_i \equiv a_j \pmod {P-1}\). By Bezout’s identity, these are exactly the pairs where \(\gcd(a_i, P-1)\) divides \(\gcd(a_j, P-1)\).

The multiplicative order of \(A_i\) modulo \(P\) is equivalent to the “additive order” of \(a_i\) modulo \(P-1\), therefore \(\large\frac{P-1}{\gcd(a_i, P-1)}= k_i\), and “\(\gcd(a_i, P-1)\) divides \(\gcd(a_j, P-1)\)” is equivalent to “\(k_j\) divides \(k_i\)”. Finding multiplicative orders is much easier computationally than finding \(\text{ind}_g(A_i)\) for all \(N\) values (which’ll cost \(O(N + \sqrt {NP})\) time and space via baby-step-giant-step), so we use those instead.