OFFSET
1,2
LINKS
Daniel Forgues, Table of n, a(n) for n=1..10000
FORMULA
(i) a(n) < n for n > 2. (ii) a(n)/n is bounded and lim sup a(n)/n must be around 0.7. (iii) Sum_{k=1..n} a(k) seems to be asymptotic to c*n^2 with c around 0.29. (iv) a(n) = 2 if n is in A070228 (proof seems self-evident), hence there is no asymptotic expression for a(n) (just the average in (iii)). - Benoit Cloitre, Oct 14 2002
EXAMPLE
a(5) = 2 because pp(5) = 16 = 2^4 (not 4^2 as we take the smallest base).
MATHEMATICA
pp = Select[ Range[5000], Apply[GCD, Last[ Transpose[ FactorInteger[ # ]]]] > 1 &]; f[n_] := Block[{b = 2}, While[ !IntegerQ[ Log[b, pp[[n]]]], b++ ]; b]; Join[{1}, Table[ f[n], {n, 2, 80}]]
(* Second program: *)
Prepend[DeleteCases[#, 0], 1] &@ Table[If[Set[e, GCD @@ #[[All, -1]]] > 1, Power[n, 1/e], 0] &@ FactorInteger@ n, {n, 4000}] (* Michael De Vlieger, Apr 25 2017 *)
PROG
(Haskell)
a025478 n = a025478_list !! (n-1) -- a025478_list defined in A001597.
-- Reinhard Zumkeller, Mar 11 2014
(Python)
from math import gcd
from sympy import mobius, integer_nthroot, factorint
def A025478(n):
if n == 1: return 1
def f(x): return int(n-2+x+sum(mobius(k)*(integer_nthroot(x, k)[0]-1) for k in range(2, 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 integer_nthroot(kmax, gcd(*factorint(kmax).values()))[0] # Chai Wah Wu, Aug 13 2024
(PARI) lista(kmax) = {my(r, e); print1(1, ", "); for(k = 1, kmax, e = ispower(k, , &r); if(e > 0, print1(r, ", "))); } \\ Amiram Eldar, Sep 07 2024
CROSSREFS
KEYWORD
easy,nonn
AUTHOR
EXTENSIONS
Definition edited and cross-reference added by Daniel Forgues, Mar 10 2009
STATUS
approved