OFFSET
1,2
COMMENTS
In general, for m >= 0 and j >= 0, Sum_{k=1..n} k^m * sigma_j(k) = Sum_{k=1..s} (k^m * F_{m+j}(floor(n/k)) + k^(m+j) * F_m(floor(n/k))) - F_{m+j}(s) * F_m(s), where s = floor(sqrt(n)) and F_m(x) are the Faulhaber polynomials defined as F_m(x) = (Bernoulli(m+1, x+1) - Bernoulli(m+1, 1)) / (m+1). - Daniel Suteu, Nov 27 2020
LINKS
Seiichi Manyama, Table of n, a(n) for n = 1..10000 (terms 1..1000 from T. D. Noe)
Project Euler, Problem 401: Sum of squares of divisors, 2012.
FORMULA
a(n) = Sum_{i=1..n} i^2 * floor(n/i). - Enrique Pérez Herrero, Sep 15 2012
G.f.: (1/(1 - x))*Sum_{k>=1} k^2*x^k/(1 - x^k). - Ilya Gutkovskiy, Jan 02 2017
a(n) ~ zeta(3) * n^3 / 3. - Vaclav Kotesovec, Sep 02 2018
a(n) = Sum_{k=1..s} (A000330(floor(n/k)) + k^2*floor(n/k)) - s*A000330(s), where s = floor(sqrt(n)). - Daniel Suteu, Nov 26 2020
MATHEMATICA
Accumulate@ Array[DivisorSigma[2, #] &, 42] (* Michael De Vlieger, Jan 02 2017 *)
PROG
(PARI) a(n) = sum(j=1, n, sigma(j, 2)); \\ Michel Marcus, Dec 15 2013
(PARI) f(n) = n*(n+1)*(2*n+1)/6; \\ A000330
a(n) = my(s=sqrtint(n)); sum(k=1, s, f(n\k) + k^2*(n\k)) - s*f(s); \\ Daniel Suteu, Nov 26 2020
(Python)
from math import isqrt
def f(n): return n*(n+1)*(2*n+1)//6
def a(n):
s = isqrt(n)
return sum(f(n//k) + k*k*(n//k) for k in range(1, s+1)) - s*f(s)
print([a(k) for k in range(1, 43)]) # Michael S. Branicky, Oct 01 2022 after Daniel Suteu
CROSSREFS
KEYWORD
nonn
AUTHOR
Labos Elemer, Sep 24 2001
STATUS
approved