login
a(n) = Sum_{k=1..n} k^n * sigma_0(k).
4

%I #25 Jan 21 2024 18:10:13

%S 1,9,71,963,9873,231749,2976863,86348423,1824883450,55584932826,

%T 1104642697680,64932555347084,1366828157222090,61273696016238014,

%U 2581786206601959958,129797968403021602450,3678372903755436314440,295835829367866540495396

%N a(n) = Sum_{k=1..n} k^n * sigma_0(k).

%H Robert Israel, <a href="/A356239/b356239.txt">Table of n, a(n) for n = 1..384</a>

%F a(n) = Sum_{k=1..n} k^n * Sum_{j=1..floor(n/k)} j^n.

%p f:= proc(n) local k; add(k^n * numtheory:-tau(k),k=1..n) end proc:

%p map(f, [$1..30]); # _Robert Israel_, Jan 21 2024

%t a[n_] := Sum[k^n * DivisorSigma[0, k], {k, 1, n}]; Array[a, 18] (* _Amiram Eldar_, Jul 30 2022 *)

%o (PARI) a(n) = sum(k=1, n, k^n*sigma(k, 0));

%o (PARI) a(n) = sum(k=1, n, k^n*sum(j=1, n\k, j^n));

%o (Python)

%o from math import isqrt

%o from sympy import bernoulli

%o def A356239(n): return (-(bernoulli(n+1, (s:=isqrt(n))+1)-(b:=bernoulli(n+1)))**2//(n+1) + sum(k**n*(bernoulli(n+1, n//k+1)-b)<<1 for k in range(1,s+1)))//(n+1) # _Chai Wah Wu_, Oct 21 2023

%Y Cf. A000005, A319194, A356129, A356243.

%K nonn

%O 1,2

%A _Seiichi Manyama_, Jul 30 2022