%I #19 Sep 21 2024 14:48:30
%S 1,15,31,40,63,127,121,156,255,600,364,511,400,1240,1023,781,1815,
%T 1093,2520,2340,2047,3751,1464,5080,5460,4836,4095,3280,2380,2801,
%U 7623,6000,3906,6240,10200,11284,9828,8191,5220,11715,15367,12400,16395,9841,7240,20440
%N a(n) = A000203(A036966(n)), the sum of divisors of the n-th cubefull number A036966(n).
%H Amiram Eldar, <a href="/A362986/b362986.txt">Table of n, a(n) for n = 1..10000</a>
%H Rafael Jakimczuk and Matilde Lalín, <a href="https://doi.org/10.7546/nntdm.2022.28.4.617-634">Asymptotics of sums of divisor functions over sequences with restricted factorization structure</a>, Notes on Number Theory and Discrete Mathematics, Vol. 28, No. 4 (2022), pp. 617-634, eq. (7).
%F Sum_{A036966(k) < x} a(k) = c * x^(4/3) + O(x^(113/96 + eps)), where c = A362985 * A362974 / 4 = 2.8912833599... (Jakimczuk and Lalín, 2022). [corrected Sep 21 2024]
%F Sum_{k=1..n} a(k) ~ c * n^4, where c = A362985 / (4 * A362974^3) = 0.006135085083... .
%t DivisorSigma[1, Select[Range[10^4], # == 1 || Min[FactorInteger[#][[;; , 2]]] > 2 &]]
%o (PARI) lista(kmax) = for(k = 1, kmax, if(k==1 || vecmin(factor(k)[, 2]) > 2, print1(sigma(k), ", ")));
%o (Python)
%o from itertools import count, islice
%o from math import prod
%o from sympy import factorint
%o def A362986_gen(): # generator of terms
%o for n in count(1):
%o f = factorint(n)
%o if all(e>2 for e in f.values()):
%o yield prod((p**(e+1)-1)//(p-1) for p,e in f.items())
%o A362986_list = list(islice(A362986_gen(),20)) # _Chai Wah Wu_, May 21 2023
%Y Cf. A000203, A036966, A180114, A362974, A362985.
%K nonn
%O 1,2
%A _Amiram Eldar_, May 12 2023