OFFSET
1,4
COMMENTS
See A001694 for the corresponding definition of powerful number.
LINKS
Antti Karttunen, Table of n, a(n) for n = 1..16383 (first 1000 terms from T. D. Noe)
Thomas Bloom, Problem 367 and Problem 935, Erdős Problems.
Erdős problems database contributors, Erdős problem database, see nos. 367 and 935.
Antti Karttunen, Data supplement: n, a(n) computed for n = 1..65537
Vaclav Kotesovec, Graph - the asymptotic ratio (100000000 terms).
Christian Krause, LODA, an assembly language, a computational model and a tool for mining integer sequences.
Victor Ufnarovski and Bo Åhlander, How to Differentiate a Number, J. Integer Seq., Vol. 6 (2003), Article 03.3.4.
FORMULA
a(n) = n / A055231(n).
Multiplicative with a(p)=1 and a(p^e)=p^e for e>1. - Vladeta Jovovic, Nov 01 2001
From Antti Karttunen, Nov 22 2017: (Start)
(End)
a(n) = gcd(n, A003415(n)^k), for all k >= 2. [This formula was found in the form k=3 by Christian Krause's LODA miner. See Ufnarovski and Åhlander paper, Theorem 5 on p. 4 for why this holds] - Antti Karttunen, Mar 09 2021
Dirichlet g.f.: zeta(s-1) * Product_{p prime} (1 + 1/p^s - 1/ p^(s-1) + 1/p^(2*s-2) - 1/p^(2*s-1)). - Amiram Eldar, Sep 18 2023
From Vaclav Kotesovec, Apr 09 2025, simplified May 11 2025: (Start)
Dirichlet g.f.: zeta(2*s-2) * Product_{p prime} (1 - 1/p^(3*s-2) + 1/p^(3*s-3) + 1/p^s).
Sum_{k=1..n} a(k) ~ c * n^(3/2) / 3, where c = Product_{p prime} (1 + 2/p^(3/2) - 1/p^(5/2)) = 3.51955505841710664719752940369857817... = A328013. (End)
EXAMPLE
a(40) = 8 since 40 = 2^3 * 5^1 so the powerful part is 2^3 = 8.
MAPLE
A057521 := proc(n)
local a, d, e, p;
a := 1;
for d in ifactors(n)[2] do
e := d[1] ;
p := d[2] ;
if e > 1 then
a := a*p^e ;
end if;
end do:
return a;
end proc: # R. J. Mathar, Jun 09 2016
MATHEMATICA
rad[n_] := Times @@ First /@ FactorInteger[n]; a[n_] := n/Denominator[n/rad[n]^2]; Table[a[n], {n, 1, 97}] (* Jean-François Alcover, Jun 20 2013 *)
f[p_, e_] := If[e > 1, p^e, 1]; a[1] = 1; a[n_] := Times @@ f @@@ FactorInteger[n]; Array[a, 100] (* Amiram Eldar, Sep 21 2020 *)
PROG
(PARI) a(n)=my(f=factor(n)); prod(i=1, #f~, if(f[i, 2]>1, f[i, 1]^f[i, 2], 1)) \\ Charles R Greathouse IV, Aug 13 2013
(PARI) a(n) = my(f=factor(n)); for (i=1, #f~, if (f[i, 2]==1, f[i, 1]=1)); factorback(f); \\ Michel Marcus, Jan 29 2021
(Python)
from sympy import factorint, prod
def a(n): return 1 if n==1 else prod(1 if e==1 else p**e for p, e in factorint(n).items())
print([a(n) for n in range(1, 51)]) # Indranil Ghosh, Jul 19 2017
(Python)
from math import prod
from sympy import factorint
def A057521(n): return n//prod(p for p, e in factorint(n).items() if e == 1) # Chai Wah Wu, Nov 14 2022
CROSSREFS
KEYWORD
nonn,mult,easy
AUTHOR
Henry Bottomley, Sep 01 2000
EXTENSIONS
Edited by Peter Munn, Nov 14 2025
STATUS
approved