%I #24 Oct 21 2023 16:52:36
%S 1,12,40,121,207,473,649,1142,1611,2401,2853,4647,5285,6879,8759,
%T 11452,12558,16739,18127,23353,27129,31171,33219,43573,47524,53210,
%U 59538,69996,73274,89694,93446,107195,116731,126545,137505,164580,169946,182244,195644,225454
%N a(n) = Sum_{k=1..n} k^2 * floor(n/k)^3.
%H Seiichi Manyama, <a href="/A350124/b350124.txt">Table of n, a(n) for n = 1..10000</a>
%F a(n) = Sum_{k=1..n} k^2 * Sum_{d|k} (d^3 - (d - 1)^3)/d^2.
%F G.f.: (1/(1 - x)) * Sum_{k>=1} (k^3 - (k - 1)^3) * x^k * (1 + x^k)/(1 - x^k)^3.
%F From _Vaclav Kotesovec_, Aug 03 2022: (Start)
%F a(n) = A064602(n) - 3*A143128(n) + 3*A319085(n).
%F a(n) ~ n^3 * (log(n) + 2*gamma + (zeta(3) - 1)/3 - Pi^2/6), where gamma is the Euler-Mascheroni constant A001620. (End)
%t Accumulate[Table[DivisorSigma[2, k] - 3*k*DivisorSigma[1, k] + 3*k^2*DivisorSigma[0, k], {k, 1, 50}]] (* _Vaclav Kotesovec_, Dec 17 2021 *)
%o (PARI) a(n) = sum(k=1, n, k^2*(n\k)^3);
%o (PARI) a(n) = sum(k=1, n, k^2*sumdiv(k, d, (d^3-(d-1)^3)/d^2));
%o (PARI) my(N=66, x='x+O('x^N)); Vec(sum(k=1, N, (k^3-(k-1)^3)*x^k*(1+x^k)/(1-x^k)^3)/(1-x))
%o (Python)
%o from math import isqrt
%o def A350124(n): return (-(s:=isqrt(n))**4*(s+1)*(2*s+1) + sum((q:=n//k)*(k*(3*(k-1))+q*(k*(9*(k-1))+q*(k*(12*k-6)+2)+3)+1) for k in range(1,s+1)))//6 # _Chai Wah Wu_, Oct 21 2023
%Y Cf. A318742, A350108, A350123, A350125.
%Y Cf. A064602, A143128, A319085.
%K nonn
%O 1,2
%A _Seiichi Manyama_, Dec 15 2021
|