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)


with javascript_semantics
include builtins\pfrac.e -- (0.8.0+)
 
function bernoulli(integer n)
    sequence a = {}
    for m=0 to n do
        a = append(a,{1,m+1})
        for j=m to 1 by -1 do
            a[j] = frac_mul({j,1},frac_sub(a[j+1],a[j]))
        end for
    end for
    if n!=1 then return a[1] end if
    return frac_uminus(a[1])
end function
 
function binomial(integer n, k)
    if n<0 or k<0 or n<k then ?9/0 end if
    if n=0 or k=0 then return 1 end if
    atom num = 1,
         denom = 1
    for i=k+1 to n do
        num *= i
    end for
    for i=2 to n-k do
        denom *= i
    end for
    return num / denom
end function
 
function faulhaber_triangle(integer p, bool asString=true)
    sequence coeffs = repeat(frac_zero,p+1)
    for j=0 to p do
        frac coeff = frac_mul({binomial(p+1,j),p+1},bernoulli(j))
        coeffs[p-j+1] = iff(asString?sprintf("%5s",{frac_sprint(coeff)}):coeff)
    end for
    return coeffs
end function
 
for i=0 to 9 do
    printf(1,"%s\n",{join(faulhaber_triangle(i),"  ")})
end for
puts(1,"\n")
 
if platform()!=JS then
    sequence row18 = faulhaber_triangle(17,false)
    frac res = frac_zero
    atom t1 = time()+1
    integer lim = 1000
    for k=1 to lim do
        bigatom nn = BA_ONE
        for i=1 to length(row18) do
            res = frac_add(res,frac_mul(row18[i],{nn,1}))
            nn = ba_mul(nn,lim)
        end for
        if time()>t1 then printf(1,"calculating, k=%d...\r",k) t1 = time()+1 end if
    end for
    printf(1,"%s        \n",{frac_sprint(res)})
end if
