allocate mem(2^30) is(n) = my (f=factor(n)); vecprod(f[,1]~)==vecprod(f[,2]~) lim = 29376344554936148249542656000+1 seq = [1] s = 1 explore(v,pp,x) = { if (v*if (#pp, vecprod(pp), 1)1, if (s++>#seq, seq = concat(seq, vector(#seq)); ); seq[s] = v; ), #pp==1, explore(v*pp[1]^x), fordiv (x, d, explore(v*pp[1]^d, pp[2..#pp], x/d); ); ); ); } sopf(n) = vecsum(factor(n)[,1]~) gpf(n) = if (n==1, 1, my (f=factor(n)); f[#f~, 1]) lpf(n) = if (n==1, 1, my (f=factor(n)); f[1, 1]) \\ square free numbers with sopf < L sf=[1] { forprime (p=2, logint(lim,2), sf = concat(sf, select(s -> s*lpf(s)^(sopf(s)-1) print (n++ " " v), Set(seq[1..s])); } quit