|
%I #25 Jun 12 2023 08:11:23
%S 0,1,5,11,25,36,71,85,145,176,260,287,455,456,649,726,961,970,1376,
%T 1331,1820,1866,2315,2301,3175,2961,3736,3830,4729,4496,5966,5457,
%U 6945,6842,8114,7890,10196,9140,11215,11126,13420,12342,15730,14191,17515,17106,19601
%N Coefficients in expansion of Sum_{n >= 1} x^n/(1-x^n)^4.
%H G. C. Greubel, <a href="/A059358/b059358.txt">Table of n, a(n) for n = 0..1000</a>
%F a(n) = (1/6)*(sigma_3(n) + 3*sigma_2(n) + 2*sigma_1(n)), i.e., this sequence is the inverse Möbius transform of tetrahedral (or pyramidal) numbers: n*(n+1)(n+2)/6 with g.f. 1/(1-x)^4 (cf. A000292). - _Vladeta Jovovic_, Aug 31 2002
%F L.g.f.: -log(Product_{k>=1} (1 - x^k)^((k+1)*(k+2)/6)) = Sum_{n>=1} a(n)*x^n/n. - _Ilya Gutkovskiy_, May 21 2018
%p a:= proc(n) option remember;
%p add(d*(d+1)*(d+2)/6, d=numtheory[divisors](n))
%p end:
%p seq(a(n), n=0..60); # _Alois P. Heinz_, Jun 12 2023
%t With[{nn=50},CoefficientList[Series[Sum[x^n/(1-x^n)^4,{n,nn}],{x,0,nn}],x]] (* _Harvey P. Dale_, May 14 2013 *)
%o (PARI) a(n) = if(n==0, 0, sumdiv(n, d, binomial(d+2, 3))); \\ _Seiichi Manyama_, Apr 19 2021
%Y Cf. A000203, A000292, A002129, A007437, A048272, A073570, A101289.
%K nonn
%O 0,3
%A _N. J. A. Sloane_, Jan 27 2001
|
