Faulhaber made easy

Introduction

Faulhaber’s formula expresses the sum of the ^th power of the  first integers as a function of  and .

…

The  are known as the triangular numbers, and the  are the square pyramidal numbers, for obvious reason.

During the last 30 years, I have had to establish it by hand for moderate  on several occasions, and found the task tedious and error prone every time. Maybe some of you who tried experienced this too. I recently discovered an easy approach that I wanted to share.

Background

As can be seen, the expression is a polynomial of degree . It must be so, because the first order finite difference  is a polynomial of degree  (namely ). And it has no independent term, given that . Also observe that the coefficient of the leading term is  (just like the primitive of ), and the sum of the coefficients is  (by ).

The coefficients are closely related to the Bernouilli numbers. A detailed description is given in Wikipedia.

The sum of the cubes

For the sake of the explanation, we shall consider the case of , not too simple, nor too complex. As said, the sum of the cubes gives rise to a fourth degree polynomial of four terms:

The trick that makes it easy is to take the difference of  and , which we know to be . We will expand the binomials  using the binomial coefficients.

,

or

After simplification, we can identify the coefficients of the powers of , and we form this system:

(You should recognize Pascal’s triangle with alternating signs and a diagonal missing.)

Divide every equation by the coefficient of its rightmost term:

This triangular system is readily solved row by row, and as announced:

Using the code

I have coded in Python the computation of the Faulhaber polynomial for arbitrary . The coefficients being fractions, I handle them as such (pairs of integers), keeping the fractions simple by means of the gcd algorithm (Euclid).

def Simplify(a, b):
    # Divide a and b by their gcd
    c, d= a, b 
    while d != 0:
        c, d= d, c % d
    return (a / c, b / c)

To carry out all computation with fractions, it is enough to implement a SAXPY operation (multiply and add):

def SAXPY(A, X, Y):
    # Compute A.X + Y, where A, X and Y are fractions
    S= Simplify(A[0] * X[0], A[1] * X[1])
    S= Simplify(S[0] * Y[1] + S[1] * Y[0], S[1] * Y[1])
    return S

For ease of programming, I precomputed Pascal’s triangle up to the desired level, in a bidimensional array.

# Maximum order of the formula
K= 10

# Setup Pascal triangle

# Apex
Pascal= [[1]]

# Rows
for i in range(1, K + 2):
    Pascal.append([1])
    for j in range(1, i):
        # Recurrence relation
        Pascal[i].append(Pascal[i-1][j-1] + Pascal[i-1][j])
    Pascal[i].append(1)

The core of the computation, solving the triangula system, is straigthforward:

# Compute the Faulhaber polynomial
for k in range(K + 1):
    # Initialize the leading coefficient
    S= [(1, k+1)]

    # Compute the next coefficients from the triangular system
    for i in range(k):
        T= (0, 1)
        for j in range(i + 1):
            # Accumulate, with alternating signs
            T= SAXPY(S[j], Simplify(Pascal[k + 1 - j][k - 1 - i], (i - k if (i + j) & 1 else k - i)), T)
        S.append(T)

Here we are:

S0(n) = n
S1(n) = n^2/2 + n/2
S2(n) = n^3/3 + n^2/2 + n/6
S3(n) = n^4/4 + n^3/2 + n^2/4
S4(n) = n^5/5 + n^4/2 + n^3/3 - n/30
S5(n) = n^6/6 + n^5/2 + 5n^4/12 - n^2/12
S6(n) = n^7/7 + n^6/2 + n^5/2 - n^3/6 + n/42
S7(n) = n^8/8 + n^7/2 + 7n^6/12 - 7n^4/24 + n^2/12
S8(n) = n^9/9 + n^8/2 + 2n^7/3 - 7n^5/15 + 2n^3/9 - n/30
S9(n) = n^10/10 + n^9/2 + 3n^8/4 - 7n^6/10 + n^4/2 - 3n^2/20
S10(n) = n^11/11 + n^10/2 + 5n^9/6 - n^7 + n^5 - n^3/2 + 5n/66
S11(n) = n^12/12 + n^11/2 + 11n^10/12 - 11n^8/8 + 11n^6/6 - 11n^4/8 + 5n^2/12
S12(n) = n^13/13 + n^12/2 + n^11 - 11n^9/6 + 22n^7/7 - 33n^5/10 + 5n^3/3 - 691n/2730
S13(n) = n^14/14 + n^13/2 + 13n^12/12 - 143n^10/60 + 143n^8/28 - 143n^6/20 + 65n^4/12 - 691n^2/420
...

Points of Interest

It looks so easy now !

History

This is the first version.