%I #19 Aug 14 2024 01:51:21

%S 2,2,3,5,3,2,7,11,5,2,13,3,17,7,19,23,29,31,11,37,41,43,2,3,13,47,53,

%T 5,59,61,67,17,71,73,79,19,83,89,2,97,101,103,107,109,23,113,127,7,

%U 131,137,139,149,151,29,157,163,167,31,173,179,181,191,193,197,199,211,223

%N Primes (in order) occurring in A053810.

%H Amiram Eldar, <a href="/A053811/b053811.txt">Table of n, a(n) for n = 1..10000</a>

%F a(n) = A006530(A053810(n)) = A020639(A053810(n)). - _David Wasserman_, Feb 17 2006

%F a(n) = A053810(n)^(1/A053812(n)). - _Amiram Eldar_, Nov 21 2020

%o (PARI) LIM = prime(80)^2; v = vector(400); count = 0; forprime (p = 2, prime(80), x = 2; while (p^x <= LIM, count++; v[count] = p^x; x = nextprime(x + 1))); v = vecsort(vector(count, i, v[i])); A = vector(count); for (i = 1, count, f = factor(v[i]); A[i] = f[1, 1]); A \\ _David Wasserman_, Feb 17 2006

%o (Python)

%o from sympy import primepi, integer_nthroot, primerange, primefactors

%o def A053811(n):

%o def f(x): return int(n-1+x-sum(primepi(integer_nthroot(x, p)[0]) for p in primerange(x.bit_length())))

%o kmin, kmax = 1,2

%o while f(kmax) >= kmax:

%o kmax <<= 1

%o while True:

%o kmid = kmax+kmin>>1

%o if f(kmid) < kmid:

%o kmax = kmid

%o else:

%o kmin = kmid

%o if kmax-kmin <= 1:

%o break

%o return primefactors(kmax)[0] # _Chai Wah Wu_, Aug 13 2024

%Y Cf. A000040, A000961, A006530, A020639, A025473, A053810, A053812.

%K easy,nonn

%O 1,1

%A _Henry Bottomley_, Mar 28 2000

%E More terms from _David Wasserman_, Feb 17 2006

%E Offset corrected by _Amiram Eldar_, Nov 21 2020