%I #39 Aug 12 2024 16:35:43
%S 1,3,4,3,6,12,8,4,18,12,12,14,24,24,18,12,20,18,32,36,24,6,42,24,30,
%T 72,32,48,54,48,12,38,60,56,42,96,44,36,24,72,48,8,18,72,42,54,72,80,
%U 90,60,72,62,96,32,84,144,68,54,96,144,72,74,114,24,60,96,168
%N The sum of squarefree divisors of the cubefree numbers.
%C The number of squarefree divisors of the n-th cubefree number is A366536(n).
%H Amiram Eldar, <a href="/A369889/b369889.txt">Table of n, a(n) for n = 1..10000</a>
%F a(n) = A048250(A004709(n)).
%F Sum_{j=1..n} a(j) ~ c * n^2, where c = zeta(3)^2/(2*zeta(5)) = 0.6967413068... .
%F In general, the formula holds for the sum of squarefree divisors of the k-free numbers with c = zeta(k)^2/(2*zeta(2*k-1))..., for k >= 2.
%t f[p_, e_] := p + 1; s[1] = 1; s[n_] := Times @@ f @@@ FactorInteger[n]; cubefreeQ[n_] := Max[FactorInteger[n][[;; , 2]]] < 3; s /@ Select[Range[100], cubefreeQ]
%t (* or *)
%t f[p_, e_] := If[e > 2, 0, p + 1]; s[1] = 1; s[n_] := Times @@ f @@@ FactorInteger[n]; Select[Array[s, 100], # > 0 &]
%o (PARI) lista(kmax) = {my(f, s, p, e); for(k = 1, kmax, f = factor(k); s = prod(i = 1, #f~, p = f[i,1]; e = f[i,2]; if(e < 3, p + 1, 0)); if(s > 0, print1(s, ", ")));}
%o (Python)
%o from math import prod
%o from sympy import mobius, integer_nthroot, primefactors
%o def A369889(n):
%o def f(x): return n+x-sum(mobius(k)*(x//k**3) for k in range(1, integer_nthroot(x,3)[0]+1))
%o m, k = n, f(n)
%o while m != k:
%o m, k = k, f(k)
%o return prod(p+1 for p in primefactors(m)) # _Chai Wah Wu_, Aug 12 2024
%Y Cf. A004709, A048250, A062822, A366440, A366536, A366537.
%Y Cf. A002117, A013663.
%K nonn,easy
%O 1,2
%A _Amiram Eldar_, Feb 15 2024