A367326 - OEIS

%I #35 Nov 14 2023 17:55:30

%S 0,0,1,4,13,16,50,36,101,100,170,100,402,144,362,440,629,256,995,324,

%T 1266,920,962,484,2578,976,1370,1576,2618,784,4180,900,3477,2408,2402,

%U 2840,7023,1296,3026,3416,7810,1600,8548,1764,6786,7496,4490,2116,14546

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

%F a(n) = n^2*A000005(n) - 2*n*A000203(n) + A001157(n) for n >= 1. - _Chai Wah Wu_, Nov 14 2023

%p divides := (k, n) -> k = n or (k > 0 and irem(n, k) = 0):

%p a := n -> local d; add(`if`(divides(n - d, n), d^2, 0), d = 0..n-1):

%p seq(a(n), n = 0..61);

%t a[0] = 0; a[n_] := DivisorSum[n, (n - #)^2 &]; Array[a, 50, 0]

%t (* _Amiram Eldar_, Nov 14 2023 *)

%t a[n_] := Sum[If[Divisible[n, n - k], k^2, 0], {k, 0, n - 1}]

%t Table[a[n], {n, 0, 48}]  (* _Peter Luschny_, Nov 14 2023 *)

%o (SageMath)

%o def A367326(n): return sum(k^2 for k in (0..n) if (n-k).divides(n))

%o print([A367326(n) for n in range(50)])

%o (Julia)

%o using AbstractAlgebra

%o function A367326(n)

%o     sum(k^2 for k in 0:n if is_divisible_by(n, n - k))

%o end

%o [A367326(n) for n in 0:48] |> println

%o (PARI)

%o a(n) = if(n == 0, 0, sumdiv(n, d, (n-d)^2)); \\ _Michel Marcus_, Nov 14 2023

%o (Python)

%o from math import prod

%o from sympy import factorint

%o def A367326(n):

%o     f = factorint(n).items()

%o     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

%o # _Chai Wah Wu_, Nov 14 2023

%Y Cf. A094471, A001157, A007434.

%K nonn

%O 0,4

%A _Peter Luschny_, Nov 14 2023