OFFSET
1,1
COMMENTS
This is essentially A005117 (the squarefree numbers) raised to the fourth power. - T. D. Noe, Mar 13 2013
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
LINKS
T. D. Noe, Table of n, a(n) for n = 1..10000
FORMULA
From Peter Munn, Oct 31 2019: (Start)
a(n) = A005117(n+1)^4.
{a(n)} = {A225546(A000351(n)) : n >= 0} \ {1}, where {a(n)} denotes the set of integers in the sequence.
(End)
Sum_{k>=1} 1/a(k) = zeta(4)/zeta(8) - 1 = 105/Pi^4 - 1. - Amiram Eldar, May 22 2020
EXAMPLE
1296 = 16*81 = 2^4*3^4 so the prime factors of 1296, 2 and 3, are raised to the fourth power.
MATHEMATICA
Select[ Range@50^4, Union[Last /@ FactorInteger@# ] == {4} &] (* Robert G. Wilson v, Jan 26 2006 *)
nn = 50; t = Select[Range[2, nn], Union[Transpose[FactorInteger[#]][[2]]] == {1} &]; t^4 (* T. D. Noe, Mar 13 2013 *)
Rest[Select[Range[100], SquareFreeQ]^4] (* Vaclav Kotesovec, May 22 2020 *)
PROG
(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
(Python)
from math import isqrt
from sympy import mobius
def A113849(n):
def f(x): return n+x-sum(mobius(k)*(x//k**2) for k in range(1, isqrt(x)+1))
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**4 # Chai Wah Wu, Aug 19 2024
CROSSREFS
KEYWORD
easy,nonn
AUTHOR
Cino Hilliard, Jan 25 2006
EXTENSIONS
More terms from Robert G. Wilson v, Jan 26 2006
STATUS
approved