First refer to the official editorial.

We consider the subproblem of handling a prime \(p \equiv 1 \pmod 4\). Let \(p = \pi \bar\pi\) be the factorization of \(p\) into a Gaussian prime and its conjugate. We want to track the distribution of \(\pi ^k\bar\pi ^{e - k} \text{ mod } C\) for integers \(0 \le k \le e\).

Consider the sequence \((\pi ^0, \pi ^1, \pi ^2, \dots) \text{ mod } C\). This is an eventually-periodic sequence, or what I myself call a rhotic sequence, out of analogy to the shape of the Greek letter \(\rho\) (rho). It’s a sequence \((a_0, a_1, a_2, \dots)\) where:

• there exists an index \(\mu\) where \((a_\mu, a_{\mu+1}, a_{\mu+2}, \dots)\) is periodic; we take the minimum of such indices,
• there exists a period \(\lambda\) where \(a_{j + \lambda} = a_j\) for all \(j \ge \mu\), we also take the minimum possible period.

\((a_0, a_1, \dots, a_{\mu - 1})\) is called the preorbit and corresponds to the “tail” of \(\rho\), while \((a_\mu, a_{\mu + 1}, \dots, a_{\mu + \lambda - 1})\) is called the orbit and corresponds to the loop of \(\rho\).

Since \(\mu + \lambda \le C^2\) in this case we could easily track the sequence and find \(\mu\) and \(\lambda\), and then define a function to find arbitrary \(\pi^k \text{ mod } C\).

Now consider the sequence \((\bar\pi ^0, \bar\pi ^1, \bar\pi ^2, \dots) \text{ mod } C\). Since \(\bar z ^k = \overline {z^k}\) for any complex number \(z\), these are just the conjugates of the first sequence, and thus is also a rhotic sequence with the same \(\mu\) and \(\lambda\).

We now consider the sequence \(s_k = \pi ^k\bar\pi ^{e - k} \text{ mod } C\) for \(0 \le k \le e\).

• If \(e < 2 \mu + \lambda\), it’s easiest to just track the entire sequence with brute force.
• Otherwise, for \(0 \le k < \mu\) and \(e - \mu < k \leq e\), one of the power sequences is in its preorbit, so we should track these by brute force. (Note that the values in the latter range are the conjugates of the values in the former range.)
• For \(\mu \le k \le e - \mu\), both power sequences are in their orbits, therefore this segment of \(s_k\) is \(\lambda\)-periodic. Thus its distribution can be tracked in \(O(\lambda)\) by obtaining the quotient and remainder of dividing the segment length \((e + 1 - 2 \mu)\) by \(\lambda\), then distributing among the first \(\lambda\) terms of the segment.