# Greetings from The On-Line Encyclopedia of Integer Sequences! http://oeis.org/

Search: id:a084371
Showing 1-1 of 1

%I A084371 #24 Sep 21 2024 14:43:16
%S A084371 1,2,2,3,2,5,3,2,6,7,2,6,3,10,6,11,5,2,6,13,14,10,6,15,3,2,6,17,6,7,
%T A084371 19,14,10,6,21,22,10,2,23,6,5,6,15,26,3,14,10,29,6,30,31,22,6,10,2,33,
%U A084371 15,6,34,35,6,21,11,26,37,14,38,39,14,10,41,6,42,30,43,22,6,10
%N A084371 Squarefree kernels of powerful numbers (A001694).
%H A084371 Amiram Eldar, <a href="/A084371/b084371.txt">Table of n, a(n) for n = 1..10000</a>
%H A084371 Rafael Jakimczuk, <a href="https://doi.org/10.12988/imf.2017.7759">The kernel of powerful numbers</a>, International Mathematical Forum, Vol. 12, No. 15 (2017), pp. 721-730, Remark 2.5, p. 729.
%F A084371 a(n) = A007947(A001694(n)).
%F A084371 From _Amiram Eldar_, May 13 2023: (Start)
%F A084371 Sum_{A001694(k) < x} a(k) = (1/2) * x + o(x) (Jakimczuk, 2017). [corrected Sep 21 2024]
%F A084371 Sum_{k=1..n} a(k) ~ c * n^2, where c = (zeta(3)/zeta(3/2))^2/2 = 0.1058641473... . (End)
%e A084371 A001694(11) = 64 = 2^6 -> a(11) = 2,
%e A084371 A001694(12) = 72 = 2^3 * 3^2 -> a(12) = 2*3 = 6,
%e A084371 A001694(13) = 81 = 3^4 -> a(13) = 3.
%t A084371 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 *)
%o A084371 (PARI) rad(n) = factorback(factorint(n)[, 1]); \\ A007947
%o A084371 lista(nn) = apply(x->rad(x), select(x->ispowerful(x), [1..nn])); \\ _Michel Marcus_, Aug 22 2019
%o A084371 (Python)
%o A084371 from math import prod, isqrt
%o A084371 from sympy import mobius, integer_nthroot, primefactors
%o A084371 def A084371(n):
%o A084371     def squarefreepi(n): return int(sum(mobius(k)*(n//k**2) for k in range(1, isqrt(n)+1)))
%o A084371     def bisection(f,kmin=0,kmax=1):
%o A084371         while f(kmax) > kmax: kmax <<= 1
%o A084371         while kmax-kmin > 1:
%o A084371             kmid = kmax+kmin>>1
%o A084371             if f(kmid) <= kmid:
%o A084371                 kmax = kmid
%o A084371             else:
%o A084371                 kmin = kmid
%o A084371         return kmax
%o A084371     def f(x):
%o A084371         c, l, j = n+x-squarefreepi(integer_nthroot(x,3)[0]), 0, isqrt(x)
%o A084371         while j>1:
%o A084371             k2 = integer_nthroot(x//j**2,3)[0]+1
%o A084371             w = squarefreepi(k2-1)
%o A084371             c -= j*(w-l)
%o A084371             l, j = w, isqrt(x//k2**3)
%o A084371         return c+l
%o A084371     return prod(primefactors(bisection(f,n,n))) # _Chai Wah Wu_, Sep 13 2024
%Y A084371 Cf. A001694, A007947, A090699.
%K A084371 nonn
%O A084371 1,2
%A A084371 _Reinhard Zumkeller_, Jun 23 2003

# Content is available under The OEIS End-User License Agreement: http://oeis.org/LICENSE