(PARI) primeslessxn2(n, p) = /*primes < p^n if prime */ { sr=0; for(x=0, n, y=primepi(p^x); if(y<>0, sr+=1.0/y); if(isprime(y), print1(y", ")); ); print(); print(sr) }
(Python)
from itertools import islice
from sympy import isprime, primepi
def
A087865_gen(): # generator of terms
k = 1
while True:
k *= 3
if isprime(p:=primepi(k)):
yield p
