[go: up one dir, main page]

login
A025475
1 and the prime powers p^m where m >= 2, thus excluding the primes.
211
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,2
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
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.
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)
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 *)
PROG
(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
(Haskell)
a025475 n = a025475_list !! (n-1)
a025475_list = filter ((== 0) . a010051) a000961_list
-- Reinhard Zumkeller, Jun 22 2011
(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
CROSSREFS
Subsequence of A000961. - Reinhard Zumkeller, Jun 22 2011
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
Sequence in context: A319163 A134611 A134612 * A246547 A195942 A125643
KEYWORD
nonn,easy,nice
EXTENSIONS
Edited by Daniel Forgues, Aug 18 2009
STATUS
approved