cp's OEIS Frontend

This is a front-end for the Online Encyclopedia of Integer Sequences, made by Christian Perfect. The idea is to provide OEIS entries in non-ancient HTML, and then to think about how they're presented visually. The source code is on GitHub.

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

This page as a plain text file.
%I A330669 #26 Aug 20 2024 14:58:40
%S A330669 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 A330669 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 A330669 6,40,41,42,43,44,45,46,47,48,49
%N A330669 The prime indices of the prime powers (A000961).
%H A330669 Michael De Vlieger, <a href="/A330669/b330669.txt">Table of n, a(n) for n = 1..10000</a>
%F A330669 a(n) = A000720(A025473(n)). - _Michel Marcus_, Dec 24 2019
%F A330669 A000040(a(n))^A025474(n) = A000961(n) for n > 0. - _Alois P. Heinz_, Feb 20 2020
%e A330669 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 A330669 b:= proc(n) option remember; local k; for k from
%p A330669       1+b(n-1) while nops(ifactors(k)[2])>1 do od; k
%p A330669     end: b(1):=1:
%p A330669 a:= n-> `if`(n=1, 0, numtheory[pi](ifactors(b(n))[2, 1$2])):
%p A330669 seq(a(n), n=1..100);  # _Alois P. Heinz_, Feb 20 2020
%t A330669 mxn = 500; Join[{0}, Transpose[ Sort@ Flatten[ Table[ {Prime@n^ex, n}, {n, PrimePi@ mxn}, {ex, Log[Prime@n, mxn]}], 1]][[2]]]
%o A330669 (PARI) lista(nn) = {print1(0); for(n=2, nn, if(isprimepower(n, &p), print1(", ", primepi(p)))); } \\ _Jinyuan Wang_, Feb 19 2020
%o A330669 (Python)
%o A330669 from sympy import primepi, integer_nthroot, primefactors
%o A330669 def A330669(n):
%o A330669     if n == 1: return 0
%o A330669     def f(x): return int(n-2+x-sum(primepi(integer_nthroot(x,k)[0]) for k in range(1,x.bit_length())))
%o A330669     kmin, kmax = 1,2
%o A330669     while f(kmax) >= kmax:
%o A330669         kmax <<= 1
%o A330669     while True:
%o A330669         kmid = kmax+kmin>>1
%o A330669         if f(kmid) < kmid:
%o A330669             kmax = kmid
%o A330669         else:
%o A330669             kmin = kmid
%o A330669         if kmax-kmin <= 1:
%o A330669             break
%o A330669     return int(primepi(primefactors(kmax)[0])) # _Chai Wah Wu_, Aug 20 2024
%Y A330669 Cf. A000040, A000961, A000720, A025473, A025474, A065515.
%K A330669 easy,nonn
%O A330669 1,3
%A A330669 Grant E. Martin and _Robert G. Wilson v_, Dec 23 2019