cp's OEIS Frontend

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

%I A383211 #24 May 28 2025 14:19:37
%S A383211 4,8,9,25,27,32,49,64,121,125,128,169,243,289,343,361,529,729,841,961,
%T A383211 1024,1331,1369,1681,1849,2048,2187,2197,2209,2809,3125,3481,3721,
%U A383211 4489,4913,5041,5329,6241,6859,6889,7921,8192,9409,10201,10609,11449,11881,12167
%N A383211 Numbers of the form p^e where p is prime and e > 1 is squarefree.
%H A383211 Michael De Vlieger, <a href="/A383211/b383211.txt">Table of n, a(n) for n = 1..10000</a>
%F A383211 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
%t A383211 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 *)
%o A383211 (SageMath)
%o A383211 def A383211List(upto: int) -> list[int]:
%o A383211     L = []
%o A383211     for p in prime_range(2, upto + 1):
%o A383211         E = A383266(upto, p)
%o A383211         for e in range(2, E+1):
%o A383211             if is_squarefree(e):
%o A383211                 n = p^e
%o A383211                 if n <= upto:
%o A383211                     L.append(n)
%o A383211     return sorted(L)
%o A383211 print(A383211List(12222))
%o A383211 (PARI) isok(k) = {my(e = isprimepower(k)); e > 1 && issquarefree(e);} \\ _Amiram Eldar_, May 28 2025
%o A383211 (Python)
%o A383211 from math import isqrt
%o A383211 from sympy import mobius, integer_log, primerange
%o A383211 def A383211(n):
%o A383211     def bisection(f,kmin=0,kmax=1):
%o A383211         while f(kmax) > kmax: kmax <<= 1
%o A383211         kmin = kmax >> 1		
%o A383211         while kmax-kmin > 1:
%o A383211             kmid = kmax+kmin>>1
%o A383211             if f(kmid) <= kmid:
%o A383211                 kmax = kmid
%o A383211             else:
%o A383211                 kmin = kmid
%o A383211         return kmax
%o A383211     def g(x): return sum(mobius(k)*(x//k**2) for k in range(1, isqrt(x)+1))
%o A383211     def f(x): return n+x-sum(g(integer_log(x,p)[0])-1 for p in primerange(isqrt(x)+1))
%o A383211     return bisection(f,n,n) # _Chai Wah Wu_, May 28 2025
%Y A383211 Cf. A005117, A383266, A144338, A053810.
%K A383211 nonn
%O A383211 1,1
%A A383211 _Peter Luschny_, Apr 21 2025