TUPLES - Editorial

[mgch]

Practice

Contest

Author: Kevin Charles Atienza

Tester: Istvan Nagy

Editorialist: Misha Chorniy

Difficulty:

Hard

Pre-Requisites:

fft, number theory, inclusion-exclusion

Problem Statement

You are given array X consisting N elements. All the elements in the input are integers. 2-dimensional array Y is defined in next way, Y_{i,j} = X_{i} * X_{j} mod 359999. Find number of integer tuples a, b, c, d, e, f, where gcd(Y_{a,b}, Y_{c,d}, Y_{e,f}) = 1 modulo 10^9+7. Denote that gcd(0, 0) = 0.

Explanation

The first part of the solution is using not 2-dimensional array Y, instead of that calculate array Z where Z_{i} denotes the frequency of element i in array Y, i will be in a range between 0 and 359998. After that problem looks in a next way: we need to find the sum of all Z_{i} * Z_{j} * Z_{k}, where gcd(i, j, k) = 1.

Subtask 1

N is less or equal than 1000. Let’s create array Y in a naive way in O(N^{2}). Also, simply calculate array Z in O(N^{2}).

Subtask 2

A problem of the solution for subtask 1 in a too slow calculation of array Y, how to improve it? Let’s consider number 359999 in more detail, 359999 = 599 * 601, both numbers 599 and 601 are primes. What can this fact give to us? Key observation: let we have a prime number P, from number theory follows next fact, any prime P has primitive root G, in other words:
For every x = 1..P-1, exists y, such G^y = x, 0 <= y < (P-1).

Assume that we have the similar problem but with modulo P(not 359999), where P is a prime number. Let’s find the primitve root for P, call it G. Problem can be divided into two subproblems: where X_{i} is zero or X_{i} isn’t zero. Because zero can’t be represented using the primitive root of number P. f(X) denoting such number that G^{f(X)} = X and 0 <= f(X) < (P-1)

Y_{i,j} = X_{i} * X_{j} = G^{f(X_{i})}*G^{f(X_{j})} = G^{(f(X_{i}) + f(X_{j}))}, in case if X_{i} != 0 and X_{j} != 0.

There are 4 subproblems:

$0 * 0 = 0$

$0 * B = 0$, $B mod P > 0$

$A * 0 = 0$, $A mod P > 0$

$A * B = G^{f(A)} * G^{f(B)} = G^{f(A)+f(B)}$, $A mod P > 0, B mod P > 0$

	zeroes = 0
	for i = 1..N
		if X[i] % P == 0
			zeroes += 1
	Z[0] = zeroes * zeroes + 
		zeroes * (N - zeroes) +
		(N - zeroes) * zeroes
	for i = 1..N
		if X[i] mod P == 0
			continue
		for j=1..N
			if X[j] mod P == 0
				continue
			da = f(X[i])	// G^da = X[i], 0<=da<(P-1)
			db = f(X[j])	// G^db = X[j], 0<=db<(P-1)
			dab = da + db	// X[i]*X[j] = G^(da+db)
			Z[(G^dab)%P] += 1

How to calculate function f in O(1). Let’s make some precomputations, ff[G^i mod P] = i, for i in the range between 0 and P-2. f(i)=ff[i]. Let’s rewrite this code a bit, namely precalculate C_{i} - frequency of the number i in the array X.

	for i=1..N
		if X[i]%P==0
			C[X[i] % P] += 1
		else
			C[f(X[i])] += 1
	Z[0] = C[0] * C[0] + 2 * C[0] * (N - C[0])
	for i = 0..P-2
		for j = 0..P-2
			Z[(G^(i+j))%P] += C[i] * C[j]

What we see, last code is very similar to multiplication of two polynomials

	for i = 1..N
		for j = 1..M
			C[i + j] += A[i] * B[j]	//Standard multiplication of polynomials

Try to change the code above a bit

	for i = 0..P-2
		for j = 0..P-2
			F[i + j] += C[i] * C[j]
	for i = 0..2*(P-2)
		Z[(G^i)%P] += F[i]

How to multiplicate two polynomials in way faster than O(N^2), there are many algorithms, let’s choose Fast Fourier Transform. Below pseudocode of multiplication using Fast Fourier Transform.

	D = FFT(C)
	for i = 0..2 * (P-2)
		E[i] = D[i] * D[i]
	F = inverseFFT(E)
	for i = 0..2*(P-2)
		Z[(G^i)%P] += F[i]

