A330669 The prime indices of the prime powers (A000961).

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

Author

Grant E. Martin and Robert G. Wilson v, Dec 23 2019

Keywords

• easy
• nonn

Examples

			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.
		

Links

• Michael De Vlieger, Table of n, a(n) for n = 1..10000

Crossrefs

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

Programs

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]]]

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

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