\\ Very crude script demonstrating a factorization-free method of computing terms of the sequence. \\ Compute the requested term of A190880. a(n,lim=25000)={ local(v=List()); foralmostprime(2,lim,n,u-> if(is(u), listput(v,vecprod(u)) ) ); v=vecsort(Vec(v)); if(#v, v[1] , print("Testing up to ", lim<<=1); a(n,lim) ) }; \\ Is k divisible by the square of the sum of its prime divisors, where v is a vector of the prime power factorization of k? is(v)={ my(t); vecprod(v)%sum(i=1,#v,t=v[i];ispower(t,,&t);t)^2==0 }; \\ Loop over prime powers (not in order). forpp(a,b,ff)={ my(t); forprime(p=2,sqrt(b),t=p;while(t<a,t*=p);while(t<=b,ff(t);t*=p)); forprime(p=max(a,floor(sqrt(b))+1),b,ff(p)) }; \\ Loop over k-almost primes (not in order). See, e.g., A007774 which corresponds to k=2. foralmostprime(a,b,k,ff,v=[])={ \\print("Called with: ",[a,b,k,ff,v]); my(newK=k-1, newV=vector(#v+1)); if (#v>1, for(i=1,#v-1, if(gcd(v[i],v[#v])!=1, return()) ) ); for(i=1,#v,newV[i]=v[i]); if (k > 1, if(#v, forpp(v[#v]+1,b^(1/k),p-> newV[#newV]=p; foralmostprime(max(p+1,a\p), b\p, newK, ff, newV) ) , forpp(2,b^(1/k),p-> newV[#newV]=p; foralmostprime(max(p+1,a\p), b\p, newK, ff, newV) ) ) , my(t); forprime(p=2,sqrt(b), for(i=1,#v,if(v[i]%p==0, next(2))); t=p; while(t<a,t*=p); while(t<=b, newV[#newV]=t; ff(newV); t*=p ) ); forprime(p=max(a,floor(sqrt(b))+1),b, newV[#newV]=p; ff(newV) ) ) };