A007947(n) = factorback(factorint(n)[, 1]) present = 0 p = 1 \\ previous term pp = 1 \\ previous previous term { for (n=1, 10 000, forbidden = gcd(p, pp); mandatory = p / gcd(p, pp); a = mandatory; while (bittest(present, a)>0 || gcd(forbidden, a)>1, a += mandatory; ); print (n " " a); present += 2^a; pp = p; p = A007947(a); ) } quit