A246549 Prime powers p^e where p is a prime and e >= 3 (prime powers without 1, the primes, or the squares of primes).

8, 16, 27, 32, 64, 81, 125, 128, 243, 256, 343, 512, 625, 729, 1024, 1331, 2048, 2187, 2197, 2401, 3125, 4096, 4913, 6561, 6859, 8192, 12167, 14641, 15625, 16384, 16807, 19683, 24389, 28561, 29791, 32768, 50653, 59049, 65536, 68921, 78125, 79507, 83521, 103823, 117649, 130321, 131072, 148877, 161051, 177147, 205379

Offset: 1

Author

Joerg Arndt, Aug 29 2014

Keywords

• nonn

Comments

Consists of 8 and the terms of A088247. - R. J. Mathar, Sep 01 2014

Links

• Jens Kruse Andersen, Table of n, a(n) for n = 1..10000

Crossrefs

Cf. A000961, A152441, A246547, A246550.

Programs

Mathematica

With[{nn=60},Take[Union[Flatten[Table[p^Range[3,nn/3],{p,Prime[ Range[ nn]]}]]],nn]] (* Harvey P. Dale, Dec 10 2015 *)

PARI

for(n=1, 10^6, if(isprimepower(n)>=3, print1(n, ", ")));

PARI

m=10^6; v=[]; forprime(p=2, m^(1/3), e=3; while(p^e<=m, v=concat(v, p^e); e++)); v=vecsort(v) \\ Faster program. Jens Kruse Andersen, Aug 29 2014

Python

from math import isqrt
from sympy import primerange, integer_nthroot, primepi
def A246549(n):
    def g(x,a,b,c,m): yield from (((d,) for d in enumerate(primerange(b+1,isqrt(x//c)+1),a+1)) if m==2 else (((a2,b2),)+d for a2,b2 in enumerate(primerange(b+1,integer_nthroot(x//c,m)[0]+1),a+1) for d in g(x,a2,b2,c*b2,m-1)))
    def f(x): return int(n+x-sum(primepi(integer_nthroot(x, k)[0]) for k in range(3, x.bit_length())))
    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 11 2024

Formula

Sum_{n>=1} 1/a(n) = Sum_{p prime} 1/(p^2*(p-1)) = A152441. - Amiram Eldar, Oct 24 2020