GCDARR2 - Editorial

PROBLEM LINK:

Practice
Contest: Division 1
Contest: Division 2
Contest: Division 3
Contest: Division 4

Author: raysh07
Tester: iceknight1093
Editorialist: iceknight1093

DIFFICULTY:

Easy

PREREQUISITES:

Dynamic programming, inclusion-exclusion

PROBLEM:

Given N and M, count the number of integer arrays A of length N with elements in [1, M] such that:

  • \gcd(A_1, \ldots, A_N) = 1, and
  • \gcd(A_L, \ldots, A_R) \gt 1 for any (L, R) \ne (1, N)

In this version, N \le 2\cdot 10^5 and \sum M \le 2\cdot 10^5.

EXPLANATION:

From the easy version, we already know that the solution is as follows:

  • First compute dp[g] as the number of sequences of length N-2 with elements in [1, M] and GCD exactly g.
  • Then, for each triple (g, A_1, A_N) such that \gcd(A_1, g) \gt 1, \gcd(A_N, g) \gt 1, \gcd(g, A_1, A_N) = 1, add dp[g] to the answer.

Let’s work on optimizing each part separately.


First, consider the computation of dp[g].
Recall that the transition was

dp[g] = (c_g)^{N-2} - \sum_{k\ge 2} dp[k\cdot g]

where c_g = \left\lfloor\frac{M}{g}\right\rfloor is the number of multiples of g in [1, M].

This DP seems like it’s \mathcal{O}(M^2) when run for every g, but it’s actually not!
Note that we do \mathcal{O}(\frac{M}{g}) work for a fixed g, so the total work done is

M + \frac{M}{2} + \frac{M}{3} + \ldots + \frac{M}{M}

which is \mathcal{O}(M\log M) overall, because it’s the Harmonic sum.

So, this part is \mathcal{O}(M\log N + M\log M) in total, easily fast enough for us.


Next, let’s look at actually computing the answer.

If we fix g, we’re interested in the number of triples (g, A_1, A_N) such that \gcd(A_1, g) \gt 1, \gcd(A_N, g) \gt 1, \gcd(g, A_1, A_N) = 1.
Knowing this count, we can multiply it by dp[g] and add it to the answer.

Let’s focus on the first condition alone: \gcd(A_1, g) \gt 1.
How many choices of A_1 exist?

One way to answer this is to note that \gcd(A_1, g) \gt 1 if and only if A_1 and g share a prime factor.
So, let p_1, p_2, \ldots, p_k be the distinct prime factors of g.
A_1 must then contain at least one of these p_i as a factor; equivalently A_1 must be a multiple of at least one of them.

Counting the number of valid choices of A_1 can now be done using inclusion-exclusion.
Specifically,

  • There are \left\lfloor \frac{M}{p_i} \right\rfloor multiples of p_i in [1, M].
    We start by adding up all these values.
  • However, this overcounts: for example a multiple of p_1 and p_2 would’ve been counted twice.
    So, for each pair of primes, we subtract out the number of multiples of their product.
  • This now undercounts: a multiple of three primes would’ve been added three times but also subtracted three times.
    So, for each triplet of primes, we add back in the number of multiples of their product.
    \vdots
  • In general, for each subset of ${p_1, p_2, \ldots, p_k},
    • If the subset has odd size, add the number of multiples of its product.
    • If the subset has even size, subtract the number of multiples of its product.

Note that each product-of-primes will still be a divisor of g, so at best we go through all the distinct divisors of g here.
Thus, by the Harmonic summation, again the complexity of doing this for all g is bounded by \mathcal{O}(M \log M).


One way of formalizing the above argument better, and not having to deal explicitly with prime masks and whatnot, is to look at the Mobius function.

The Mobius function \mu(n) is defined as:

\mu(n) = \begin{cases} 1 & \text{if } n = 1 \\ (-1)^k & \text{if } n \text{ is the product of } k \text{ distinct primes} \\ 0 & \text{otherwise} \end{cases}

Note that this is exactly what we’re looking for: products of primes get +1/-1 signs based on parity of number of terms, and everything else is essentially ignored.
The only difference is that we invert the signs’ parity (we want +1 for odd), which is easy to handle by just multiplying -1.

So, for a fixed g, the value we’re computing can be written as just

- \sum_{\substack{d\mid g \\ d\gt 1}} \mu(d) \cdot\left\lfloor \frac{M}{d} \right\rfloor

Since we’re still just going through all divisors of g, the complexity remains \mathcal{O}(M\log M) across all g.
Writing it this way saves us from having to prime factorize g and then go through all masks of its divisors - the mobius function does that automatically as we go through divisors.

Let the above quantity be denoted S.
We now know that there are S choices for A_1.
Similarly, there are S choices for A_N since it has the exact same constraints.
Thus, there are S^2 choices in total.

However, not all of these choices are good - we didn’t account for the condition \gcd(g, A_1, A_N) = 1.

Luckily, this is easily fixable - we can simply subtract out cases where \gcd(g, A_1, A_N) \gt 1 holds.
This can, once again, be done using the Mobius function and going over divisors of g, because \gcd(g, A_1, A_N) \gt 1 if and only if g, A_1, A_N all share a common prime factor.
By the same inclusion-exclusion + Mobius logic as above, this turns into just

- \sum_{\substack{d\mid g \\ d\gt 1}} \mu(d) \cdot\left\lfloor \frac{M}{d} \right\rfloor^2

The complexity of this part hence remains \mathcal{O}(M\log M) and we’re done.

TIME COMPLEXITY:

\mathcal{O}(M\log M + M\log N) per testcase.

CODE:

Editorialist's code (PyPy3)
import math
mod = 998244353

M = 200005
mu = [1]*M
divs = [ [1] for _ in range(M) ]
for i in range(2, M):
    for j in range(i, M, i):
        divs[j].append(i)
    
    if len(divs[i]) == 2:
        for j in range(i, M, i):
            mu[j] *= -1
        for j in range(i*i, M, i*i):
            mu[j] = 0

for _ in range(int(input())):
    n, m = map(int, input().split())
    
    dp = [0]*(m+1)
    for x in reversed(range(1, m+1)):
        dp[x] = pow(m//x, n-2, mod)
        for y in range(2*x, m+1, x):
            dp[x] -= dp[y]
        dp[x] %= mod
        
    ans = 0
    for x in range(1, m+1):
        ct = 0
        for d in divs[x]:
            ct += mu[d] * (m // d)
        ct = m - ct
        
        # ct^2 ways to choose a[1], a[n] at base
        # subtract ways where gcd(a[1], a[n], x) > 1
        sub = 0
        for d in divs[x]:
            sub += mu[d] * (m // d) * (m // d)
        
        ans += dp[x] * (ct*ct - (m*m - sub))
    print(ans % mod)