A064602 Partial sums of A001157: Sum_{j=1..n} sigma_2(j).

1, 6, 16, 37, 63, 113, 163, 248, 339, 469, 591, 801, 971, 1221, 1481, 1822, 2112, 2567, 2929, 3475, 3975, 4585, 5115, 5965, 6616, 7466, 8286, 9336, 10178, 11478, 12440, 13805, 15025, 16475, 17775, 19686, 21056, 22866, 24566, 26776, 28458, 30958

Offset: 1

Author

Labos Elemer, Sep 24 2001

Keywords

• nonn

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.

Crossrefs

Cf. A001157, A064605.

Cf. A064603, A064604, A248076.

Programs

Mathematica

Accumulate@ Array[DivisorSigma[2, #] &, 42] (* Michael De Vlieger, Jan 02 2017 *)

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

Formula

a(n) = a(n-1) + A001157(n) = Sum_{j=1..n} sigma_2(j) where sigma_2(j) = A001157(j).

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