Year-end appeal: Please make a donation to the OEIS Foundation to support ongoing development and maintenance of the OEIS. We are now in our 61st year, we have over 378,000 sequences, and we’ve reached 11,000 citations (which often say “discovered thanks to the OEIS”).
%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