%I #86 Mar 18 2022 22:55:41
%S 1,5,36,354,4425,67171,1200304,24684612,574304985,14914341925,
%T 427675990236,13421957361110,457593884876401,16841089312342855,
%U 665478473553144000,28101527071305611528,1262899292504270591313,60182438244917445266889,3031284048960901518840700
%N a(n) = Sum_{k=1..n} k^n.
%C From _Alexander Adamchuk_, Jul 21 2006: (Start)
%C p^(3k - 1) divides a(p^k) for prime p > 2 and k = 1, 2, 3, 4, ... or p^2 divides a(p) for prime p > 2. p^5 divides a(p^2) for prime p > 2. p^8 divides a(p^3) for prime p > 2. p^11 divides a(p^4) for prime p > 2.
%C p^2 divides a(2p) for prime p > 3. p^3 divides a(3p) for prime p > 2. p^2 divides a(4p) for prime p > 5. p^3 divides a(5p) for prime p > 3. p^2 divides a(6p) for prime p > 7.
%C p divides a(2p - 1) for all prime p. p^3 divides a(2p^2 - 1) for all prime p. p^5 divides a(2p^3 - 1) for all prime p.
%C p divides a((p - 1)/2) for p = 5, 13, 17, 29, 37, 41, 53, 61, ... = A002144 Pythagorean primes: primes of form 4n + 1.
%C (End)
%C If p prime then a(p-1) == -1 (mod p) [see De Koninck & Mercier reference]. Example: for p = 7, a(6) = 67171 = 7 * 9596 - 1. - _Bernard Schott_, Mar 06 2020
%D J.-M. De Koninck et A. Mercier, 1001 Problèmes en Théorie Classique des Nombres, Problème 327 pp. 48-200, Ellipses, Paris (2004).
%H T. D. Noe, <a href="/A031971/b031971.txt">Table of n, a(n) for n = 1..100</a>
%H J. M. Grau, A. M. Oller-Marcen, and J. Sondow, <a href="http://arxiv.org/abs/1309.7941">On the congruence 1^m + 2^m + ... + m^m == n (mod m) with n|m</a>, arXiv 1309.7941 [math.NT], 2013-2014.
%F a(n) is asymptotic to (e/(e - 1))*n^n. - _Benoit Cloitre_, Dec 17 2003
%F a(n) = zeta(-n) - zeta(-n, n + 1), where zeta(s) is the Riemann zeta function and zeta(s, a) is the Hurwitz zeta function, a generalization of the Riemann zeta function. - _Alexander Adamchuk_, Jul 21 2006
%F a(n) == 1 (mod n) <==> n is in A014117 = 1, 2, 6, 42, 1806 (see the link "On the congruence ..."). - _Jonathan Sondow_, Oct 18 2013
%F a(A054377(n)) = A233045(n). - _Jonathan Sondow_, Dec 11 2013
%F a(n) = n! * [x^n] exp(x)*(exp(n*x) - 1)/(exp(x) - 1). - _Ilya Gutkovskiy_, Apr 07 2018
%F a(n) ~ ((e*n+1)/((e-1)*(n+1))) * n^n. - _N. J. A. Sloane_, Oct 13 2018, based on email from Claude F. Leibovici who claims this is slightly better than Cloitre's version when n is small.
%p a := n->sum('i'^n,'i'=1..n);
%p # alternative
%p A031971 := proc(n)
%p (bernoulli(n+1,n+1)-bernoulli(n+1))/(n+1) ;
%p end proc: # _R. J. Mathar_, May 10 2013
%t Table[Zeta[-n] - Zeta[-n, n + 1], {n, 25}] (* _Alexander Adamchuk_, Jul 21 2006 *)
%t Table[Total[Range[n]^n], {n,25}] (* _T. D. Noe_, Apr 19 2011 *)
%t Table[HarmonicNumber[n, -n], {n, 1, 25}] (* _Jean-François Alcover_, Apr 09 2015 *)
%o (Magma) [&+[(k)^n: k in [0..n]]: n in [1..30]]; // _Vincenzo Librandi_, Apr 18 2011
%o (Haskell)
%o a031971 = sum . a089072_row -- _Reinhard Zumkeller_, Mar 18 2013
%o (PARI) a(n)=sum(k=1,n,k^n) \\ _Charles R Greathouse IV_, Jun 05 2015
%o (Python)
%o from sympy import harmonic
%o def A031971(n):
%o return harmonic(n,-n) # _Chai Wah Wu_, Feb 15 2020
%Y A diagonal of array A103438.
%Y Cf. A002144, A014117, A089072, A229303, A256016, A261490.
%Y For a(n) mod n see A182398.
%K nonn,nice,easy
%O 1,2
%A Chris du Feu (chris(AT)beckingham0.demon.co.uk)