|
%I #54 Aug 16 2024 04:49:19
%S 1,8,16,27,32,64,81,125,128,216,243,256,343,512,625,729,1000,1024,
%T 1296,1331,1728,2048,2187,2197,2401,2744,3125,3375,4096,4913,5832,
%U 6561,6859,7776,8000,8192,9261,10000,10648,12167,13824,14641,15625,16384,16807
%N Perfect powers m^k where m is a positive integer and k >= 3.
%C If p|n with p prime then p^3|n.
%H Alois P. Heinz, <a href="/A076467/b076467.txt">Table of n, a(n) for n = 1..10000</a> (terms n = 1..250 from Reinhard Zumkeller)
%F For n > 1: GCD(exponents in prime factorization of a(n)) > 2, cf. A124010. - _Reinhard Zumkeller_, Apr 13 2012
%F Sum_{n>=1} 1/a(n) = 2 - zeta(2) + Sum_{k>=2} mu(k)*(2 - zeta(k) - zeta(2*k)) = 1.3300056287... - _Amiram Eldar_, Jul 02 2022
%p N:= 10^5: # to get all terms <= N
%p S:= {1, seq(seq(m^k, m = 2 .. floor(N^(1/k))),k=3..ilog2(N))}:
%p sort(convert(S,list)); # _Robert Israel_, Sep 30 2015
%t a = {1}; Do[ If[ Apply[ GCD, Last[ Transpose[ FactorInteger[n]]]] > 2, a = Append[a, n]; Print[n]], {n, 2, 17575}]; a
%t (* Second program: *)
%t n = 10^5; Join[{1}, Table[m^k, {k, 3, Floor[Log[2, n]]}, {m, 2, Floor[n^(1/k)]}] // Flatten // Union] (* _Jean-François Alcover_, Feb 13 2018, after _Robert Israel_ *)
%o (Haskell)
%o a076467 n = a076467_list !! (n-1)
%o a076467_list = 1 : filter ((> 2) . foldl1 gcd . a124010_row) [2..]
%o -- _Reinhard Zumkeller_, Apr 13 2012
%o (Haskell)
%o import qualified Data.Set as Set (null)
%o import Data.Set (empty, insert, deleteFindMin)
%o a076467 n = a076467_list !! (n-1)
%o a076467_list = 1 : f [2..] empty where
%o f xs'@(x:xs) s | Set.null s || m > x ^ 3 = f xs $ insert (x ^ 3, x) s
%o | m == x ^ 3 = f xs s
%o | otherwise = m : f xs' (insert (m * b, b) s')
%o where ((m, b), s') = deleteFindMin s
%o -- _Reinhard Zumkeller_, Jun 18 2013
%o (PARI) is(n)=ispower(n)>2||n==1 \\ _Charles R Greathouse IV_, Sep 03 2015, edited for n=1 by _M. F. Hasler_, May 26 2018
%o (PARI) A076467(lim)={my(L=List(1),lim2=logint(lim,2),m,k);for(k=3,lim2, for(m=2,sqrtnint(lim,k),listput(L, m^k);));listsort(L,1);L}
%o b076467(lim)={my(L=A076467(lim)); for(i=1,#L,print(i ," ",L[i]));} \\ _Anatoly E. Voevudko_, Sep 29 2015, edited by _M. F. Hasler_, May 25 2018
%o (PARI) A076467_vec(LIM,S=List(1))={for(x=2,sqrtnint(LIM,3),for(k=3, logint(LIM, x), listput(S, x^k))); Set(S)} \\ _M. F. Hasler_, May 25 2018
%o (Python)
%o from sympy import mobius, integer_nthroot
%o def A076467(n):
%o def f(x): return int(n-1+x-integer_nthroot(x,4)[0]+sum(mobius(k)*(integer_nthroot(x,k)[0]+integer_nthroot(x,k<<1)[0]-2) for k in range(3,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 kmax # _Chai Wah Wu_, Aug 14 2024
%Y Subsequence of A036966.
%Y Cf. A001597, A076468, A076469, A076470.
%K nonn,changed
%O 1,2
%A _Robert G. Wilson v_, Oct 14 2002
%E Edited by _Robert Israel_, Sep 30 2015
|
