A076467 Perfect powers m^k where m is a positive integer and k >= 3.

1, 8, 16, 27, 32, 64, 81, 125, 128, 216, 243, 256, 343, 512, 625, 729, 1000, 1024, 1296, 1331, 1728, 2048, 2187, 2197, 2401, 2744, 3125, 3375, 4096, 4913, 5832, 6561, 6859, 7776, 8000, 8192, 9261, 10000, 10648, 12167, 13824, 14641, 15625, 16384, 16807

Offset: 1

Author

Robert G. Wilson v, Oct 14 2002

Keywords

• nonn

Comments

If p|n with p prime then p^3|n.

Links

• Alois P. Heinz, Table of n, a(n) for n = 1..10000 (terms n = 1..250 from Reinhard Zumkeller)

Crossrefs

Subsequence of A036966.

Cf. A001597, A076468, A076469, A076470.

Programs

Haskell

a076467 n = a076467_list !! (n-1)
a076467_list = 1 : filter ((> 2) . foldl1 gcd . a124010_row) [2..]
-- Reinhard Zumkeller, Apr 13 2012

Haskell

import qualified Data.Set as Set (null)
import Data.Set (empty, insert, deleteFindMin)
a076467 n = a076467_list !! (n-1)
a076467_list = 1 : f [2..] empty where
   f xs'@(x:xs) s | Set.null s || m > x ^ 3 = f xs $ insert (x ^ 3, x) s
                  | m == x ^ 3  = f xs s
                  | otherwise = m : f xs' (insert (m * b, b) s')
                  where ((m, b), s') = deleteFindMin s
-- Reinhard Zumkeller, Jun 18 2013

Maple

N:= 10^5: # to get all terms <= N
S:= {1, seq(seq(m^k, m = 2 .. floor(N^(1/k))),k=3..ilog2(N))}:
sort(convert(S,list)); # Robert Israel, Sep 30 2015

Mathematica

a = {1}; Do[ If[ Apply[ GCD, Last[ Transpose[ FactorInteger[n]]]] > 2, a = Append[a, n]; Print[n]], {n, 2, 17575}]; a
(* Second program: *)
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 *)

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

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

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

Python

from sympy import mobius, integer_nthroot
def A076467(n):
    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())))
    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 kmax # Chai Wah Wu, Aug 14 2024

Formula

For n > 1: GCD(exponents in prime factorization of a(n)) > 2, cf. A124010. - Reinhard Zumkeller, Apr 13 2012

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

Extensions

Edited by Robert Israel, Sep 30 2015