login
A053810
Prime powers of prime numbers.
47
4, 8, 9, 25, 27, 32, 49, 121, 125, 128, 169, 243, 289, 343, 361, 529, 841, 961, 1331, 1369, 1681, 1849, 2048, 2187, 2197, 2209, 2809, 3125, 3481, 3721, 4489, 4913, 5041, 5329, 6241, 6859, 6889, 7921, 8192, 9409, 10201, 10609, 11449, 11881, 12167
OFFSET
1,1
FORMULA
a(n) = A053811(n)^A053812(n). - David Wasserman, Feb 17 2006
A010055(a(n)) * A010051(A100995(a(n))) = 1. - Reinhard Zumkeller, Jun 05 2013
Sum_{n>=1} 1/a(n) = Sum_{p prime} P(p) = 0.6716752222..., where P is the prime zeta function. - Amiram Eldar, Nov 21 2020
MATHEMATICA
pp={}; Do[if=FactorInteger[n]; If[Length[if]==1&&PrimeQ[if[[1, 1]]]&&PrimeQ[if[[1, 2]]], pp=Append[pp, n]], {n, 13000}]; pp
Sort[ Flatten[ Table[ Prime[n]^Prime[i], {n, 1, PrimePi[ Sqrt[12800]]}, {i, 1, PrimePi[ Log[ Prime[n], 12800]]}]]]
PROG
(PARI) is(n)=isprime(isprimepower(n)) \\ Charles R Greathouse IV, Mar 19 2013
(Haskell)
a053810 n = a053810_list !! (n-1)
a053810_list = filter ((== 1) . a010051 . a100995) $ tail a000961_list
-- Reinhard Zumkeller, Jun 05 2013
(Python)
from sympy import primepi, integer_nthroot, primerange
def A053810(n):
def f(x): return int(n-1+x-sum(primepi(integer_nthroot(x, p)[0]) for p in primerange(x.bit_length())))
kmin, kmax = 1, 2
while f(kmax) >= kmax:
kmax <<= 1
while True:
kmid = kmax+kmin>>1
if f(kmid) < kmid:
kmax = kmid
else:
kmin = kmid
if kmax-kmin <= 1:
break
return kmax # Chai Wah Wu, Aug 13 2024
CROSSREFS
Cf. A203967; subsequence of A000961.
Cf. A113877 (similar for semiprimes).
Sequence in context: A371223 A051676 A114129 * A076702 A051761 A153326
KEYWORD
easy,nonn
AUTHOR
Henry Bottomley, Mar 28 2000
EXTENSIONS
More terms from David Wasserman, Feb 17 2006
STATUS
approved