But we have two primes 599 and 601, what to do in this case? This solution can be modified for the case with two primes, find primitive roots for 599 and 601, let’s call them G and H.

Now every number between 0 and 359998 can be in one from the four types:

(0, 0), only 0 satisfied this condition

($G^i mod 599$, 0), only 598 numbers satisfied this condition

(0, $H^i mod 601$), only 600 numbers satisfied this condition

($G^i mod 599$, $H^j mod 601$)

How to multiply two numbers from these types, we can multiply first and second components independently from each other. In other words: (A,B)*(C,D)=(A*C mod 599, B*D mod 601) from Chinese Remainder Theorem Consider all possible multiplications from types:

$(0, 0) = (0, 0) * (A, B)$ OR $(A, B) * (0, 0)$, $0 <= A <= 598, 0 <= B <= 600$

$(0, 0) = (0, A) * (B, 0)$ OR $(C, 0) * (0, D)$, $0 < A, D <= 598, 0 < B, C <= 600$

$(0, E) = (0, A) * (0, B)$ OR $(0, A) * (C, B)$ OR $(C, A) * (0, B)$, $0 < C <= 598, 0 < A, B, E <= 600$

$(E, 0) = (A, 0) * (B, 0)$ OR $(A, 0) * (B, C)$ OR $(A, C) * (B, 0)$, $0 < A, B, E <= 598$, $0 < C <= 600$

$(X, Y) = (A, B) * (C, D)$, $0 < A, C, X <= 598, 0 < B, D, Y <= 600$

First four convolutions can be processed in time O(599 * 601) = O(359999), consider last convolution(it will be between two multipliers of fourth type):

(X, Y) = (A, B) * (C, D) = (G^{i1} mod 599, H^{j1} mod 601) * (G^{i2} mod 599, H^{j2} mod 601) = (G^{i1+i2} mod 599, H^{j1+j2} mod 601)

Let’s create polynomial P, where P[i * 2 * 601 + j] is the number of numbers G^{i}*H^{j} mod 359999, if we’ll multiply it by itself with FFT, we can take data of (X, Y)

(i1 * 2 * 601 + j1) * (i2 * 2 * 601 + j2) = (i1 + i2) * 2 * 601 + (j1 + j2), 0 <= (j1 + j2) < 2 * 601

for i = 0..598
	for j = 0..600
		P[i * 2 * 601 + j] += A[i][j]	// number of numbers which represented in type G^i * H^j % 359999
fft(P)
for i = 0..(1<<21)-1
	D[i] = P[i] * P[i]	//Multiplying with FFT
Q[i] = inverseFFT(D)
for i = 0..4*359999
	i12 = i / (2 * 601)
	j12 = i % (2 * 601)
	Z[G^(i12)%599 * H^(j12)%601] += Q[i]

After we’ll find array Z in fast way. How to speed up the following code:

	ans=0
	for i = 0..359999-1
		for j = 0..359999-1
			for k = 0..359999-1
				if gcd(i, j, k) = 1
					ans += Z[i] * Z[j] * Z[k] 
					ans %= 1000 * 1000 * 1000 + 7

Let’s use some inclusion-exclusion, more precisely mobius-function, you can see this tutorial for similar problem:
But in this problem, we don’t need to order i, j, k, therefore

for i = 1..359999
	for j = i;j <= 359999; j+=i
		D[i] += Z[j];
	ans += mu[i] * D[i] * D[i] * D[i]	//number of positive integer triplets, where each number is divisible by i multiplied by mobius-function from i

We forget about zeroes. How to count it:

$i = j = k = 0$, $gcd(i, j, k) = 0$, skip it

$i * j = 0$ OR $j * k = 0$ OR $i * k = 0$, we need to add $Z_{1} * Z_{0} * Z_{0}$

$i = 0$ OR $j = 0$ OR $k = 0$, we need to find the number of pairs $(a, b)$, where $gcd(Z_{a}, Z_{b}) = 1$ and multiply it by value of $Z_{0}$

ans = 3 * Z[1] * Z[0] * Z[0]
for i=1..359999
	for j=i;j<=359999; j+=i
		D[i]+=Z[j];
	ans += mu[i] * D[i] * D[i] * D[i]	//number of positive integer triplets, where each number is divisible by i multiplied by mobius-function from i
	ans += 3 * mu[i] * D[i] * D[i]		//Don't forget about multipling by 3

