A383211 Numbers of the form p^e where p is prime and e > 1 is squarefree.

4, 8, 9, 25, 27, 32, 49, 64, 121, 125, 128, 169, 243, 289, 343, 361, 529, 729, 841, 961, 1024, 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

Author

Peter Luschny, Apr 21 2025

Keywords

• nonn

Links

• Michael De Vlieger, Table of n, a(n) for n = 1..10000

Crossrefs

Cf. A005117, A383266, A144338, A053810.

Programs

Mathematica

lmt = 12500; Sort[ Select[ Flatten[ Table[ Prime[p]^If[ SquareFreeQ@ exp, exp, 0], {p, PrimePi@ Sqrt@ lmt}, {exp, 2, Log[Prime@ p, lmt]} ]], # != 1 &]] (* Robert G. Wilson v, May 05 2025 *)

PARI

isok(k) = {my(e = isprimepower(k)); e > 1 && issquarefree(e);} \\ Amiram Eldar, May 28 2025

Python

from math import isqrt
from sympy import mobius, integer_log, primerange
def A383211(n):
    def bisection(f,kmin=0,kmax=1):
        while f(kmax) > kmax: kmax <<= 1
        kmin = kmax >> 1		
        while kmax-kmin > 1:
            kmid = kmax+kmin>>1
            if f(kmid) <= kmid:
                kmax = kmid
            else:
                kmin = kmid
        return kmax
    def g(x): return sum(mobius(k)*(x//k**2) for k in range(1, isqrt(x)+1))
    def f(x): return n+x-sum(g(integer_log(x,p)[0])-1 for p in primerange(isqrt(x)+1))
    return bisection(f,n,n) # Chai Wah Wu, May 28 2025

SageMath

def A383211List(upto: int) -> list[int]:
    L = []
    for p in prime_range(2, upto + 1):
        E = A383266(upto, p)
        for e in range(2, E+1):
            if is_squarefree(e):
                n = p^e
                if n <= upto:
                    L.append(n)
    return sorted(L)
print(A383211List(12222))

Formula

Sum_{n>=1} 1/a(n) = Sum_{n>=2} P(A005117(n)) = 0.68983147577186859321..., where P(s) is the prime zeta function. - Amiram Eldar, May 28 2025