A330669 - OEIS

%I #26 Aug 20 2024 14:58:40

%S 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,

%T 20,21,22,2,23,24,25,26,27,28,29,30,5,3,31,1,32,33,34,35,36,37,38,39,

%U 6,40,41,42,43,44,45,46,47,48,49

%N The prime indices of the prime powers (A000961).

%H Michael De Vlieger, <a href="/A330669/b330669.txt">Table of n, a(n) for n = 1..10000</a>

%F a(n) = A000720(A025473(n)). - _Michel Marcus_, Dec 24 2019

%F A000040(a(n))^A025474(n) = A000961(n) for n > 0. - _Alois P. Heinz_, Feb 20 2020

%e 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.

%p b:= proc(n) option remember; local k; for k from

%p       1+b(n-1) while nops(ifactors(k)[2])>1 do od; k

%p     end: b(1):=1:

%p a:= n-> `if`(n=1, 0, numtheory[pi](ifactors(b(n))[2, 1$2])):

%p seq(a(n), n=1..100);  # _Alois P. Heinz_, Feb 20 2020

%t mxn = 500; Join[{0}, Transpose[ Sort@ Flatten[ Table[ {Prime@n^ex, n}, {n, PrimePi@ mxn}, {ex, Log[Prime@n, mxn]}], 1]][[2]]]

%o (PARI) lista(nn) = {print1(0); for(n=2, nn, if(isprimepower(n, &p), print1(", ", primepi(p)))); } \\ _Jinyuan Wang_, Feb 19 2020

%o (Python)

%o from sympy import primepi, integer_nthroot, primefactors

%o def A330669(n):

%o     if n == 1: return 0

%o     def f(x): return int(n-2+x-sum(primepi(integer_nthroot(x,k)[0]) for k in range(1,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 int(primepi(primefactors(kmax)[0])) # _Chai Wah Wu_, Aug 20 2024

%Y Cf. A000040, A000961, A000720, A025473, A025474, A065515.

%K easy,nonn

%O 1,3

%A Grant E. Martin and _Robert G. Wilson v_, Dec 23 2019