Read problems statements in Mandarin Chinese, Russian and Vietnamese as well. 

There are N integers X_{1}, X_{2}, ..., X_{N}.

Let's define Y_{i,j} = X_{i} × X_{j} mod 359999.

How many integer 6-tuples (a, b, c, d, e, f) are there such that:

1 ≤ a, b, c, d, e, f ≤ N
gcd(Y_{a, b}, Y_{c, d}, Y_{e, f}) = 1

We define gcd(0, 0) = 0.

------ Input ------ 

The first line of input contains an integer T denoting the number of test cases. The description of T test cases follows.
The first line of each test case contains a single integer N.
The second line contains N integers separated by single spaces: X_{1}, X_{2}, ..., X_{N}.

------ Output ------ 

For each test case, output a single line containing the answer. Since the answer can be very large, only output it modulo 10^{9} + 7.

------ Constraints ------ 

$1 ≤ T ≤ 3$
$1 ≤ X_{i} ≤ 10^{6}$

------ Subtasks ------ 

Subtask #1 (47 points):

$1 ≤ N ≤ 10^{3}$
 
Subtask #2 (53 points):

$1 ≤ N ≤ 10^{6}$
$The sum of the Ns is ≤ 10^{6}$
 
------ Example ------ 

Input:
1
3
300 3000 30000

Output:
234

Source: codechef | Difficulty: VERY_HARD
Tags: Polynomials and generating functions

Number theory

Mathematics


Solutions:
import sys
import numpy as np
p1 = 599
p2 = 601
pp = p1 * p2
MOD = 10 ** 9 + 7

def getLogPows(p, root=7):
	v = 1
	ret = np.zeros(p, dtype=np.int64)
	pw = np.zeros(p, dtype=np.int64)
	pw[0] = 1
	for i in range(p - 1):
		ret[v] = i
		v = v * root % p
		pw[i + 1] = v
	return (ret, pw)
(LOG1, POW1) = getLogPows(p1)
(LOG2, POW2) = getLogPows(p2)
CRT = np.zeros((p1, p2), dtype=np.int64)
for i in range(pp):
	CRT[i % p1, i % p2] = i
primes = np.ones(pp, dtype=np.int64)
MU = np.ones(pp, dtype=np.int64)
MU[0] = 0
primes[:2] = 0
for i in range(2, pp):
	if primes[i] == 1:
		primes[2 * i::i] = 0
		MU[i::i] *= -1
		if i * i < len(MU):
			MU[i * i::i * i] = 0

def getLogArray(xx):
	arr = np.zeros((p1 - 1, p2 - 1))
	arr10 = np.zeros(p2 - 1)
	arr20 = np.zeros(p1 - 1)
	zeros = 0
	x1 = xx % p1
	x2 = xx % p2
	for (x1, x2) in zip(xx % p1, xx % p2):
		if x1 == 0:
			if x2 == 0:
				zeros += 1
			else:
				arr10[LOG2[x2]] += 1
		elif x2 == 0:
			arr20[LOG1[x1]] += 1
		else:
			arr[LOG1[x1], LOG2[x2]] += 1
	return (arr, arr10, arr20, zeros)

def getTotals(xx):
	ret = np.zeros(pp, dtype=np.int64)
	(arr, arr10, arr20, zeros) = getLogArray(xx)
	aff = np.fft.rfft2(arr)
	aff = aff * aff
	r0 = np.fft.irfft2(aff)
	for i in range(p1 - 1):
		for j in range(p2 - 1):
			ret[CRT[POW1[i], POW2[j]]] = int(np.round(r0[i, j]))
	arrs1 = np.sum(arr, 0)
	arrs2 = np.sum(arr, 1)
	aff = np.fft.rfft(arr10)
	aff2 = np.fft.rfft(arrs1)
	r0 = np.fft.irfft(aff * aff) + 2 * np.fft.irfft(aff * aff2)
	for j in range(p2 - 1):
		ret[CRT[0, POW2[j]]] = int(np.round(r0[j]))
	aff = np.fft.rfft(arr20)
	aff2 = np.fft.rfft(arrs2)
	r0 = np.fft.irfft(aff * aff) + 2 * np.fft.irfft(aff * aff2)
	for i in range(p1 - 1):
		ret[CRT[POW1[i], 0]] = int(np.round(r0[i]))
	rr = 0
	rrr = ret % MOD
	ret = np.zeros(pp, dtype=np.int64)
	for i in range(1, pp):
		ret[i] = np.sum(rrr[i::i])
	ret %= MOD
	ret0 = ret * ret % MOD
	ret1 = ret * ret0 % MOD
	ret1 %= MOD
	rval = np.sum(ret1 * MU) % MOD
	zpairs = (zeros * len(xx) * 2 - zeros * zeros) % MOD
	zpairs += int(np.sum(arr10)) % MOD * (int(np.sum(arr20)) % MOD) * 2
	zpairs %= MOD
	ret0 = ret * ret % MOD
	rval += np.sum(ret0 * MU) % MOD * zpairs * 3 % MOD
	rval %= MOD
	rval += rrr[1] % MOD * (zpairs * zpairs % MOD) * 3 % MOD
	return rval % MOD
f = sys.stdin
t = int(f.readline())
for i in range(t):
	f.readline()
	x = np.array(list(map(int, f.readline().split())), dtype=np.int64)
	print(getTotals(x))