A367326 a(n) = Sum_{(n - k)|n} k^2.

0, 0, 1, 4, 13, 16, 50, 36, 101, 100, 170, 100, 402, 144, 362, 440, 629, 256, 995, 324, 1266, 920, 962, 484, 2578, 976, 1370, 1576, 2618, 784, 4180, 900, 3477, 2408, 2402, 2840, 7023, 1296, 3026, 3416, 7810, 1600, 8548, 1764, 6786, 7496, 4490, 2116, 14546

Offset: 0

Author

Peter Luschny, Nov 14 2023

Keywords

• nonn

Crossrefs

Cf. A094471, A001157, A007434.

Programs

Julia

using AbstractAlgebra
function A367326(n)
    sum(k^2 for k in 0:n if is_divisible_by(n, n - k))
end
[A367326(n) for n in 0:48] |> println

Maple

divides := (k, n) -> k = n or (k > 0 and irem(n, k) = 0):
a := n -> local d; add(`if`(divides(n - d, n), d^2, 0), d = 0..n-1):
seq(a(n), n = 0..61);

Mathematica

a[0] = 0; a[n_] := DivisorSum[n, (n - #)^2 &]; Array[a, 50, 0]
(* Amiram Eldar, Nov 14 2023 *)
a[n_] := Sum[If[Divisible[n, n - k], k^2, 0], {k, 0, n - 1}]
Table[a[n], {n, 0, 48}]  (* Peter Luschny, Nov 14 2023 *)

PARI

a(n) = if(n == 0, 0, sumdiv(n, d, (n-d)^2)); \\ Michel Marcus, Nov 14 2023

Python

from math import prod
from sympy import factorint
def A367326(n):
    f = factorint(n).items()
    return n*(n*prod(e+1 for p,e in f) - (prod((p**(e+1)-1)//(p-1) for p,e in f)<<1))+prod((p**(e+1<<1)-1)//(p**2-1) for p,e in f) if n else 0
# Chai Wah Wu, Nov 14 2023

SageMath

def A367326(n): return sum(k^2 for k in (0..n) if (n-k).divides(n))
print([A367326(n) for n in range(50)])

Formula

a(n) = n^2*A000005(n) - 2*n*A000203(n) + A001157(n) for n >= 1. - Chai Wah Wu, Nov 14 2023