OFFSET
0,4
FORMULA
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 *)
PROG
(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)])
(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
(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
CROSSREFS
KEYWORD
nonn
AUTHOR
Peter Luschny, Nov 14 2023
STATUS
approved