Sums of Powers of the Natural Numbers

Click the arrow to go up a level to the Math Index.

Follow the above link or click the graphic below to visit the Homepage.

Sums of Powers
of the Natural Numbers
Jim Cullen

The Harmonic Sum as an Introduction to Advanced Functions

The Harmonic Sum is a deceivingly simple concept that quickly leads to devilish infinities and then, upon closer inspection, we are introduced to some elegant solutions when the sum is expressed in terms of some advanced functions. Few mathematical constructions are as rich in the deeper meanings of number theory. Some of the most profound, such as the Riemann Hypothesis, are tied to theorems involving the prime numbers and remain unsolved to this day.

So many important interconnected mathematical ideas are drawn into the analysis of such a simple concept as the Harmonic Sum that I'm surprised it isn't required study at the high school level or freshman college level. Even the calculators commonly used by these students lack ANY of the advanced functions necessary to evaluate the Harmonic Sum or related mathematical problems. The one exception would be the Hewlett-Packard HP48/49/50 series scientific graphing calculators, the latest incarnation of which is the remarkable HP50g. Oddly enough it is not very likely that a high school or college student in the United States could be found using an HP calculator in their math studies. Whether or not there is an inference to be made here is a matter of opinion and, ultimately, depends on the skill of the user. To Texas Instruments' credit, TI89/92/V200 devices may now be armed to the teeth with powerful calculator software suites such as Bhuvanesh Bhatt's excellent Mathtools, closing the gap with HP calculators in the areas of advanced special functions and number theory applications. Also, very capable and efficient functions may also be implemented by the individual user on either brand of calculator.

Preference of calculating devices aside, we begin with a cursory examination of the Harmonic Sum. We will cover the necessary areas in some detail, step by step for some of the more important concepts but we will quickly gain momentum after that. Where required, more detail will be given but I will avoid excess detail as much as is possible in order to preserve an unobstructed view of how the related concepts are interwoven into one mathematical tapestry. Standard texts should be consulted or online resources visited if more detailed information is desired in any one area. Links at the bottom of this page to articles located at Wolfram's Mathworld are provided as a starting point for further research.

The mathematical concept of the Harmonic Sum is not a difficult one. It is defined as the sum of the inverses of the natural numbers from one to infinity. This sum can be expressed as H = 1/1 + 1/2 + 1/3 + 1/4 + 1/5 + ^{. . .} or " one, plus one half, plus one third, plus one fourth, plus one fifth, etc ". Each term is smaller than the last and, as we approach infinity, the increments to the sum become infinitely small. Nothing could be simpler . . . and nothing could be further from the truth! The surprise here is that the sum does not converge - there is no finite number as a solution to the problem - the sum continues to increase towards infinity and the ultimate sum, as we approach infinity, is itself infinite in value. This can be expressed as a summation in the following way:

	inf
H =		1 i	= infinity
	i = 1

Not only is the infinite sum H infinite in value but the finite sums H_n, which are the n'th harmonic numbers have no algebraic closed-form expressions as solutions. The finite sum is calculated by summing the inverses of the integers from one to some fixed number n instead of to infinity. Though the final sum is now some finite value, there is no closed-form algebraic solution to express the sum. In other words:

	n
H_n =		1 i
	i = 1

has no simple algebraic solution. Just about anytime one encounters reciprocal expressions within a summation it can be expected that there will be difficulties in developing closed-form algebraic solutions for them. This doesn't mean that there can't be some other solution to be found in a simple expression involving some higher functions. In this case there is just such a solution to be found in the Psi Function ( pronounced 'sigh' ). For our purposes, this function will be first introduced as a solution to the Harmonic Sum rather than in the terms of its own definition. As we will use it and in its simplest form, the Psi Function requires an integer argument and returns a real number. It can be shown that:

	n
H_n =		1 i	= Psi( n + 1 ) - Psi( 1 )
	i = 1

where ( - ) Psi( 1 ) is numerically equal to and is the basic definition of a well known mathematical constant called the Euler-Mascheroni Constant. The value of this constant is approximately equal to 0.577215664902. In the above summation however, Psi( 1 ) is present only because our summation begins at one. If we were to redefine our summation as the sum of the inverses of the natural numbers from m to n then the solution would be:

	n
