OFFSET
1,1
LINKS
Jens Kruse Andersen, Table of n, a(n) for n = 1..10000
FORMULA
Sum_{n>=1} 1/a(n) = Sum_{p prime} 1/(p^3*(p-1)) = 0.1461466097... - Amiram Eldar, Oct 24 2020
MAPLE
N:= 10^7: # to get all terms <= N
{seq(seq(p^m, m=4..floor(log[p](N))), p = select(isprime, [2, seq(2*i+1, i=1..floor(N^(1/4)))]))}; # Robert Israel, Aug 29 2014
MATHEMATICA
With[{max = 10^6}, Sort @ Flatten @ Table[p^Range[4, Floor[Log[p, max]]], {p, Select[Range[Surd[max, 4]], PrimeQ]}]] (* Amiram Eldar, Oct 24 2020 *)
PROG
(PARI) m=10^7; v=[]; forprime(p=2, m^(1/4), e=4; while(p^e<=m, v=concat(v, p^e); e++)); v=vecsort(v) \\ Jens Kruse Andersen, Aug 29 2014
(Python)
from math import isqrt
from sympy import primerange, integer_nthroot, primepi
def A246550(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(4, 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 12 2024
CROSSREFS
KEYWORD
nonn
AUTHOR
Joerg Arndt, Aug 29 2014
STATUS
approved