# Greetings from The On-Line Encyclopedia of Integer Sequences! http://oeis.org/

Search: id:a036454
Showing 1-1 of 1

%I A036454 #51 Sep 12 2024 14:38:51
%S A036454 4,9,16,25,49,64,81,121,169,289,361,529,625,729,841,961,1024,1369,
%T A036454 1681,1849,2209,2401,2809,3481,3721,4096,4489,5041,5329,6241,6889,
%U A036454 7921,9409,10201,10609,11449,11881,12769,14641,15625,16129,17161,18769,19321
%N A036454 Prime powers with special exponents: q^(p-1) where p > 2 and q are prime numbers.
%C A036454 Composite numbers with a prime number of divisors.
%H A036454 T. D. Noe, <a href="/A036454/b036454.txt">Table of n, a(n) for n = 1..1000</a>
%F A036454 d(d(a(n))) = 2, where d(x) = tau(x) = sigma_0(x) is the number of divisors of x.
%F A036454 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
%F A036454 (1 - A010051(a(n))) * A010055(a(n)) * A010051(A100995(a(n))+1) = 1. - _Reinhard Zumkeller_, Jun 05 2013
%F A036454 A036459(a(n)) = 2. - _Ivan Neretin_, Jan 25 2016
%F A036454 a(n) = A283262(n)^2. - _Amiram Eldar_, Jul 04 2022
%F A036454 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
%e A036454 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.
%p A036454 N:= 10^5:
%p A036454 P1:= select(isprime,[2,seq(2*i+1,i=1..floor((sqrt(N)-1)/2))]):
%p A036454 P2:= select(`<=`,P1,1+ilog2(N))[2..-1]:
%p A036454 S:= {seq(seq(p^(q-1), q = select(`<=`,P2,1+floor(log[p](N)))),p=P1)}:
%p A036454 sort(convert(S,list)); # _Robert Israel_, May 18 2015
%t A036454 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 *)
%t A036454 Select[Range[20000], ! PrimeQ[#] && PrimeQ[DivisorSigma[0, #]] &] (* _Carlos Eduardo Olivieri_, May 18 2015 *)
%o A036454 (PARI) for(n=1,34000, if(isprime(n), n++,x=numdiv(n); if(isprime(x),print(n))))
%o A036454 (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))
%o A036454 \\ _Charles R Greathouse IV_, Apr 26 2012
%o A036454 (Haskell)
%o A036454 a009087 n = a009087_list !! (n-1)
%o A036454 a009087_list = filter ((== 1) . a010051 . (+ 1) . a100995) a000961_list
%o A036454 -- _Reinhard Zumkeller_, Jun 05 2013
%o A036454 (Magma) [n: n in [1..20000] | not IsPrime(n) and IsPrime(DivisorSigma(0, n))]; // _Vincenzo Librandi_, May 19 2015
%o A036454 (Python)
%o A036454 from sympy import primepi, integer_nthroot, primerange
%o A036454 def A036454(n):
%o A036454     def f(x): return int(n+x-sum(primepi(integer_nthroot(x, p-1)[0]) for p in primerange(3,x.bit_length()+1)))
%o A036454     def bisection(f,kmin=0,kmax=1):
%o A036454         while f(kmax) > kmax: kmax <<= 1
%o A036454         while kmax-kmin > 1:
%o A036454             kmid = kmax+kmin>>1
%o A036454             if f(kmid) <= kmid:
%o A036454                 kmax = kmid
%o A036454             else:
%o A036454                 kmin = kmid
%o A036454         return kmax
%o A036454     return bisection(f,n,n) # _Chai Wah Wu_, Sep 12 2024
%Y A036454 Intersection of A002808 and A009087.
%Y A036454 Cf. A000005, A010051, A010553, A010055, A100995, A036450, A036452, A036459, A283262.
%K A036454 nonn,easy,nice
%O A036454 1,1
%A A036454 _Labos Elemer_

# Content is available under The OEIS End-User License Agreement: http://oeis.org/LICENSE