Interesting stuff. I've seen that thing before. I think it is called Freddy's factoring trick or FFT for short. What is the probability that competitive programming has two different FFTs?

By the way, since you this is your own problem, why don't you lower the constraints to $$$10^3$$$ and just do two for loops?

For the specific case $$$K = 1$$$, we can solve it as follows: For each $$$i$$$: if $$$\gcd(a_i, m) > 1$$$, we cannot find a corresponding $$$j$$$, so just continue. If $$$\gcd(a_i, m) = 1$$$, then let $$$m = p_1^{\alpha_1}p_2^{\alpha_2} \cdots p_k^{\alpha_k}$$$, and we get $$$k$$$ congruences, where the $$$\ell^{\text{th}}$$$ congruence gives us the value of the corresponding $$$a_j$$$ modulo $$$p_{\ell}^{\alpha_{\ell}}$$$. We can solve these congruences using CRT to get the corresponding $$$a_j$$$ modulo $$$m$$$.

The problem you proposed is more difficult than the specific case $$$K = 1$$$; In this case, both $$$arr[i]$$$ and $$$arr[j]$$$ must be coprime to $$$M$$$, so the same solution for prime $$$M$$$ works. Anyway, here's the general solution.

Let's consider all values modulo $$$M$$$. We can do some algebra to get a solution similar to the case of $$$M$$$ being prime. Algebra is not intuitive to me so here's the intuition first:

Fix $$$a = arr[i]$$$. If $$$a$$$ and $$$M$$$ are coprime, then $$$b = arr[j]$$$ is fixed modulo $$$M$$$; it is $$$Ka^{-1}$$$. Now suppose $$$M$$$ is even and $$$a = 2$$$. Then if $$$K$$$ is odd there is no solution — but otherwise $$$b$$$ has two solutions modulo $$$M$$$, since the solutions repeat cyclicly with period $$$\frac{M}{2}$$$. This would also be true if $$$gcd(a, M) = 2$$$.

Now let's do it formally with some algebra. Observe that if $$$x \equiv y \pmod z$$$, and all are divisible by some $$$c$$$, then it is true iff $$$\frac{x}{c} \equiv \frac{y}{c} \pmod{\frac{z}{c}}$$$.

Again, if we fix $$$a = arr[i]$$$, denote $$$g = gcd(a, M)$$$, then we are looking for $$$b = arr[j]$$$ such that $$$ab \equiv K \pmod M \iff \frac{a}{g}b \equiv \frac{K}{g} \pmod{\frac{M}{g}}$$$. Denote $$$a' = \frac{a}{g}$$$, same for others; $$$a'b \equiv K' \pmod {M'}$$$. Now $$$a', M'$$$ are coprime, so it is true iff $$$b \equiv K'a'^{-1} \pmod {M'}$$$. Note that if $$$K$$$ isn't divisible by $$$g$$$ there is no solution.

If we compute $$$c_i$$$ as the number of values in the array equal to $$$i$$$ mod $$$M$$$, then for each $$$a$$$ we can compute the answer in $$$O(g)$$$ by summing up the subset of relevant $$$b$$$'s. If you compute these subset sums lazily (only when needed) and cache their results, you can show an $$$O(N + M \sqrt{M})$$$ time complexity upperbound.

We can also do better; denote $$$dp_{i, M'} = \sum_{k = 0}^{g - 1} c_{i + kM'}$$$, where $$$g = \frac{M}{M'}$$$ as above. These are exactly the subset sums you want to compute. We can compute this in a dp manner, in decreasing order of $$$M'$$$. Let $$$p$$$ be the smallest prime divisor of $$$g$$$, then $$$dp_{i, M'} = \sum_{k = 0}^{p - 1} dp_{i + kM', pM'}$$$.

If we denote $$$pr_x$$$ as the smallest prime factor of $$$x$$$, then the time complexity is $$$\sum_{g = 1, g | M} \frac{M}{g} \cdot pr_g$$$. It is at least $$$\Omega(M \log \log M)$$$, due to sum of divisors asymptotics. For an upperbound, note that the inner summand is $$$\frac{M}{g/pr_g}$$$. The term $$$g/pr_g$$$ is a divisor of $$$M$$$, and in addition it can only appear a number of times bounded by the number of distinct prime divisors of $$$M$$$. If $$$M$$$ has $$$k$$$ distinct prime divisors, then the following is an upperbound: $$$k \cdot \sum_{g = 1, g | M} \frac{M}{g}$$$. Since $$$k = O(\frac{\log M}{\log \log M})$$$, we get an upperbound of $$$O(M \log M)$$$.