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