Click here to Skip to main content
Click here to Skip to main content

Tagged as

Faulhaber made easy

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

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.

License

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

Share

About the Author

YvesDaoust
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

 
GeneralMy vote of 3 PinprofessionalCatchExAs14-Jul-14 12:40 
GeneralRe: My vote of 3 [modified] PinmemberYvesDaoust14-Jul-14 12:47 
GeneralRe: My vote of 3 PinprofessionalCatchExAs14-Jul-14 12:55 
GeneralRe: My vote of 3 PinmemberYvesDaoust14-Jul-14 13:06 
GeneralRe: My vote of 3 PinprofessionalCatchExAs14-Jul-14 13:11 

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

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

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