A025475 1 and the prime powers p^m where m >= 2, thus excluding the primes.

1, 4, 8, 9, 16, 25, 27, 32, 49, 64, 81, 121, 125, 128, 169, 243, 256, 289, 343, 361, 512, 529, 625, 729, 841, 961, 1024, 1331, 1369, 1681, 1849, 2048, 2187, 2197, 2209, 2401, 2809, 3125, 3481, 3721, 4096, 4489, 4913, 5041, 5329, 6241, 6561, 6859, 6889, 7921, 8192

Offset: 1

Author

David W. Wilson

Keywords

• nonn
• easy
• nice

Comments

Also nonprime n such that sigma(n)*phi(n) > (n-1)^2. - Benoit Cloitre, Apr 12 2002

If p is a term of the sequence, then the index n for which a(n) = p is given by n := b(p) := 1 + Sum_{k>=2} PrimePi(p^(1/k)). Here, the sum has floor(log_2(p)) positive terms. For any m > 0, the greatest number n such that a(n) <= m is also given by b(m), thus, b(m) is the number of such prime powers <= m. - Hieronymus Fischer, May 31 2013

That 8 and 9 are the only two consecutive integers in this sequence is known as Catalan's Conjecture and was proved in 2002 by Preda Mihăilescu. - Geoffrey Critzer, Nov 15 2015

Links

• T. D. Noe, Table of n, a(n) for n = 1..10000
• Romeo Meštrović, Generalizations of Carmichael numbers I, arXiv:1305.1867v1 [math.NT], May 4, 2013.
• Preda Mihăilescu, On Catalan's Conjecture, Kuwait Foundation Lecture 30 - April 28, 2003.
• Eric Weisstein's World of Mathematics, Prime Power.

Crossrefs

Subsequence of A000961. - Reinhard Zumkeller, Jun 22 2011

Cf. A001597, A000720, A068982, A136141, A193166.

Differences give A053707.

Cf. A076048 (number of terms < 10^n).

There are four different sequences which may legitimately be called "prime powers": A000961 (p^k, k >= 0), A246655 (p^k, k >= 1), A246547 (p^k, k >= 2), A025475 (p^k, k=0 and k >= 2). When you refer to "prime powers", be sure to specify which of these you mean. Also A001597 is the sequence of nontrivial powers n^k, n >= 1, k >= 2. - N. J. A. Sloane, Mar 24 2018

Programs

Haskell

a025475 n = a025475_list !! (n-1)
a025475_list = filter ((== 0) . a010051) a000961_list
-- Reinhard Zumkeller, Jun 22 2011

Maple

isA025475 := proc(n)
    if n < 1 then
        false;
    elif n = 1 then
        true;
    elif isprime(n) then
        false;
    elif nops(numtheory[factorset](n)) = 1 then
        true;
    else
        false;
    end if;
end proc:
A025475 := proc(n)
    option remember;
    local a;
    if n = 1 then
        1;
    else
        for a from procname(n-1)+1 do
            if isA025475(a) then
                return a;
            end if;
        end do:
    end if;
end proc:
# R. J. Mathar, Jun 06 2013
# alternative:
N:= 10^5: # to get all terms <= N
Primes:= select(isprime, [2,(2*i+1 $ i = 1 .. floor((sqrt(N)-1)/2))]):
sort([1,seq(seq(p^i, i=2..floor(log[p](N))),p=Primes)]); # Robert Israel, Jul 27 2015

Mathematica

A025475 = Select[ Range[ 2, 10000 ], ! PrimeQ[ # ] && Mod[ #, # - EulerPhi[ # ] ] == 0 & ]
A025475 = Sort[ Flatten[ Table[ Prime[n]^i, {n, 1, PrimePi[ Sqrt[10^4]]}, {i, 2, Log[ Prime[n], 10^4]}]]]
{1}~Join~Select[Range[10^4], And[! PrimeQ@ #, PrimePowerQ@ #] &] (* Michael De Vlieger, Jul 04 2016 *)
Join[{1},Select[Range[100000],PrimePowerQ[#]&&!PrimeQ[#]&]] (* Harvey P. Dale, Oct 29 2023 *)

PARI

for(n=1,10000,if(sigma(n)*eulerphi(n)*(1-isprime(n))>(n-1)^2,print1(n,",")))

PARI

is_A025475(n)={ ispower(n,,&p) && isprime(p) || n==1 }  \\ M. F. Hasler, Sep 25 2011

PARI

list(lim)=my(v=List([1]),L=log(lim+.5));forprime(p=2,(lim+.5)^(1/3),for(e=3,L\log(p),listput(v,p^e))); vecsort(concat(Vec(v), apply(n->n^2,primes(primepi(sqrtint(lim\1)))))) \\ Charles R Greathouse IV, Nov 12 2012

PARI

list(lim)=my(v=List([1])); for(m=2,logint(lim\=1,2), forprime(p=2,sqrtnint(lim,m), listput(v, p^m))); Set(v) \\ Charles R Greathouse IV, Aug 26 2015

Python

from sympy import primerange
A025475_list, m = [1], 10*2
m2 = m**2
for p in primerange(1,m):
    a = p**2
    while a < m2:
        A025475_list.append(a)
        a *= p
A025475_list = sorted(A025475_list) # Chai Wah Wu, Sep 08 2014

Python

from sympy import primepi, integer_nthroot
def A025475(n):
    if n==1: return 1
    def f(x): return int(n-2+x-sum(primepi(integer_nthroot(x,k)[0]) 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 kmax # Chai Wah Wu, Aug 13 2024

Formula

The number of terms <= N is O(sqrt(N)*log N). [See Weisstein link] - N. J. A. Sloane, May 27 2022

A005171(a(n))*A010055(a(n)) = 1. - Reinhard Zumkeller, Nov 01 2009

A192280(a(n)) = 0 for n > 1. - Reinhard Zumkeller, Aug 26 2011

A014963(a(n)) - A089026(a(n)) = A014963(a(n)) - 1. - Eric Desbiaux, May 18 2013

From Hieronymus Fischer, May 31 2013: (Start)

The greatest number n such that a(n) <= m is given by 1 + Sum_{k>=2} A000720(floor(m^(1/k))).

Example 1: m = 10^10 ==> n = 10085;

Example 2: m = 10^11 ==> n = 28157;

Example 3: m = 10^12 ==> n = 80071;

Example 4: m = 10^15 ==> n = 1962690. (End)

Sum_{n>=2} 1/a(n) = Sum_{p prime} 1/(p*(p-1)) = A136141. - Amiram Eldar, Oct 11 2020

From Amiram Eldar, Jan 28 2021: (Start)

Product_{n>=2} (1 + 1/a(n)) = Product_{k>=2} zeta(k)/zeta(2*k) = 2.0729553047...

Product_{n>=2} (1 - 1/a(n)) = A068982. (End)

Extensions

Edited by Daniel Forgues, Aug 18 2009