**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.**