OFFSET
1,2
LINKS
Amiram Eldar, Table of n, a(n) for n = 1..10000
Rafael Jakimczuk, The kernel of powerful numbers, International Mathematical Forum, Vol. 12, No. 15 (2017), pp. 721-730, Remark 2.5, p. 729.
FORMULA
From Amiram Eldar, May 13 2023: (Start)
Sum_{A001694(k) < x} a(k) = (1/2) * x + o(x) (Jakimczuk, 2017). [corrected Sep 21 2024]
Sum_{k=1..n} a(k) ~ c * n^2, where c = (zeta(3)/zeta(3/2))^2/2 = 0.1058641473... . (End)
EXAMPLE
MATHEMATICA
s = {1}; Do[f = FactorInteger[n]; If[Min @ f[[;; , 2]] > 1, AppendTo[s, Times @@ f[[;; , 1]]]], {n, 2, 10^4}]; s (* Amiram Eldar, Aug 22 2019 *)
PROG
(PARI) rad(n) = factorback(factorint(n)[, 1]); \\ A007947
lista(nn) = apply(x->rad(x), select(x->ispowerful(x), [1..nn])); \\ Michel Marcus, Aug 22 2019
(Python)
from math import prod, isqrt
from sympy import mobius, integer_nthroot, primefactors
def A084371(n):
def squarefreepi(n): return int(sum(mobius(k)*(n//k**2) for k in range(1, isqrt(n)+1)))
def bisection(f, kmin=0, kmax=1):
while f(kmax) > kmax: kmax <<= 1
while kmax-kmin > 1:
kmid = kmax+kmin>>1
if f(kmid) <= kmid:
kmax = kmid
else:
kmin = kmid
return kmax
def f(x):
c, l, j = n+x-squarefreepi(integer_nthroot(x, 3)[0]), 0, isqrt(x)
while j>1:
k2 = integer_nthroot(x//j**2, 3)[0]+1
w = squarefreepi(k2-1)
c -= j*(w-l)
l, j = w, isqrt(x//k2**3)
return c+l
return prod(primefactors(bisection(f, n, n))) # Chai Wah Wu, Sep 13 2024
CROSSREFS
KEYWORD
nonn,changed
AUTHOR
Reinhard Zumkeller, Jun 23 2003
STATUS
approved