A036454 Prime powers with special exponents: q^(p-1) where p > 2 and q are prime numbers.

4, 9, 16, 25, 49, 64, 81, 121, 169, 289, 361, 529, 625, 729, 841, 961, 1024, 1369, 1681, 1849, 2209, 2401, 2809, 3481, 3721, 4096, 4489, 5041, 5329, 6241, 6889, 7921, 9409, 10201, 10609, 11449, 11881, 12769, 14641, 15625, 16129, 17161, 18769, 19321

Offset: 1

Author

Labos Elemer

Keywords

• nonn
• easy
• nice

Comments

Composite numbers with a prime number of divisors.

Examples

			From powers of 2: 4,16,64,1024,4096,65536 are in the sequence since exponent+1 is also prime. The same powers of any prime base are also included.
		

Links

• T. D. Noe, Table of n, a(n) for n = 1..1000

Crossrefs

Intersection of A002808 and A009087.

Cf. A000005, A010051, A010553, A010055, A100995, A036450, A036452, A036459, A283262.

Programs

Haskell

a009087 n = a009087_list !! (n-1)
a009087_list = filter ((== 1) . a010051 . (+ 1) . a100995) a000961_list
-- Reinhard Zumkeller, Jun 05 2013

Magma

[n: n in [1..20000] | not IsPrime(n) and IsPrime(DivisorSigma(0, n))]; // Vincenzo Librandi, May 19 2015

Maple

N:= 10^5:
P1:= select(isprime,[2,seq(2*i+1,i=1..floor((sqrt(N)-1)/2))]):
P2:= select(`<=`,P1,1+ilog2(N))[2..-1]:
S:= {seq(seq(p^(q-1), q = select(`<=`,P2,1+floor(log[p](N)))),p=P1)}:
sort(convert(S,list)); # Robert Israel, May 18 2015

Mathematica

specialPrimePowerQ[n_] := With[{f = FactorInteger[n]}, Length[f] == 1 && PrimeQ[f[[1, 1]]] && f[[1, 2]] > 1 && PrimeQ[f[[1, 2]] + 1]]; Select[Range[20000], specialPrimePowerQ]  (* Jean-François Alcover, Jul 02 2013 *)
Select[Range[20000], ! PrimeQ[#] && PrimeQ[DivisorSigma[0, #]] &] (* Carlos Eduardo Olivieri, May 18 2015 *)

PARI

for(n=1,34000, if(isprime(n), n++,x=numdiv(n); if(isprime(x),print(n))))

PARI

list(lim)=my(v=List(),t);lim=lim\1+.5;forprime(p=3,log(lim)\log(2) +1, t=p-1; forprime(q=2,lim^(1/t),listput(v,q^t))); vecsort(Vec(v))
\\ Charles R Greathouse IV, Apr 26 2012

Python

from sympy import primepi, integer_nthroot, primerange
def A036454(n):
    def f(x): return int(n+x-sum(primepi(integer_nthroot(x, p-1)[0]) for p in primerange(3,x.bit_length()+1)))
    def bisection(f,kmin=0,kmax=1):
        while f(kmax) > kmax: kmax <<= 1
        while kmax-kmin > 1:
            kmid = kmax+kmin>>1
            if f(kmid) <= kmid:
                kmax = kmid
            else:
                kmin = kmid
        return kmax
    return bisection(f,n,n) # Chai Wah Wu, Sep 12 2024

Formula

d(d(a(n))) = 2, where d(x) = tau(x) = sigma_0(x) is the number of divisors of x.

a(n) = (n log n)^2 + 2n^2 log n log log n + O(n^2 log n). - Charles R Greathouse IV, Apr 26 2012

(1 - A010051(a(n))) * A010055(a(n)) * A010051(A100995(a(n))+1) = 1. - Reinhard Zumkeller, Jun 05 2013

A036459(a(n)) = 2. - Ivan Neretin, Jan 25 2016

a(n) = A283262(n)^2. - Amiram Eldar, Jul 04 2022

Sum_{n>=1} 1/a(n) = Sum_{k>=2} P(prime(k)-1) = 0.54756961912815344341..., where P is the prime zeta function. - Amiram Eldar, Jul 10 2022