H_n - H_m-1 =		1 i	= Psi( n + 1 ) - Psi( m )
	i = m

which is a surprisingly elegant solution considering the alternative if the Psi Function was not utilized. So what is the Psi Function? It is the derivative of the natural logarithm of the Gamma Function. Basically, the Gamma Function is very similar to the factorial except the Gamma Function allows real numbers as arguments instead of integers only. For positive integers, Gamma( n ) = ( n - 1 )!. It is here, with the introduction of the Gamma Function within the definition of the Psi Function, that we first begin to overlap the domain of number theory and take our first steps towards the deeper meanings surrounding the Prime Number Theorem and the Riemann Hypothesis. Some may think that to be a silly statement but consider that the Harmonic Sum is just a special case of sums of inverses of powers of the natural numbers - which in turn evokes thoughts of the Riemann Zeta Function, the solution of which can have a form very similar to the solution given for the Harmonic Sum. If anything is intimately related to unsolved problems in analytic number theory then it would have to be the Riemann Zeta Function which has its roots with the Harmonic Sum...

Reciprocal Sums of Higher Powers - Stepping Off Into The Deep End

We'll go through that again, slower this time for the rest of us! Consider the same type of summation we considered in the above Harmonic Sum except that we will square each term in the sum - like this:

H_n,2

( 1 )²

( 2 )²

( 3 )²

( 4 )²

( 5 )²

+ ^{. . .}

In this form the Harmonic Sum would be written as:

H_n,1

( 1 )¹

( 2 )¹

( 3 )¹

( 4 )¹

( 5 )¹

+ ^{. . .}

One notices immediately that the Harmonic Sum H_n,1 has the same form as H_n,2. The only difference is the power in the denominator. One may as well take the sums and rename them H_n,p, the subscript denoting the value of the power in the denominator. At first glance one would guess that H_n,2 would also result in a sum with an infinite value since we have squared each term in the summation. In fact the opposite occurs - we have squared the inverse of each term such that each term grows smaller much faster than in the Harmonic Sum. The end result is that the sum of the inverses of the squares of the natural numbers, from one to infinity, is a finite value - the summation converges. We will use a shorthand notation S_p to denote a summation of powers of natural numbers to explore some of the convergent values in the sums. Here is the first result for the sums of the inverses of the natural numbers squared, which can be written as:

	inf
S₂ =		1 i²	=	pi^2 6
	i = 1

where pi is the well-known ratio of the circumference to the diameter of a circle, roughly equal to 3.14159265359. The value of pi squared over six is approximately 1.64493406685. This is a curious result; that squaring the natuaral numbers and adding their inverses should result in a sum with the value of pi as a factor. A little curiosity and experimentation will get you everywhere so we try other integer values for the powers in the denominator. One soon discovers that the even powers result in summations whose limit is a simple expression involving pi. Here are the next three results, for powers in the denominator equal to four, six, and eight:

	inf
S₄ =		1 i⁴	=	pi^4 90
	i = 1

	inf
S₆ =		1 i⁶	=	pi^6 945
	i = 1

	inf
S₈ =		1 i⁸	=	pi^8 9450
	i = 1

This pattern will continue for higher even powers; the limit will be equal to that power of pi multiplied by some rational number or fraction. At this point we begin to drag in new mathematical concepts such as the Bernoulli Numbers, which are partially responsible for the rational numbers that the powers of pi are multiplied by in the above summations. Take note at this point that the odd powers are not nearly as elegant as the even powers, at least not in terms of the powers of pi and the Bernoulli Numbers. We might however attempt to revisit the Psi Function to see if perhaps the summations may be better behaved in terms of this function. The answer is yes but we should first discuss a more complete definition of the Psi Function.

