login
a(n) = Sum_{k=1..n} k^2*sigma_2(k), where sigma_2 is A001157.
2

%I #30 Oct 22 2023 00:46:19

%S 1,21,111,447,1097,2897,5347,10787,18158,31158,45920,76160,104890,

%T 153890,212390,299686,383496,530916,661598,879998,1100498,1395738,

%U 1676108,2165708,2572583,3147183,3744963,4568163,5276285,6446285,7370767,8768527,10097107

%N a(n) = Sum_{k=1..n} k^2*sigma_2(k), where sigma_2 is A001157.

%H Wikipedia, <a href="https://en.wikipedia.org/wiki/Faulhaber%27s_formula">Faulhaber's formula</a>.

%F a(n) = Sum_{k=1..n} k^4 * A000330(floor(n/k)).

%F a(n) ~ (zeta(3)/5) * n^5. - _Amiram Eldar_, Oct 20 2023

%t Accumulate[Table[n^2*DivisorSigma[2, n], {n, 1, 33}]] (* _Amiram Eldar_, Oct 20 2023 *)

%o (PARI) f(n, m) = (subst(bernpol(m+1, x), x, n+1)-subst(bernpol(m+1, x), x, 0))/(m+1);

%o a(n, s=2, t=2) = sum(k=1, n, k^(s+t)*f(n\k, s));

%o (Python)

%o def A364268(n): return sum(k**4*(m:=n//k)*(m+1)*((m<<1)+1)//6 for k in range(1,n+1)) # _Chai Wah Wu_, Oct 20 2023

%o (Python)

%o from math import isqrt

%o def A364268(n): return (((s:=isqrt(n))*(s+1)*(2*s+1))**2*(1-3*s*(s+1))//6 + sum((q:=n//k)*(q+1)*(2*q+1)*k**2*(5*k**2+3*q*(q+1)-1) for k in range(1,s+1)))//30 # _Chai Wah Wu_, Oct 21 2023

%Y Cf. A356125, A364269.

%Y Cf. A000330, A001157.

%K nonn,easy

%O 1,2

%A _Seiichi Manyama_, Oct 20 2023