login
A330669
The prime indices of the prime powers (A000961).
1
0, 1, 2, 1, 3, 4, 1, 2, 5, 6, 1, 7, 8, 9, 3, 2, 10, 11, 1, 12, 13, 14, 15, 4, 16, 17, 18, 1, 19, 20, 21, 22, 2, 23, 24, 25, 26, 27, 28, 29, 30, 5, 3, 31, 1, 32, 33, 34, 35, 36, 37, 38, 39, 6, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49
OFFSET
1,3
LINKS
FORMULA
a(n) = A000720(A025473(n)). - Michel Marcus, Dec 24 2019
A000040(a(n))^A025474(n) = A000961(n) for n > 0. - Alois P. Heinz, Feb 20 2020
EXAMPLE
a(16) is 2 since A000961(16) is 27 which is 3^3 = (p_2)^3, i.e., the prime index of 3 is 2.
MAPLE
b:= proc(n) option remember; local k; for k from
1+b(n-1) while nops(ifactors(k)[2])>1 do od; k
end: b(1):=1:
a:= n-> `if`(n=1, 0, numtheory[pi](ifactors(b(n))[2, 1$2])):
seq(a(n), n=1..100); # Alois P. Heinz, Feb 20 2020
MATHEMATICA
mxn = 500; Join[{0}, Transpose[ Sort@ Flatten[ Table[ {Prime@n^ex, n}, {n, PrimePi@ mxn}, {ex, Log[Prime@n, mxn]}], 1]][[2]]]
PROG
(PARI) lista(nn) = {print1(0); for(n=2, nn, if(isprimepower(n, &p), print1(", ", primepi(p)))); } \\ Jinyuan Wang, Feb 19 2020
(Python)
from sympy import primepi, integer_nthroot, primefactors
def A330669(n):
if n == 1: return 0
def f(x): return int(n-2+x-sum(primepi(integer_nthroot(x, k)[0]) for k in range(1, x.bit_length())))
kmin, kmax = 1, 2
while f(kmax) >= kmax:
kmax <<= 1
while True:
kmid = kmax+kmin>>1
if f(kmid) < kmid:
kmax = kmid
else:
kmin = kmid
if kmax-kmin <= 1:
break
return int(primepi(primefactors(kmax)[0])) # Chai Wah Wu, Aug 20 2024
CROSSREFS
KEYWORD
easy,nonn
AUTHOR
Grant E. Martin and Robert G. Wilson v, Dec 23 2019
STATUS
approved