This is a front-end for the Online Encyclopedia of Integer Sequences, made by Christian Perfect. The idea is to provide OEIS entries in non-ancient HTML, and then to think about how they're presented visually. The source code is on GitHub.
%I A086787 #40 Mar 20 2023 16:50:50
%S A086787 1,8,56,494,5699,82200,1419760,28501116,651233661,16676686696,
%T A086787 472883843992,14705395791306,497538872883727,18193397941038736,
%U A086787 714950006521386976,30046260016074301944,1344648068888240941017
%N A086787 a(n) = Sum_{i=1..n} ( Sum_{j=1..n} i^j ).
%C A086787 p divides a(p+1) for all prime p except 3. p^2 divides a(p+1) for prime p in A123374.
%C A086787 2 divides a(n) for n = {2, 3, 4, 6, 7, 8, 10, 11, 12, 14, 15, 16, 18, 19, 20, 22, 23, 24, 26, 27, 28, 30, 31, 32, 34, 35, 36, 38, 39, 40, 42, 43, 44, 46, 47, 48, 50, ...}.
%C A086787 2^2 divides a(n) for n = {2, 3, 6, 7, 8, 10, 11, 14, 15, 16, 18, 19, 22, 23, 24, 26, 27, 30, 31, 32, 34, 35, 38, 39, 40, 42, 43, 46, 47, 48, 50, ...}.
%C A086787 2^3 divides a(n) for n = {2, 3, 6, 7, 10, 11, 14, 15, 16, 18, 19, 22, 23, 26, 27, 30, 31, 32, 34, 35, 38, 39, 42, 43, 46, 47, 48, 50, ...}.
%C A086787 2^4 divides a(n) for n = {7, 14, 15, 18, 23, 30, 31, 32, 34, 39, 46, 47, 50, ...}.
%C A086787 2^5 divides a(n) for n = {15, 30, 31, 34, 47, 62, 63, 64, 66, 79, 94, 95, 98, ...}.
%C A086787 2^6 divides a(n) for n = {31, 62, 63, 66, 95, ...}.
%C A086787 2^7 divides a(n) for n = {63, 126, 127, 130, ...}.
%C A086787 It appears that for k > 2 the least few n such that a(n) is divisible by 2^(k+1) are n = {(2^k-1), 2*(2^k-1), 2*(2^k-1)+1, 2*(2^k-1)+3, 3*(2^k-1)+2, 4*(2^k-1)+2, 4*(2^k-1)+3, 4*(2^k-1)+4, 4*(2^k-1)+6, 5*(2^k-1)+4, 6*(2^k-1)+4, 6*(2^k-1)+5, 6*(2^k-1)+8, ...}. - _Alexander Adamchuk_, Oct 08 2006
%C A086787 Numbers n that divide a(n) are listed in A014741. - _Alexander Adamchuk_, Nov 03 2006
%H A086787 G. C. Greubel, <a href="/A086787/b086787.txt">Table of n, a(n) for n = 1..385</a>
%F A086787 1 - Psi(n) - gamma + Sum_{i=2..n} (i^(n+1)/(i-1)), where Psi(n) is the digamma function and gamma is Euler's constant.
%F A086787 a(n) = Sum[ i^j, {i,1,n}, {j,1,n} ] = n + Sum[ i*(i^n - 1)/(i - 1), {i,2,n} ]. - _Alexander Adamchuk_, Nov 03 2006
%F A086787 a(n) = Sum_{k=1..n} (B(k+1, n+1) - B(k+1, 1))/(k+1), where B(n, x) are the Bernoulli polynomials. - _Daniel Suteu_, Jun 25 2018
%F A086787 a(n) ~ c * n^n, where c = 1 / (1 - exp(-1)) = A185393 = 1.58197670686932642438... - _Vaclav Kotesovec_, Nov 06 2021
%e A086787 a(2) = 8 = 1 + 1 + 2 + 4 = 1^1 + 1^2 + 2^1 + 2^2.
%p A086787 seq(1-Psi(n)-gamma+sum(i^(n+1)/(i-1),i = 2 .. n),n=1..20);
%t A086787 Table[Sum[i^j,{i,1,n},{j,1,n}],{n,1,24}] (* _Alexander Adamchuk_, Oct 08 2006 *)
%t A086787 Table[ n + Sum[ i*(i^n-1)/(i-1), {i,2,n} ], {n,1,17} ] (* _Alexander Adamchuk_, Nov 03 2006 *)
%o A086787 (PARI) a(n)=sum(i=1,n,sum(j=1,n,i^j)) \\ _Charles R Greathouse IV_, Jul 19 2013
%o A086787 (PARI) a(n)=round(1-psi(n)-Euler+sum(i=2,n,i^(n+1)/(i-1))) \\ _Charles R Greathouse IV_, Jul 19 2013
%o A086787 (Python)
%o A086787 def A086787(n): return sum(i**j for i in range(1,n+1) for j in range(1,n+1)) # _Chai Wah Wu_, Jan 08 2022
%o A086787 (Python)
%o A086787 from sympy import digamma, EulerGamma
%o A086787 from fractions import Fraction
%o A086787 def A086787(n): return 1-digamma(n)-EulerGamma + sum(Fraction(i**(n+1),i-1) for i in range(2,n+1)) # _Chai Wah Wu_, Jan 08 2022
%Y A086787 Cf. A000295, A014741, A215084, A349117.
%K A086787 nonn
%O A086787 1,2
%A A086787 Klaus Strassburger (strass(AT)ddfi.uni-duesseldorf.de), Aug 04 2003
%E A086787 Edited by _Max Alekseyev_, Jan 29 2012