login
Numbers whose prime factors are raised to the fourth power.
20

%I #41 Aug 19 2024 13:16:44

%S 16,81,625,1296,2401,10000,14641,28561,38416,50625,83521,130321,

%T 194481,234256,279841,456976,707281,810000,923521,1185921,1336336,

%U 1500625,1874161,2085136,2313441,2825761,3111696,3418801,4477456,4879681,6765201,7890481

%N Numbers whose prime factors are raised to the fourth power.

%C This is essentially A005117 (the squarefree numbers) raised to the fourth power. - _T. D. Noe_, Mar 13 2013

%C All positive integers have a unique factorization into powers of squarefree numbers with distinct exponents that are powers of two. So every positive number is a product of at most one squarefree number (A005117), at most one square of a squarefree number (A062503), at most one 4th power of a squarefree number (term of this sequence), at most one 8th power of a squarefree number, and so on. - _Peter Munn_, Mar 12 2020

%H T. D. Noe, <a href="/A113849/b113849.txt">Table of n, a(n) for n = 1..10000</a>

%F From _Peter Munn_, Oct 31 2019: (Start)

%F a(n) = A005117(n+1)^4.

%F {a(n)} = {A225546(A000351(n)) : n >= 0} \ {1}, where {a(n)} denotes the set of integers in the sequence.

%F (End)

%F Sum_{k>=1} 1/a(k) = zeta(4)/zeta(8) - 1 = 105/Pi^4 - 1. - _Amiram Eldar_, May 22 2020

%e 1296 = 16*81 = 2^4*3^4 so the prime factors of 1296, 2 and 3, are raised to the fourth power.

%t Select[ Range@50^4, Union[Last /@ FactorInteger@# ] == {4} &] (* _Robert G. Wilson v_, Jan 26 2006 *)

%t nn = 50; t = Select[Range[2, nn], Union[Transpose[FactorInteger[#]][[2]]] == {1} &]; t^4 (* _T. D. Noe_, Mar 13 2013 *)

%t Rest[Select[Range[100], SquareFreeQ]^4] (* _Vaclav Kotesovec_, May 22 2020 *)

%o (PARI) allpwrfact(n,p) = { local(x,j,ln,y,flag); for(x=4,n, y=Vec(factor(x)); ln = length(y[1]); flag=0; for(j=1,ln, if(y[2][j]==p,flag++); ); if(flag==ln,print1(x",")); ) } \\ All prime factors are raised to the power p

%o (Python)

%o from math import isqrt

%o from sympy import mobius

%o def A113849(n):

%o def f(x): return n+x-sum(mobius(k)*(x//k**2) for k in range(1, isqrt(x)+1))

%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**4 # _Chai Wah Wu_, Aug 19 2024

%Y Proper subset of A000583.

%Y Other powers of squarefree numbers: A005117(1), A062503(2), A062838(3), A113850(5), A113851(6), A113852(7), A072774(all).

%Y Cf. A000351, A225546.

%K easy,nonn

%O 1,1

%A _Cino Hilliard_, Jan 25 2006

%E More terms from _Robert G. Wilson v_, Jan 26 2006