Read problems statements in Mandarin Chinese and Vietnamese as well. 

It is well-known that ∑sqrt(a_{i}), a_{i} ∈ N is a root of some integer-coefficient polynomial. For example: sqrt(2) is a root of polynomial: x^{2} − 2. Now, your task is to find not only such polynomial but also the minimal one. When comparing two polynomials, firstly, we consider their degree and then the coefficient of the highest term, and then the second highest term and so on. (Note that we consider only polynomial with the coefficient of highest term is positive)

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

First line contains an integer T denotes the number of test cases. Each test case is described by an integer n and followed by n space-separated integers a_{i}.

------ Output: ------ 

Each test case print an integer k denotes the degree of polynomial in a single line. Then next line print k + 1  space-separated integers modulo 10^{9} + 7, coefficients from lowest to highest term.

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

1 ≤ T ≤ 5, 1 ≤ n ≤ 15
a_{i} are n distinct primes, 1 ≤ a_{i} ≤ 10^{9}

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

$Subtask #1: (10 points):  n ≤ 3$
$Subtask #2: (20 points):  n ≤ 5$
$Subtask #3: (30 points):  n ≤ 10$
$Subtask #4: (40 points):  n ≤ 15$

------ Example ------ 

Input:

2
1
2
2
2 3

Output:

2
1000000005 0 1
4
1 0 999999997 0 1

------ Explanation ------ 

The first polynomial is x^{2} − 2, and the second one is x^{4} − 10x^{2} + 1.

Difficulty: VERY_HARD
Tags: Polynomials and generating functions

Mathematics

Starter code:

Solutions:
P = 10 ** 9 + 7
from math import ceil, log

def addpolymod(p, q):
	l = max(len(p), len(q))
	r = p + [0] * (l - len(p))
	s = q + [0] * (l - len(q))
	return [(r[i] + s[i]) % P for i in range(l)]

def multpolymod(p, q):
	l = len(p) + len(q) - 1
	r = [0] * l
	if len(p) <= len(q):
		for i in range(len(p)):
			for j in range(i, i + len(q)):
				r[j] = (r[j] + p[i] * q[j - i] % P) % P
	else:
		for i in range(len(q)):
			for j in range(i, i + len(p)):
				r[j] = (r[j] + q[i] * p[j - i] % P) % P
	return r

