A025476 - OEIS

%I #35 Aug 16 2024 08:37:11

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

%T 43,2,3,13,47,7,53,5,59,61,2,67,17,71,73,79,3,19,83,89,2,97,101,103,

%U 107,109,23,113,11,5,127,2,7,131,137,139,3,149,151,29,157,163,167,13,31,173,179

%N Prime root of n-th nontrivial prime power (A025475, A246547).

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

%p cvm := proc(n, level) local f,opf; if n < 2 then RETURN() fi;

%p f := ifactors(n); opf := op(1,op(2,f)); if nops(op(2,f)) > 1 or

%p op(2,opf) <= level then RETURN() fi; op(1,opf) end:

%p A025476_list := n -> seq(cvm(i,1),i=1..n); # n is search limit

%p A025476_list(30000);  # _Peter Luschny_, Sep 21 2011

%p # Alternative:

%p isA246547 := n -> n > 1 and not isprime(n) and type(n, 'primepower'):

%p seq(ifactors(p)[2][1][1], p in select(isA246547, [$1..30000])); # _Peter Luschny_, Jul 15 2023

%t Transpose[ Flatten[ FactorInteger[ Select[ Range[2, 30000], !PrimeQ[ # ] && Mod[ #, # - EulerPhi[ # ]] == 0 &]], 1]][[1]] (* _Robert G. Wilson v_ *)

%o (PARI) forcomposite(n=4,10^5,if( ispower(n, , &p) && isprime(p), print1(p,", "))) \\ _Joerg Arndt_, Sep 11 2021

%o (Python)

%o from sympy import primepi, integer_nthroot, primefactors

%o def A025476(n):

%o     def f(x): return int(n-1+x-sum(primepi(integer_nthroot(x,k)[0]) for k in range(2,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 15 2024

%Y Cf. A025473, A025475, A048148, A246547.

%K easy,nonn

%O 1,1

%A _David W. Wilson_