It has already been mentioned that the Psi Function is the derivative of the natural logarithm of the Gamma Function. This is the base definition, often notated as Psi₀( x ) or the zero'th derivative of Psi( x ). We may continue taking derivatives of the Psi Function and increase our subscript appropriately, such that Psi_g( x ) denotes the g'th derivative of Psi( x ). This is the generalization of the Psi Function that is commonly called the Polygamma Function. Using this new definition, we will compose the solutions of the first eight sums of the inverses of the powers of the natural numbers; the powers from one through eight. In addition, the summations will be evaluated from a lower limit ( m ) to an upper limit ( n ). The odd powers are on the left side of the table and the even powers are on the right. All solutions, the values to which the summations converge, are given in terms of Psi_g( x ).

	n
S_p =		1 i^p
	i = m

Power
( p )

Sum

Power
( p )

Sum

Psi₀( n + 1 ) - Psi₀( m )

0 !

Psi₁( m ) - Psi₁( n + 1 )

1 !

Psi₂( n + 1 ) - Psi₂( m )

2 !

Psi₃( m ) - Psi₃( n + 1 )

3 !

Psi₄( n + 1 ) - Psi₄( m )

4 !

Psi₅( m ) - Psi₅( n + 1 )

5 !

Psi₆( n + 1 ) - Psi₆( m )

6 !

Psi₇( m ) - Psi₇( n + 1 )

7 !

Note that in the above table, when Power = 1, we have the same result as the Harmonic Sum. They are one and the same since Psi₀( x ) is just Psi( x ) and, by definition, 0! = 1 in the denominator ( the fact that 0! = 1 is itself a hotly debated topic ). Note also the difference between the even powers and the odd powers: the numerators are negated for the even powers as compared to the Harmonic Sum which is an odd power of one. Now that the solutions are in table form we can play with the limits of the summations to see how it changes the solutions. We may learn a few things. First we'll push the upper limit of the summations to infinity as follows:

	inf
S_p =		1 i^p
	i = m

Power
( p )

Sum

Power
( p )

Sum

- Psi₀( m )

0 !

Psi₁( m )

1 !

- Psi₂( m )

2 !

Psi₃( m )

3 !

- Psi₄( m )

4 !

Psi₅( m )

5 !

- Psi₆( m )

6 !

Psi₇( m )

7 !

What has happened is that, as the upper limit ( n ) approaches infinity, the value of Psi_g( n ) becomes vanishingly small, and we are left with the term for the lower limit. This form of summation, when the lower limit is variable, is actually the definition of what is known as the Hurwitz Zeta Function, what we could call a generalization of the Riemann Zeta Function where we can begin our summation at some other point besides one. From the table above we can also deduce a well known relationship between the Hurwitz Zeta Function and the Polygamma Function. We could also let the lower limit ( m ) of the summation defining the Riemann Zeta Function assume the value of one but the results will be uninteresting unless we set the upper limit ( n ) to infinity as in the above table. In this case we also have a set of vanished terms and the results are as follows:

	inf
S_p =		1 i^p
	i = 1

Power
( p )

Sum

Power
( p )

Sum

- Psi₀( 1 )

0 !

Psi₁( 1 )

1 !

- Psi₂( 1 )

2 !

Psi₃( 1 )

3 !

- Psi₄( 1 )

4 !

Psi₅( 1 )

5 !

- Psi₆( 1 )

6 !

Psi₇( 1 )

7 !

Using the results in the last few tables, it would be easy enough for just about anyone to immediately write down the resulting value for the summation of inverse natural numbers to any integer power. This is a function in itself: one supplies the argument, in this case the power ( p ), and the expression for the value to which the sum converges is returned by the function. This is in fact the purest form of the Riemann Zeta Function. It is also clear from the above table that the appearance of powers of pi in the sums of the even powers is due to the Psi Function. In fact there is a well known expression for the even powers that includes this power of pi in the formula:

For even integer p > 0
Z_R( p ) =	2^p-1 abs( B_p ) pi^p p !

where abs( B_p ) is the absolute value of the p'th Bernoulli Number. There are no such elegant formulas known for odd values of p. There are some very complicated expressions that do contain powers of pi for odd p but it would not be practical to discuss them here since, once the power of pi is extracted, the remainder of the expression is nothing like the rational numbers that accompany the powers of pi for even values of p.