Overall time complexity of this approach is O(2^{M} * M + X log X), where M = 21, X = 359999.

Solution:

Setter’s solution can be found here

Tester’s solution can be found here

Please feel free to post comments if anything is not clear to you.

[likecs]

Please upload the full editorial.

[likecs]

Very nice Editorial.

[swetankmodi]

A very good problem and editorial! Though i could not solve it during the contest, I can definitely solve it now on my own

[kevinsogo]

There’s a slight difference in my solution from the one described in this editorial, particularly in the part where you have to count the occurrences of Y_{i,j} coprime to m = 359999 = 599\cdot 601.

I’ll introduce the notation I’ll use. Let x_v be the number of indices i \in [1, N] such that X_i = v. Let y_v be the number of index pairs (i,j) \in [1, N]^2 such that Y_{i,j} = v. Our goal is to compute the values y_v for all v \in \left(\mathbb{Z}/m\mathbb{Z}\right)^\times. Here, \left(\mathbb{Z}/m\mathbb{Z}\right)^\times denotes the multiplicative group of integers mod m and consists of numbers coprime to m. This group has \phi(m) elements.

I’ll also mention that this explanation requires some standard abstract algebra facts. Nothing too deep, though.

---

If there were a number g that can generate all numbers in \left(\mathbb{Z}/m\mathbb{Z}\right)^\times, then we can use FFT to compute the y_v by rearranging the elements of \left(\mathbb{Z}/m\mathbb{Z}\right)^\times according to their logarithms base g. Specifically, every element of \left(\mathbb{Z}/m\mathbb{Z}\right)^\times is of the form g^i for i \in [0, \phi(m)), and we have the equality:

y_{g^i} = \sum_{j=0}^{\phi(m)-1} x_{g^j} \cdot x_{g^{i-j}}

for all i \in [0, \phi(m)). This is a convolution and can be computed in O(\phi(m) \log \phi(m)) time, and the details of this convolution has already been described in the editorial.

But the generator g doesn’t exist. In fact, it only exists if m is one of the forms 2, 4, p^k, or 2p^k, for k > 0 and p an odd prime.

Our case m = 359999 = 599\cdot 601 isn’t of that form. But we can modify the algorithm above to work on our m by using a near generator instead. To explain what I mean by this, note that the most that any number g can generate mod m = 359999 is \phi(m)/2, or exactly half the numbers. Let’s call such a number a near generator. An example of a near generator is g = 7. This number generates exactly half of the numbers in \left(\mathbb{Z}/m\mathbb{Z}\right)^\times. Note that g^{\phi(m)/2} \equiv 1 \pmod{m}.

An example of a number that g = 7 cannot generate is 17. Let’s define h = 17.

Let’s denote the set of numbers that g can generate as (g). Notice that:

• (g) is a subgroup of \left(\mathbb{Z}/m\mathbb{Z}\right)^\times, so if we form the cosets of (g), we get two sets: (g) itself and h\cdot (g). This implies that every element of \left(\mathbb{Z}/m\mathbb{Z}\right)^\times is exactly one of the following forms: g^i or hg^i, for i \in [0, \phi(m)/2).
• (g) is a normal subgroup of \left(\mathbb{Z}/m\mathbb{Z}\right)^\times, which means we can form the quotient \left(\mathbb{Z}/m\mathbb{Z}\right)^\times/(g). Furthermore, it must be isomorphic to \mathbb{Z}/2\mathbb{Z}, the only group of size 2 (up to isomorphism). This implies that (h\cdot (g))\cdot (h\cdot (g)) = (g). This is the same as saying h^2 \in (g). This is the same as saying h^2 = g^H for some H \in [0, \phi(m)/2). (In fact, H = 21370.)

Using these facts, we can now compute y_v for various kinds of v using the following:

\begin{aligned} y_{hg^i} &= \sum_{j=0}^{\phi(m)/2-1} x_{g^j} \cdot x_{hg^{i-j}} + \sum_{j=0}^{\phi(m)/2-1} x_{hg^j} \cdot x_{g^{i-j}} \\\ y_{g^i} &= \sum_{j=0}^{\phi(m)/2-1} x_{g^j} \cdot x_{g^{i-j}} + \sum_{j=0}^{\phi(m)/2-1} x_{hg^j} \cdot x_{hg^{i-j-H}} \end{aligned}

for all i \in [0, \phi(m)/2). Now, we get 4 convolutions which can all be computed with FFT!