A007947(n) = factorback(factorint(n)[, 1])

present = 0	\\ mask of present values
w = -1	\\ previous term
p = 1	\\ radical of previous term
pp = 1	\\ radical of previous previous term

{
	for (n=1, 10 001,
		forbidden = gcd(p, pp);
		mandatory = p / gcd(p, pp);

		a = mandatory;
		while (bittest(present, a)>0 || gcd(forbidden, a)>1,
			a += mandatory;
		);

		if (n > 1,
			print (n-1 " " gcd(w, a));
		);

		w = a;

		present += 2^a;
		pp = p;
		p = A007947(a);
	)
}

quit