Click here to Skip to main content
13,088,846 members (54,597 online)
Click here to Skip to main content
Add your own
alternative version

Tagged as


6 bookmarked
Posted 2 Jul 2014

Faulhaber made easy

, 2 Jul 2014
Rate this:
Please Sign up or sign in to vote.
Computation of the Faulhaber polynomials coefficients


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.


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.



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):
    for j in range(1, i):
        # Recurrence relation
        Pascal[i].append(Pascal[i-1][j-1] + Pascal[i-1][j])

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)

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 !


This is the first version.


This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)


About the Author

Belgium Belgium
I fell into applied algorithmics at the age of 16 or so. This eventually brought me to develop machine vision software as a professional. This is Dreamland for algorithm lovers.

You may also be interested in...


Comments and Discussions

QuestionVote of 4 Pin
Kenneth Haugland16-Apr-15 6:15
professionalKenneth Haugland16-Apr-15 6:15 
AnswerRe: Vote of 4 Pin
YvesDaoust16-Apr-15 7:17
memberYvesDaoust16-Apr-15 7:17 
GeneralRe: Vote of 4 Pin
Kenneth Haugland16-Apr-15 9:23
professionalKenneth Haugland16-Apr-15 9:23 
GeneralMy vote of 3 Pin
CatchExAs14-Jul-14 11:40
professionalCatchExAs14-Jul-14 11:40 
GeneralRe: My vote of 3 Pin
YvesDaoust14-Jul-14 11:47
memberYvesDaoust14-Jul-14 11:47 
GeneralRe: My vote of 3 Pin
CatchExAs14-Jul-14 11:55
professionalCatchExAs14-Jul-14 11:55 
GeneralRe: My vote of 3 Pin
YvesDaoust14-Jul-14 12:06
memberYvesDaoust14-Jul-14 12:06 
GeneralRe: My vote of 3 Pin
CatchExAs14-Jul-14 12:11
professionalCatchExAs14-Jul-14 12:11 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

Permalink | Advertise | Privacy | Terms of Use | Mobile
Web01 | 2.8.170813.1 | Last Updated 2 Jul 2014
Article Copyright 2014 by YvesDaoust
Everything else Copyright © CodeProject, 1999-2017
Layout: fixed | fluid