Named after Johann Faulhaber, the rows of Faulhaber's triangle are the coefficients of polynomials that represent sums of integer powers, which are extracted from Faulhaber's formula:








∑

k
=
1


n



k

p


=


1

p
+
1




∑

j
=
0


p





(



p
+
1

j


)




B

j



n

p
+
1
−
j




{\displaystyle \sum _{k=1}^{n}k^{p}={1 \over p+1}\sum _{j=0}^{p}{p+1 \choose j}B_{j}n^{p+1-j}}




where 




B

n




{\displaystyle B_{n}}

 is the nth-Bernoulli number.



The first 5 rows of Faulhaber's triangle, are:

    1
  1/2  1/2
  1/6  1/2  1/3
    0  1/4  1/2  1/4
-1/30    0  1/3  1/2  1/5



Using the third row of the triangle, we have:






∑

k
=
1


n



k

2


=


1
6


n
+


1
2



n

2


+


1
3



n

3




{\displaystyle \sum _{k=1}^{n}k^{2}={1 \over 6}n+{1 \over 2}n^{2}+{1 \over 3}n^{3}}





 Task
 show the first 10 rows of Faulhaber's triangle.
 using the 18th row of Faulhaber's triangle, compute the sum: 




∑

k
=
1


1000



k

17




{\displaystyle \sum _{k=1}^{1000}k^{17}}

 (extra credit).
 See also

 Bernoulli numbers
 Evaluate binomial coefficients
 Faulhaber's formula (Wikipedia)
 Faulhaber's triangle (PDF)


function bernoulli(n)
    A = Vector{Rational{BigInt}}(undef, n + 1)
    for i in 0:n
        A[i + 1] = 1 // (i + 1)
        for j = i:-1:1
            A[j] = j * (A[j] - A[j + 1])
        end
    end
    return n == 1 ? -A[1] : A[1]
end
 
function faulhabercoeffs(p)
    coeffs = Vector{Rational{BigInt}}(undef, p + 1)
    q = Rational{BigInt}(1, p + 1)
    sign = -1
    for j in 0:p
        sign *= -1
        coeffs[p - j + 1] = bernoulli(j) * (q * sign) * Rational{BigInt}(binomial(p + 1, j), 1)
    end
    coeffs
end
 
faulhabersum(n, k) = begin coe = faulhabercoeffs(k); mapreduce(i -> BigInt(n)^i * coe[i], +, 1:k+1) end
 
prettyfrac(x) = (x.num == 0 ? "0" : x.den == 1 ? string(x.num) : replace(string(x), "//" => "/"))
 
function testfaulhaber()
    for i in 0:9
        for c in faulhabercoeffs(i)
            print(prettyfrac(c), "\t")
        end
        println()
    end
    println("\n", prettyfrac(faulhabersum(1000, 17)))
end
 
testfaulhaber()
 