Congratulations, you now have a seat-of-the-pants understanding of one of the most important special functions of all time in physics and mathematics. The Riemann Zeta Function is deeply rooted in fundamental mathematics and it turns up in studies of definite integrals, number theory, and the Prime Number Theorem. The notation I will use for the Riemann Zeta Function on this page is Z_R( p ). Now we can state the basic definition of the Riemann Zeta Function:

	inf
Z_R( p ) =		1 i^p	=	(-1)^p ^. Psi_p-1( 1 ) ( p - 1 ) !
	i = 1

and the Hurwitz Zeta Function:

	inf
Z_H( p , a ) =		1 ( i + a )^p	=	(-1)^p ^. Psi_p-1( a ) ( p - 1 ) !
	i = 0

where I have used ( -1 )^p to control the sign of the answer depending on whether ( p ) is an odd or even power. I could have instead just included an absolute value around the Psi Function to ensure that the result would always be positive. The Hurwitz Zeta Function is one generalization of the Riemann Zeta Function where the summation begins at some number, other than one, and then the summation carries on to infinity from there, in the same spirit that the Incomplete Gamma Function serves the Gamma Function. Using this definition, one can say that Z_H( p , 0 ) is the same as Z_R( p ). The Riemann Zeta Function may also be expressed as an infinite product involving the prime numbers. This is an intriguing expression but not altogether unexpected; due to the fundamental theorem of arithmetic there are many other expressions that contain the natural numbers in their terms that can be reformulated as expressions involving the prime numbers. In this case, the infinite product definition of the Riemann Zeta Function is known as the Euler Product. It was uncovered in the 18'th century by Leonhard Euler and is:

	inf
Z_R( s ) =		1 1 - ( P_k )^{- s}
	k = 1

where P_k is the k'th Prime Number. Though the definition based on the prime numbers is correct, it is not a very practical method for calculating the Riemann Zeta Function, especially for the lower arguments. As an example, while calculating Z_R( 2 ), ten primes will give us just two digits of accuracy. One hundred primes gives us just four digits of accuracy and one thousand primes still only gives us a five digits of accuracy - and that's if we are willing to wait for the calculation results. There are ways to accelerate these kinds of calculations but the infinite product is still not a very practical way to evaluate the Riemann Zeta Function. These infinite products involving the prime numbers do serve a greater purpose in number theory and also stand as examples of the elegance and consistency often encountered in advanced mathematics.

Notice that in the infinite sum formula given for the even powers which involve the Bernoulli Numbers, it seems that the alternating sign nature of the Psi Function should be attributed to the presence of the Bernoulli Numbers. However, earlier on we attributed this alternating sign to the Gamma Function within the definition of the Psi Function. So which is correct, or is this a sign that the Gamma Function is somehow related to the Bernoulli Numbers? Now we have an infinite product formula that has no alternating signs at all! It is after all a matter of definition; whatever you use to express another quantity must take into account exactly how it is related to that quantity and these interrelations are not always on equal terms.

Though we know it as the Riemann Zeta Function, Bernhard Riemann didn't actually enter the picture until about mid-19'th century when he extended the use of the Zeta Function to the complex plane. A century earlier, a now famous Swiss mathematician by the name of Leonhard Euler forged the link between the Harmonic Sum, the Zeta Function, and the Prime Numbers. Euler, knowing the Harmonic Sum was divergent ( has an infinite sum ), wondered if the Harmonic Sum of Primes was also divegent. He managed to split the Harmonic Sum into two series, one prime and one not prime, using a clever workaround in order to avoid the rules of infinite sums. In the end Euler had another proof of an infinite number of primes, he found that the Harmonic Sum of Primes was divergent, and he was able to express the Zeta Function as an infinite product involving only the Prime Numbers.

Going Beyond The Infinite for the Other Side of the Story

Having progressed from Harmonic Sums to the Riemann Zeta Function, and having touched on the concept of the Prime Number Theorem as it relates to the Riemann Hypothesis, we find we are at a good point to pause and take a closer look at some of the unsettled issues with the nature of the summations we have considered thus far. During our observations, we've managed to draw in so many concepts and advanced functions that it's easy to lose sight of what the basic driving force is behind the form and behavior of the resulting expressions. As is often the case in mathematics - when you've arrived at some advanced stage of finding solutions without finding reasons - going back to inspect a generalization of the problem by reformulating with inverses or negatives of key values will often provide much needed insight into the nature of the solutions you are finding.

