rad(n) = vecprod(factor(n)[,1]~)

explore(n, seq, m, r, w) = {
	if (w==1 && #best < #seq,
		best = seq;
	);

	for (v = 1, n,
		if (!bittest(m, v) && gcd(r, v)==w,
			my (rv = rad(v));
			explore(n, concat(seq, v), m + 2^v, rv, rv/w);
		);
	);
}

row(n) = {
	my (r);
	best = [];

	for (v = 1, n,
		explore(n, [v], 2^v, r = rad(v), r);
	);

	return (best);
}

{
	k = 0;
	for (n = 1, 20,
		if (!isprime(n),
			r = row(n);
		);
		apply (v -> print (k++ " " v), r);
	);
}

quit