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)


 
ft_rows(Lz) :-
    lazy_list(ft_row, [], Lz).
 
ft_row([], R1, R1) :- R1 = [1].
ft_row(R0, R2, R2) :-
    length(R0, P),
    Jmax is 1 + P, numlist(2, Jmax, Qs),
    maplist(term(P), Qs, R0, R1),
    sum_list(R1, S), Bk is 1 - S, % Bk is Bernoulli number
    R2 = [Bk | R1].
 
term(P, Q, R, S) :- S is R * (P rdiv Q).
 
show(N) :-
    ft_rows(Rs),
    length(Rows, N), prefix(Rows, Rs),
    forall(
        member(R, Rows),
            (format(string(S), "~w", [R]),
             re_replace(" rdiv "/g, "/", S, T),
             re_replace(","/g, ", ", T, U),
             write(U), nl)).
 
sum(N, K, S) :- % sum I=1,N (I ** K)
    ft_rows(Rows), drop(K, Rows, [Coefs|_]),
    reverse([0|Coefs], Poly),
    foldl(horner(N), Poly, 0, S).
 
horner(N, A, S0, S1) :-
    S1 is N*S0 + A.
 
drop(N, Lz1, Lz2) :-
    append(Pfx, Lz2, Lz1), length(Pfx, N), !.
 