A084371 Squarefree kernels of powerful numbers (A001694).

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, 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, 15, 6, 34, 35, 6, 21, 11, 26, 37, 14, 38, 39, 14, 10, 41, 6, 42, 30, 43, 22, 6, 10

Offset: 1

Author

Reinhard Zumkeller, Jun 23 2003

Keywords

• nonn

Examples

			A001694(11) = 64 = 2^6 -> a(11) = 2,
A001694(12) = 72 = 2^3 * 3^2 -> a(12) = 2*3 = 6,
A001694(13) = 81 = 3^4 -> a(13) = 3.
		

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.

Crossrefs

Cf. A001694, A007947, A090699.

Programs

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 *)

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

Formula

a(n) = A007947(A001694(n)).

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)