M = 10 000 p = primes([1,2*M]) odious = select(p -> hammingweight(p)%2==1, p) evil = select(p -> hammingweight(p)%2==0, p) a(n) = { my (f=factor(n), i); for (k=1, #f~, i = setsearch(odious, f[k,1]); if (i, f[k,1] = evil[i], f[k,1] = odious[setsearch(evil, f[k,1])] ); ); return (factorback(f)); } for (n=1, M, print (n " " a(n))) quit