def xgcd(b, n):
	(x0, x1, y0, y1) = (1, 0, 0, 1)
	while n != 0:
		(q, b, n) = (b // n, n, b % n)
		(x0, x1) = (x1, x0 - q * x1)
		(y0, y1) = (y1, y0 - q * y1)
	return (b, x0, y0)

def modinv(a, m):
	(g, x, y) = xgcd(a, m)
	return x % m
p1 = 380 * 2 ** 15 + 1
p2 = 404 * 2 ** 15 + 1
p3 = 756 * 2 ** 15 + 1
BIG = p1 * p2 * p3
primess = [p1, p2, p3]
unitroot = {p1: 623, p2: 498, p3: 770}
PP = [0, 0, 0]
PP[0] = modinv(p2 * p3, p1) * p2 * p3 % BIG
PP[1] = modinv(p1 * p3, p2) * p1 * p3 % BIG
PP[2] = modinv(p1 * p2, p3) * p1 * p2 % BIG
w = [[], [], []]
wi = [[], [], []]
m = []
mi = [[], [], []]
for i in range(17):
	m.append(1 << i)
	for j in range(3):
		w[j].append(pow(unitroot[primess[j]], 1 << 16 - i, primess[j]))
		wi[j].append(modinv(w[j][i], primess[j]))
		mi[j].append(modinv(m[i], primess[j]))
fact = [1]
infact = [1]
MAX = m[14] + 1
temp = 1
intemp = 1
for i in range(1, MAX):
	temp = temp * i % P
	intemp = intemp * modinv(i, P) % P
	fact += [temp]
	infact += [intemp]

def nspot(l):
	return int(ceil(log(l, 2)))

def chinese_remainder(K):
	ans = 0
	for i in range(3):
		ans = (ans + K[i] * PP[i] % BIG) % BIG
	return ans % P

def fft(t, p, i):
	if len(p) == 1:
		return [p[0]] * m[t]
	a = fft(t - 1, p[0::2], i)
	b = fft(t - 1, p[1::2], i)
	ret = [0] * m[t]
	x = 1
	for k in range(m[t - 1]):
		ret[k] = (a[k] + x * b[k]) % primess[i]
		ret[k + m[t - 1]] = (a[k] - x * b[k]) % primess[i]
		x = x * w[i][t] % primess[i]
	return ret

def ifft(t, p, i):
	if len(p) == 1:
		return [p[0]] * m[t]
	a = ifft(t - 1, p[0::2], i)
	b = ifft(t - 1, p[1::2], i)
	ret = [0] * m[t]
	x = 1
	for k in range(m[t - 1]):
		ret[k] = (a[k] + x * b[k]) % primess[i]
		ret[k + m[t - 1]] = (a[k] - x * b[k]) % primess[i]
		x = x * wi[i][t] % primess[i]
	return ret

def square(t, p1, i):
	a = fft(t, p1, i)
	p3 = [a[k] * a[k] % primess[i] for k in range(m[t])]
	final = ifft(t, p3, i)
	return [x * mi[i][t] % primess[i] for x in final]

def squarefast(p):
	l = 2 * len(p) - 1
	k = nspot(l)
	Y = [square(k, p, j) for j in range(3)]
	return [chinese_remainder([Y[0][i], Y[1][i], Y[2][i]]) for i in range(l)]

def squarepoly(p):
	l = 2 * len(p) - 1
	if len(p) < 17:
		return multpolymod(p, p)
	return squarefast(p)

def ftrans(f, t):
	lf = len(f)
	fn = [fact[lf - 1 - i] * f[lf - 1 - i] % P for i in range(lf)]
	gn = [infact[i] * pow(t, i // 2, P) % P for i in range(lf)]
	dlf = 2 * (lf - 1)
	ResH = [0] * (lf - 1) + [(fn[-1] * gn[0] % P + gn[-1] * fn[0] % P) % P]
	del fn[-1]
	del gn[-1]
	H = [0] * (lf + 1)
	if lf < 63:
		H1 = multpolymod(fn[0::2], gn[0::2]) + [0] * 2
		H2 = multpolymod(fn[0::2], gn[1::2]) + [0] * 2
		for i in range(1 + (lf - 1) // 2):
			H[2 * i] = H1[i]
			H[2 * i + 1] = H2[i]
	else:
		l = lf - 2
		k = nspot(l)
		fftfn = [fft(k, fn[0::2], j) for j in range(3)]
		fftgne = [fft(k, gn[0::2], j) for j in range(3)]
		fftgno = [fft(k, gn[1::2], j) for j in range(3)]
		pointwise = [[fftfn[j][i] * fftgne[j][i] % primess[j] for i in range(m[k])] for j in range(3)]
		pointwiso = [[fftfn[j][i] * fftgno[j][i] % primess[j] for i in range(m[k])] for j in range(3)]
		ifftpwe = [[mi[j][k] * x % primess[j] for x in ifft(k, pointwise[j], j)] for j in range(3)]
		ifftpwo = [[mi[j][k] * x % primess[j] for x in ifft(k, pointwiso[j], j)] for j in range(3)]
		for i in range(1 + (lf - 1) // 2):
			H[2 * i] = chinese_remainder([ifftpwe[0][i], ifftpwe[1][i], ifftpwe[2][i]])
			H[2 * i + 1] = chinese_remainder([ifftpwo[0][i], ifftpwo[1][i], ifftpwo[2][i]])
	hn = addpolymod(ResH, H)
	s = [infact[i] * hn[lf - 1 - i] % P for i in range(lf)]
	g = s[0::2]
	h = s[1::2]
	ResA = addpolymod([0] * (lf - 1) + [g[-1] * g[-1] % P], [0] * ((lf - 1) // 2) + [2 * g[i] * g[-1] % P for i in range((lf - 1) // 2)])
	del g[-1]
	a = addpolymod(ResA, squarepoly(g))
	b = [0] + squarepoly(h)
	c = [0] * (2 * len(a) - 1)
	b += [0] * (len(a) - len(b))
	for i in range(len(a)):
		c[2 * i] = (a[i] - t * b[i] % P) % P
	return c

def minpolymod(n):
	f = [-n[0] % P, 0, 1]
	for i in range(1, len(n)):
		f = ftrans(f, n[i])
	return f
T = int(input())
for i in range(T):
	n = int(input())
	a = list(map(int, input().split()))
	print(m[n])
	print(' '.join((str(x) for x in minpolymod(a))))