{ mx = 10 000; sf = select(issquarefree, [1..mx]); for (n = 1, mx, f = factor(n); i = prod(k = 1, #f~, f[k, 1]^(f[k, 2] \ 2)); j = setsearch(sf, n / i^2); v = j^2 * sf[i]; print (n " " v); ); } quit