In our case, note that we began with the Harmonic Sum and worked our way through to the Riemann Zeta Function, looking at specific as well as general cases of these forms of expressions. We have neglected to fully investigate one key parameter . . . before you read on you may want to try scanning quickly through everything again and seeing if you can spot it. Give up yet? It's hard to blame anyone for missing it since the concept doesn't come naturally but - we also want to look at negative powers! In our Harmonic Sums, Psi Function, and the Zeta Function, we've inspected the 'sums of the powers of natural numbers' where the 'power' has always been greater than zero.

Appendix 1: Numeric Evaluation of the Bernoulli Numbers

The Bernoulli Numbers are rational numbers that are very important in number theory and analytics. They also appear in series expansions of trigonometric functions and special forms of summations. The Bernoulli Numbers are difficult to evaluate quickly - for the same basic reasons that integer factoring is difficult to do quickly. The Bernoulli Numbers seem almost random in nature even though there are precise regular patterns that can be constructed to obtain their exact values. The most elegant of these constructions is based on the divisors of the natural numbers which in turn is determined by the distribution of the Prime Numbers - explaining the importance of the Bernoulli Numbers in number theory. We will use B_n to notate the n'th Bernoulli Number. which can be evaluated by using the following summation involving the Prime Numbers:

		inf
B_{2 n} =	2 ^. ( -1 )^n-1^. ( 2 n ) ! ( 2 pi )^{2 n}		P^-2n
		P = 1

where P are the Prime Numbers. The process is tediously slow however and there are better methods for the task. These usually involve some advanced functions, recursive methods, or double sums that are nearly as time consuming to evaluate so some shortcuts need to be considered to quickly return precise estimates of the Bernoulli Numbers. Modern methods for quickly evaluating exact values by calculating the numerator and denominator of B_n separately have been developed in recent years. However, I will again assume evaluation of the Bernoulli Numbers is to be performed on a handheld calculating device and so I will again present routines for exact or at least 11 or 12 digits of accuracy to be utilized as a floating-point result. There are other routines in the literature, again tediously slow for high values of B_n, that will calculate exact rationals for all B_n.

Review again what we discussed regarding the definition of the Riemann Zeta Function - recall that the Bernoulli Numbers are part of one formula for calculating Z_R( p ) for even values of p. Here n will be our value of p. This works out just fine since we don't need an expression for odd values of p since the odd Bernoulli Numbers are all equal to zero by definition. When n is very large then Z_R( n ) is very very close to equalling one. So, for a large enough value of n, we can assume that Z_R( n ) = 1 and rearrange the formula to solve for the value of B_n. Some years ago Simon Plouffe made use of just such a shortcut, along with several other innovations, in his record-breaking work with the calculation of the Bernoulli Numbers. In our calculation eleven or twelve significant digits will suffice and Z_R( n ) = 1 to this level of precision roughly when n > 38. Calculating the Bernoulli Numbers, in floating-point format, is then a very quick calculation that is simply:

for even integer n > 38
B_n =	2 ^. n ! ^. ( - 1 )^{( n/2-1 )} ( 2 pi )ⁿ

For B_n where n < 40, I've opted to store the rational values of the even Bernoulli Numbers in a list to be used as a look-up table as needed. Odd Bernoulli Numbers are again equal to zero by definition. This was done as a tradeoff between memory and speed. In this way, nearly instantaneous results at a precision of 12 significant digits for the Bernoulli Numbers can be obtained. When n < 40 the exact rational expression from the list is returned. When n > 38 the value returned is a real number in floating-point format.

Appendix 2: Numeric Evaluation of the Psi Function Psi₀( x )

One of the first things to take note of is that the Psi Function approaches the Natural Logarithm as the argument approaches infinity. The convergence is very gradual however so using ln ( x ) is no substitute for Psi ( x ), even at relatively high values substituted for x. As an example, ln( 1000 ) = 6.90775527898 while Psi₀( 1000 ) = 6.90725519565 so there is barely four digits of agreement between the two, even at x = 1000.

