# Problem Link:

# Difficulty:

Easy-Medium

# Pre-requisites:

Number Theory: Principle of Inclusion-Exclusion, Moebius Function, Euler’s Totient Function, Modular Arithmetic

# Problem:

You are given K intervals [A_{0}, B_{0}], [A_{1}, B_{1}], …, [A_{K-1}, B_{K-1}]. From each interval i, you randomly pick a number X_{i}. What is the expected value of the GCD of the picked numbers (X_{0}, X_{1}, …, X_{K-1})? If the answer is expressed as a fraction P/Q, then output an integer N in [0, MOD-1] (MOD = 10^9+7) such that P + Q*N is divisible by MOD.

# Preliminary Definitions / Lemmas etc:

You may skip whichever you’re familiar with.

Euler’s phi function. Phi(N) = |{n: n <= N and n co-prime to N, i.e. gcd(n, N) = 1}|

Formula for phi(N):

if N = ∏(p_{i}^{ai}) (after prime factorization), then

phi(N) = ∏ p_{i}^{ai-1}*(p_{i}-1)

Lemma: ∑_{d|N} (phi(d)) = N for all N.

A simple proof of the above:

Consider N fractions : 1/N, 2/N, …, N/N.

Rewrite these in reduced form: 1/N, …, 1/1.

Now, for each d|N, there will be phi(d) of the above fractions having denominator = d.

Hence the lemma follows, by counting in two ways, the same N fractions.

Moebius function mu(N):

mu(N) = 0 if p^{2} | N for some prime p (i.e. if N is not “square-free”)

