OFFSET
1,1
COMMENTS
LINKS
Michael De Vlieger, Table of n, a(n) for n = 1..3351
FORMULA
Sum_{n>=1} 1/a(n) = Sum_{n>=1} 1/A088730(n) = 0.372116188498... . - Amiram Eldar, Jan 20 2024
EXAMPLE
This sequence contains prime powers of the following form:
2^2, 2^4, i.e., 2^k such that k is even.
3^3, 3^6, 3^9, i.e., 3^k such that 3 | k.
5^5, 5^10, 5^15, i.e., 5^k such that 5 | k, etc.
MAPLE
N:= 10^13: # for terms <= N
R:= NULL:
for i from 1 do
p:= ithprime(i);
if p^p > N then break fi;
R:= R, seq(p^k, k=p..floor(log[p](N)), p);
od:
sort([R]); # Robert Israel, Jan 16 2024
MATHEMATICA
nn = 10^12; i = 1; p = 2; While[p^p <= nn, p = NextPrime[p] ];
MapIndexed[Set[S[First[#2]], #1] &, Prime@ Range@ PrimePi[p] ];
Union@ Reap[
While[j = S[i];
While[S[i]^j < nn,
Sow[S[i]^j]; j += S[i] ]; j > 2,
i++] ][[-1, 1]]
PROG
(Python)
import heapq
from itertools import islice
from sympy import nextprime
def agen(): # generator of terms
v, h, m, nextp = 4, [(4, 2)], 4, 3
while True:
v, p = heapq.heappop(h)
yield v
if v >= m:
m = nextp**nextp
heapq.heappush(h, (m, nextp))
nextp = nextprime(nextp)
heapq.heappush(h, (v*p**p, p))
print(list(islice(agen(), 31))) # Michael S. Branicky, Jan 16 2024
(Python)
from sympy import integer_nthroot, primefactors
def A368107(n):
def f(x):
c = n+x
for k in range(1, x.bit_length()):
m = integer_nthroot(x, k)[0]
c -= sum(1 for p in primefactors(k) if p<=m)
return c
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 12 2024
CROSSREFS
KEYWORD
nonn,easy
AUTHOR
Michael De Vlieger, Jan 15 2024
STATUS
approved