OFFSET
1,1
COMMENTS
Composite numbers with a prime number of divisors.
LINKS
T. D. Noe, Table of n, a(n) for n = 1..1000
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
EXAMPLE
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.
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 *)
PROG
(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
(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
(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
CROSSREFS
KEYWORD
nonn,easy,nice
AUTHOR
STATUS
approved