Reminder: The OEIS is hiring a new managing editor, and the application deadline is January 26.
%I #15 May 17 2022 17:47:17
%S 1,0,1,0,-1,2,0,4,-8,6,0,-36,72,-60,24,0,600,-1200,1020,-480,120,0,
%T -16584,33168,-28320,13776,-4200,720,0,705600,-1411200,1206240,
%U -591360,188160,-40320,5040,0,-43751232,87502464,-74813760,36747648,-11813760,2661120,-423360,40320
%N Triangle read by rows. The Faulhaber numbers. F(0, k) = 1 and otherwise F(n, k) = (n + 1)!*(-1)^(k+1)*Sum_{j=0..floor((k-1)/2)} C(2*k-2*j, k+1)*C(2*n+1, 2*j+1) * Bernoulli(2*n-2*j) / (k - j).
%C I. Gessel and X. Viennot call the rational numbers F(n, k)/(n + 1)! 'Faulhaber numbers'. However, for our purposes it is more convenient to define the integers F(n, k). For the Faulhaber polynomials see A335951/A335952.
%C Let S(r, m) = Sum_{k=0..m} k^r, with 0^0 = 1 and S(0, m) = m + 1. Faulhaber's theorem (the sums of powers formula) is:
%C S(2*n+1, m) = (1/(n+1)!)*(1/2)*Sum_{k=0..n} F(n, k)*(m*(m + 1))^(k + 1).
%C Gessel and Viennot give two combinatorial interpretations for the Faulhaber numbers, for this see A354043.
%H I. M. Gessel and X. G. Viennot, <a href="https://people.brandeis.edu/~gessel/homepage/papers/pp.pdf">Determinants, Paths, and Plane Partitions</a>, 1989 preprint.
%F F(n,1) = (2*n +1)*Bernoulli(2*n)*(n+1)! for n >= 1.
%F F(n,2) = -(4*n+2)*Bernoulli(2*n)*(n+1)! for n >= 2.
%F F(n,3) = ((10*n+5)*Bernoulli(2*n) + binomial(2*n+1,3)*Bernoulli(2*n-2)/2)*(n+1)! for n >= 3.
%e Triangle starts:
%e 0: 1
%e 1: 0, 1
%e 2: 0, -1, 2
%e 3: 0, 4, -8, 6
%e 4: 0, -36, 72, -60, 24
%e 5: 0, 600, -1200, 1020, -480, 120
%e 6: 0, -16584, 33168, -28320, 13776, -4200, 720
%e 7: 0, 705600, -1411200, 1206240, -591360, 188160, -40320, 5040
%e 8: 0, -43751232, 87502464, -74813760, 36747648, -11813760, 2661120, -423360, 40320
%e .
%e Let n = 4 and m = 3, then S(2*n + 1, m) = S(9, 3) = 20196. Faulhaber's formula gives this as (0*12 + (-36)*144 + 72*1728 + (-60)*20736 + 24*248832) / (2*120).
%p F := (n, k) -> ifelse(n = 0, 1, (n + 1)!*(-1)^(k + 1)*add(binomial(2*k - 2*j, k + 1)*binomial(2*n + 1, 2*j + 1)*bernoulli(2*n - 2*j) / (k - j), j = 0..(k - 1)/2)): for n from 0 to 8 do seq(F(n, k), k = 0..n) od;
%Y Cf. A335951/A335952, A000367/A002445, A354043, A263445.
%K sign,tabl
%O 0,6
%A _Peter Luschny_, May 17 2022