%I #85 Sep 17 2024 10:36:47
%S 4,8,9,16,25,27,32,49,64,81,121,125,128,169,243,256,289,343,361,512,
%T 529,625,729,841,961,1024,1331,1369,1681,1849,2048,2187,2197,2209,
%U 2401,2809,3125,3481,3721,4096,4489,4913,5041,5329,6241,6561,6859,6889,7921,8192,9409,10201,10609,11449,11881,12167,12769,14641
%N Prime powers p^e where p is a prime and e >= 2 (prime powers without the primes or 1).
%C These are sometimes called the proper prime powers.
%C A proper subset of A001597.
%C Equals A000961 \ A008578 = { x in A001597 | A001221(x)=1 }. - _M. F. Hasler_, Aug 29 2014
%H Jens Kruse Andersen, <a href="/A246547/b246547.txt">Table of n, a(n) for n = 1..10000</a>
%H Chai Wah Wu, <a href="https://arxiv.org/abs/2409.05844">Algorithms for complementary sequences</a>, arXiv:2409.05844 [math.NT], 2024.
%F a(n) = A025475(n+1). - _M. F. Hasler_, Aug 29 2014
%F Sum_{n>=1} 1/a(n) = Sum_{p prime} 1/(p*(p-1)) = A136141. - _Amiram Eldar_, Dec 21 2020
%p isA246547 := proc(n)
%p local ifs;
%p ifs := ifactors(n)[2] ;
%p if nops(ifs) <> 1 then
%p false;
%p else
%p is(op(2, op(1, ifs)) > 1);
%p end if;
%p end proc:
%p for n from 2 do
%p if isA246547(n) then
%p print(n) ;
%p end if;
%p end do: # _R. J. Mathar_, Feb 01 2016 # Or:
%p isA246547 := n -> not isprime(n) and nops(numtheory:-factorset(n)) = 1:
%p select(isA246547, [$1..10000]); # _Peter Luschny_, May 01 2018
%t With[{upto=15000},Complement[Select[Range[upto],PrimePowerQ],Prime[ Range[ PrimePi[ upto]]]]] (* _Harvey P. Dale_, Nov 28 2014 *)
%t Select[ Range@ 15000, PrimePowerQ@# && !SquareFreeQ@# &] (* _Robert G. Wilson v_, Dec 01 2014 *)
%t With[{n = 15000}, Union@ Flatten@ Table[Array[p^# &, Floor@ Log[p, n] - 1, 2], {p, Prime@ Range@ PrimePi@ Sqrt@ n}] ] (* _Michael De Vlieger_, Jul 06 2018, faster program *)
%o (PARI) for(n=1,10^5,if(isprimepower(n)>=2,print1(n,", ")));
%o (PARI) m=10^5; v=[]; forprime(p=2, sqrtint(m), e=2; while(p^e<=m, v=concat(v, p^e); e++)); v=vecsort(v) \\ Faster program. _Jens Kruse Andersen_, Aug 29 2014
%o (SageMath)
%o def A246547List(n):
%o return [p for p in srange(2, n) if p.is_prime_power() and not p.is_prime()]
%o print(A246547List(14642)) # _Peter Luschny_, Sep 16 2023
%o (Python)
%o from sympy import primepi, integer_nthroot
%o def A246547(n):
%o def f(x): return int(n-1+x-sum(primepi(integer_nthroot(x,k)[0]) for k in range(2,x.bit_length())))
%o kmin, kmax = 1,2
%o while f(kmax) >= kmax:
%o kmax <<= 1
%o while True:
%o kmid = kmax+kmin>>1
%o if f(kmid) < kmid:
%o kmax = kmid
%o else:
%o kmin = kmid
%o if kmax-kmin <= 1:
%o break
%o return kmax # _Chai Wah Wu_, Aug 14 2024
%Y Essentially the same as A025475.
%Y Cf. A000961, A001597, A025528, A051953, A136141, A246655.
%Y There are four different sequences which may legitimately be called "prime powers": A000961 (p^k, k >= 0), A246655 (p^k, k >= 1), A246547 (p^k, k >= 2), A025475 (p^k, k=0 and k >= 2). When you refer to "prime powers", be sure to specify which of these you mean. Also A001597 is the sequence of nontrivial powers n^k, n >= 1, k >= 2. - _N. J. A. Sloane_, Mar 24 2018
%K nonn,easy
%O 1,1
%A _Joerg Arndt_, Aug 29 2014