There are several ways to get a closer estimate of Psi₀( x ). One that I dreamt up for myself and had used for many years is one simple formula that gives six or seven digits of accuracy at x = 10, eight digits of accuracy at x = 25, about nine digits of accuracy at x = 50, and nearly the full 12 digits of accuracy by the time one reaches x = 100. This simple formula, taken from my own notes written during the Summer of 2000, is the following:

Psi₀( x ) =

ln ( x² - x )

6 x² - 6 x +1

The formula is very simple, easy to remember, and calculates extremely fast for either real or complex arguments. For more digits of accuracy it is generally no problem to then apply the recurrence relation to lift a low argument to a higher one where accuracy is better, then backstep the results to the lower argument. This recurrence formula defines how the value of Psi₀( x + 1 ) is related to the value of Psi₀( x ) and is as follows:

Psi₀( x + 1 ) = Psi₀( x ) + x^-1

which is simply a matter of subtracting a series of reciprocals. To apply the recurrence relation you subtract the reciprocal of the elevated argument ( minus one ) and the result is then the value of Psi for the elevated argument ( minus one ). Example: Psi( 17.3 ) - 1/16.3 = Psi( 16.3 ). It is even possible to evaluate negative values of x if the reflection formula for Psi₀( x ) is used. This property states:

Psi₀( - x ) = Psi₀( x ) + pi / tan [ pi ( x + 1 ) ]

where pi = 3.14159265358979 approximately. For more precise evaluation of the Psi Function, professional math software packages and calculating devices usually employ the Asymptotic Expansion expression for Psi₀( x ). For most purposes, only the first few terms of the expansion are required. The expression again involves the Bernoulli Numbers but, since only a few terms are usually needed, these can be simply squirreled away as constants within the calculation itself. The expression is:

Psi₀

ln ( x )

2 x

12 x²

120 x⁴

252 x⁶

240 x⁸

132 x¹⁰

+ ^{. . .}

which will be rewritten in code as nested multiplications of x² and constants to avoid having to calculate time-consuming higher powers of x. Generally one will find that only the terms up to the one involving x⁸ are used and the Asymptotic Expansion is also good for real or complex numbers. The expression is then very accurate as long as the argument x is greater than about 16. If x = 16 then the term ( 240 x⁸ )^-1 = 9.7 E -13, representing just about twelve digits of accuracy. For arguments below this limiting value of 16, the recurrence property is again utilized to maintain precision in the numerical results. For my own implementation, which was written specifically for small arguments, I've included terms up to x¹⁴ so that only 9 iterations of the recurrence relation are required instead of 16. If the derivative of the above expansion is taken multiple times it is also possible to calculate higher derivatives of the Psi Function and so is also a handy formula for evaluating the Polygamma Function though actual practice of this method is easier said than done, especially for higher derivatives. The Asymptotic Expansion is based on the following summation which converges fairly quickly to the required value of Psi₀( x ):

			inf
Psi₀( x ) = ln ( x ) -	1 2 x	-		B_{2 k} 2 k^.x^{2 k}
			k = 1

Appendix 3: Numeric Evaluation of the Polygamma Function Psi_n( x )

End Notes And Links To More Information On Relevant Mathematical Topics

Besides the standard texts, perhaps the greatest online mathematical resource is Wolfram's Mathworld, from the makers of Mathematica software. Hardware support for much of the work on this page includes: Texas Instruments' Voyage-200 PLT ( 12MHz M68000 CPU ) for summations, algorithmic work, numerical analysis, and some algebraic operations with the help of Bhuvanesh Bhatt's Mathtools calculator software; Hewlett-Packard's HP50G ( 75MHz ARM CPU ) for summations and some algebraic operations involving advanced functions; Microsoft Quick-Basic Ver4.5 for generation of tables of high-speed numeric sums and cross-checking. The best stuff, as always, is worked out on pencil and paper ( CPU supplied by user ).

For more information, please refer to the following articles at Wolfram's Mathworld.

Use your Back Button or click here to go to the Math Index