Click here to Skip to main content
15,868,164 members
Articles / Programming Languages / Python
Tip/Trick

Faulhaber made easy

Rate me:
Please Sign up or sign in to vote.
4.47/5 (5 votes)
2 Jul 2014CPOL2 min read 18.6K   6   8
Computation of the Faulhaber polynomials coefficients

Introduction

Faulhaber’s formula expresses the sum of the Image 1th power of the Image 2 first integers as a function of Image 3 and Image 4.

Image 5
Image 6
Image 7
Image 8

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

During the last 30 years, I have had to establish it by hand for moderate Image 11 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 Image 12. It must be so, because the first order finite difference Image 13 is a polynomial of degree Image 14 (namely Image 15). And it has no independent term, given that Image 16. Also observe that the coefficient of the leading term is Image 17 (just like the primitive of Image 18), and the sum of the coefficients is Image 19 (by Image 20).

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 Image 21, not too simple, nor too complex. As said, the sum of the cubes gives rise to a fourth degree polynomial of four terms:

Image 22

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

Image 27,

or

Image 28

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

Image 30

 

 

 

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

Divide every equation by the coefficient of its rightmost term:

Image 31

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

Image 32

Using the code

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

C++
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):

C++
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.

C++
# 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:

C++
# 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:

C++
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.

License

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


Written By
CEO VISION fOr VISION
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.

Comments and Discussions

 
QuestionVote of 4 Pin
Kenneth Haugland16-Apr-15 6:15
mvaKenneth Haugland16-Apr-15 6:15 
AnswerRe: Vote of 4 Pin
YvesDaoust16-Apr-15 7:17
YvesDaoust16-Apr-15 7:17 
GeneralRe: Vote of 4 Pin
Kenneth Haugland16-Apr-15 9:23
mvaKenneth 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
YvesDaoust14-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
YvesDaoust14-Jul-14 12:06 
GeneralRe: My vote of 3 Pin
CatchExAs14-Jul-14 12:11
professionalCatchExAs14-Jul-14 12:11 
Thanks for the explanation. Will change my vote Smile | :)

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.