login
Prime powers p^e where p is a prime and e >= 2 (prime powers without the primes or 1).
107

%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