else (-1)^(#primes dividing N) otherwise.

Hence, both phi and mu can be calculated efficiently using the Sieve of Eratosthenes.

Lemma: Moebius inversion formula:

∑ _{ d|N } (d * mu(N/d)) = phi(N).

Lemma: gcd(x1, x2, x3, …, xn) = g <=> gcd(x1/g, x2/g, x3/g, …, xn/g) = 1.

This can be seen by taking prime factorizations of xi’s and comparing LHS and RHS of both equations.

Inverse modulo prime:

Given a prime number ‘M’, to every number x in the range [1, M-1], you have a unique number y, denoted by x^{-1}, such that x * y = 1 (mod M).

Given prime M, and value x, x^{-1} can be efficiently found using either of the following methods:

- Euclidean algorithm: given x and M, find y and n such that x * y + M * n = 1. Now, by taking modulo M on both sides, we get that x
^{-1}= y (mod M) - Fermat’s little theorem: for any x not divisible by prime M, x
^{M-1}= 1 (mod M). This shows that x^{M-2}= x^{-1}(mod M). Further x^{M-2}is calculated efficiently by repeated squaring in log(M) time.

# Quick Explanation:

We try to find P = sum of gcd(x_{0}, x_{1}, …, x_{K-1}) over all possible values of x_{i}.

Also find Q = number of tuples (x_{0}, x_{1}, …, x_{K-1}) = product(B_{i} - (A_{i} - 1)).

With such P and Q, output N = -P * Q^{-1} (mod 10^9+7).

The question boils down to finding P. Using principle of inclusion-exclusion (and concept of moebius function), get that P = ∑_{R, m, mR <= min(Bi)}(R * mu(m) * ∏_{0<=i < K }(floor(B/(m * R)) - floor((A*-1)/(m * R))

This is equivalent (using moebius inversion formula) to P = ∑_{Y <= min(B_i)} (phi(Y) * ∏_{0<=i < K }(floor(B*/Y) - floor((A*-1)/Y)).

# Detailed Explanation:

Consider a K-dimensional array (lets call it GCD) where the value in the (x_{0}, x_{1}, x_{2}, …, x_{K-1})'th cell is gcd(x_{0}, x_{1}, x_{2}, …, x_{K-1}). We are interested in finding the expected value of GCD[x_{0}][x_{1}]…[x_{K-1}] over all possible values of x_{0}, x_{1}, …, x_{K-1}. Since all x_{i}’s are picked randomly, the expected value is equal to the arithmetic mean. Finally, this question is equivalent to finding the sum of the values in this K-dimensional GCD array, since the mean is then just dividing this sum by the number of elements.

In fact, looking at P/Q, we can get P = sum of elements, and Q = number of elements. Let MOD = 10^9 + 7. We need to find N such that P + Q * N is divisible by MOD. In other words, P + Q * N = 0 (mod MOD). Since MOD is prime, we can find “Q^{-1}” modulo MOD, which gives us that N = -P * Q^{-1} (mod MOD). Hence, finding P = sum of elements in GCD array modulo MOD, and Q = number of elements = ∏ (B_{i} - (A_{i}-1)) (modulo MOD), we would have enough information to get N. Here, Q can be calculated in O(K) time per test-case.

P is calculated in the following manner.

Firstly, let us try to solve the following intermediate question:

Given any arbitrary set of K intervals (not necessarily the A_{i}, B_{i} as desired in the question), denoted by (U_{i}, V_{i}] (intervals exclusive of U_{i}, inclusive of V_{i}), let us ask: How many points are there in the corresponding K-dim space : i’th coordinate lies in (U_{i}, V_{i}], whose gcd(point’s coordinates) is 1. Let us call this set of points as C(U, V). The answer to our intermediate question is now |C(U, V)|, where U and V specify the intervals’ endpoints we are interested in.

This helps us because

the number of points in (A_{i} - 1, B_{i}] intervals whose coordinates’ gcd is R (for arbitrary R) =

the number of points in (floor((A_{i}-1)/R), floor(B_{i}/R)] whose coordinates’ gcd is 1.

The above follows from the fact that gcd(x_{0}, x_{1}, x_{2}, …, x_{K-1}) = R <=> gcd(x_{0}/R, x_{1}/R, x_{2}/R, …, x_{K-1}/R) = 1. And the fact that

A_{i} - 1 < x_{i} <= B_{i} <=> floor((A_{i} - 1)/R) < x_{i}/R <= floor(B_{i}/R)

Thus, P = required sum = ∑ _{R > 0} (R * |C(floor((A_{i}-1)/R), floor(B_{i}/R))|)

We now proceed to calculate |C(U, V)|. Qualitatively speaking, from the universal set of points ARR(U, V), we remove those whose coordinates’ gcd > 1, to get C(U, V).

Let ARR(U, V) = {(x_{0}, x_{1}, …, x_{K-1}) : U* < x_i <= V*}.

Let G(U, V, D) = {(x_{0}, x_{1}, …, x_{K-1}) : D | gcd(x_{0}, x_{1}, …, x_{K-1}), U* < xi <= V*}

From the above 2 definitions, it is easily seen that

C(U, V) = ARR(U, V) \ UNION_{D > 1} G(U, V, D)

: from the set of all possible points, remove points whose gcd is divisible by any D > 1.

In fact, we need consider only D as being prime.

Indeed, UNION _{ D > 1} G(U, V, D) = UNION _{prime D} G(U, V, D) since if the gcd of a point is D (> 1), then some prime p|D and that point belongs to the set G(U, V, p). Similarly if p|gcd(point’s coordinates), then that point belongs to the set G(U, V, p). Also, we only consider prime D <= min(V) since for prime D > min(V), G(U, V, D) is ∅.

Now, using the fact that if co-prime m and n are such that m|x and n|x, it is equivalent to (m * n)|x, we get that

given co-prime D1, and D2: G(U, V, D1) ∩ G(U, V, D2) = G(U, V, D1*D2).

The good thing about the set G(U, V, D) is that it is easy to count its size. In fact, |G(U, V, D)| = ∏ (floor(V*/D) - floor(U*/D)).

This is because (x_{0}, x_{1}, …, x_{K-1}) ∈ G(U, V, D)

<=> D|x_{i}, and U_{i} < x_{i} <= V_{i}.

|{x : (D|x and x <= n)}| = floor(n/D), which when applied above gives

|{x_{i} : D|x_{i} and U_{i} < x_{i} <= V_{i}}| = |{x_{i} : D|x_{i} and x_{i} <= V_{i}} \ {x_{i}: D|x_{i} and x_{i} <= U_{i}}| = |{x_{i} : D|x_{i} and x_{i} <= V_{i}}| - |{x_{i}: D|x_{i} and x_{i} <= U_{i}}| = floor(V_{i}/D) - floor (U_{i}/D).

Since x_{i}’s are chosen independently of other x_{i}’s, |G(U, V, D)| = ∏(floor(V_{i}/D) - floor(U_{i}/D)).

And now we are well equipped to apply Principle of Inclusion-Exclusion which gives us:

|C(U, V)| = ∑_{M = set of distinct primes} (-1)^|M| * |G(U, V, ∏ _{p ∈ M} p)|)

= ∑_{M = set of distinct primes} (-1)^|M| * (floor(V_{i}/∏ _{p ∈ M} p) - floor(U_{i}/∏ _{p ∈ M} p))

This gives us our first possible algorithm.

P = 0

For R from 1 to N(= min(B_{0}, B_{1}, …, B_{K-1})

V_{i} = B_{i}/R, U_{i} = (A_{i}-1)/R

for all m = product of distinct primes s.t. m <= min(floor(N/R))

P += R * (-1)^(number of distinct primes in m) * ∏(floor(V_{i}/m) - floor(U_{i}/m))

The above algorithm takes time O(∑(floor(N/j)) * K) per test-case, which is O(N * logN * K) per testcase. This is too slow however. It can be sped up by noticing that product(floor(V_{i}/m) - floor(U_{i}/m)) = product (floor(B_{i}/(m * R)) - floor((A_{i}-1)/(m * R))). This product can be calculated once for all m * R product pairs in O(NK) time. This now gives us an algorithm which is O(N * logN + N * K) per test-case. This is fast enough to pass.

Also, note that in the above, (-1)^(number of distinct primes in m) = mu(m) where m = product of distinct primes. Also, since mu(m) = 0 if m is not the product distinct primes, the above can also be written as

P = ∑ _{m, R s.t. m * R <= N} (R * mu(m) * ∏(floor(B_{i}/(m * R)) - floor((A_{i}-1)/(m * R)))

mu(m) can be precomputed for all m <= 2*10^5 using Sieve.

Also, if you rewrite m * R = Y, then you can write ∑(R * mu(Y/R)) as coeff(Y), giving P = sum ( coeff(Y) * ∏(floor(B_{i}]/Y) - floor((A_{i}-1)/Y))).

Using the moebius inversion formula, you get that coeff(Y) is actually phi(Y). So using sieve, we can precompute in time O(N * log log N) the values of phi(i) for all i <= 2*10^5. This gives us a O(N * log log N + T * N * K) time algorithm.

## Alternate approach:

This approach gives us directly P = ∑ ( phi(Y) * ∏(floor(B_{i}/Y) - floor((A_{i}-1)/Y))) using some clever symbolic manipulation.

Here, we use that n = ∑ _{d|n} (phi(d)).

sum(over A_{i} <= x_{i} <= B_{i}) GCD(x_{0}, x_{1}, x_{2}, …, x_{K-1})

= sum(over A_{i} <= x_{i} <= B_{i}) of (sum(over d | GCD(x_{0}, x_{1}, x_{2}, …, x_{K-1})) phi(d)

= sum(over A_{i} <= x_{i} <= B_{i}) of (sum(over d | all of x_{0}, x_{1}, x_{2}, …, x_{K-1})) phi(d)

= sum(over d <= min(B_{i})) of(sum (over d | x_{0}, and d | x_{1}, and d | x_{2}, …, and d | x_{K-1} and A_{i} <= x_{i} <= B_{i})) phi(d)

= sum(over d <= min(B_{i})) of ∏(floor (B_{i}/d) - floor((A_{i}-1)/d)) * phi(d)

# Author’s Solution:

Can be found here

# Tester’s Solution:

Can be found here