«  
Home
  »

Sum of Integer Powers (Part 3)

By , on January 19th, 2018 (12,277 views) | 2 comments - Leave a comment!
Topic: Mathematics, Maths--Technical
Spread the love

(Discrete Mathematics Techniques III)

1st ed. Apr 2nd, 2010

Abstract
This is the last in the 3-part series of articles on finding for oneself the solution to the sum of integer power problem, and in the process discovering the Bernoulli numbers. In Part 3 (this paper), we find a direct closed-form solution, i.e. one that does not require iteration, for the general case of the finite-summation-of-integer-powers problem S_p(N) = \sum_{k=1}^{N} k^p. Having established in Part 2 that the closed-form solution is a polynomial, the summation is here rewritten as the sum of the p+1 independent monomials a_j N^j (1 \leq j \leq p+1), where the a_j are unknown coefficients. Using the recurrence relation S_p(N+1) = S_p(N) + (N+1)^p, we obtain a linear combination of the monomials, which reduces to an easily solvable (p+1)-by-(p+1) triangular linear system in the unknown coefficients a_j of the closed-form polynomial solution. Maxima and Octave/Matlab codes for directly computing the closed-form solutions are included in the Appendices.

A lovely paper by Bearden (March 1996, American Mathematical Monthly), which was shared with me by a reader, tells the mathematical story nicely, with much of the history filled in.

NOTE: This paper contains many formulas. As such, you may prefer reading the paper as a typeset PDF.

Finding a Closed-Form solution for S_p(N) = \sum_{k=1}^N k^p using a Direct Method

Motivation

Our goal is to obtain, without using iteration, a closed form solution for the general case of the finite-summation-of-integer-powers problem:

\displaystyle S_p(N) = \sum_{k=1}^N k^p, \mbox{ (where } p \in \mathbb{N} = \{0,1,2,\ldots\}; N\geq0) \ \ \ \ \ (1)

Recall that in Part 2 of this paper , we used a method of recurrence relations to obtain the following results about S_p(N):

  1. S_p(N) may be expressed as a linear recurrence relation that uses the p-1 lower order formulas S_j(N) and S_j(N-1), \ \ (j = 1,\ldots,p-1) as follows:

    \displaystyle S_p(N) = \frac{1}{1+\alpha_p}\left\{N^2 + N^p + \sum_{j=1}^{p-1} C_j(N) S_j(N-1) - \sum_{j=1}^{p-1} \alpha_j S_j(N) \right\}, \ \ \ \ \ (2).


    The \alpha_j in (2) are the coefficients from the (p-1)-st polynomial solution S_{p-1}(N), while the terms C_j(N) are defined in terms of binomial coefficients1 as:

    \displaystyle C_j(N) = \left[\binom{p-1}{j}N - \binom{p-1}{j-1}\right].

  2. The closed form solution of S_p(N) is a polynomial in N of degree p+1 with rational coefficients and a constant term equal to zero.

While the expression for S_p(N) in (2) is indeed accurate, it requires repeated iteration to obtain the closed form expression for any particular p. Such an approach without the assistance of a computer algebra system such as Maxima2 would be prohibitively time-consuming and prone to error.

We are hence left with the following question:

Can we find a closed-form solution for the coefficients of the general solution polynomial that can be obtained directly, i.e. that does not require iteration?

This paper uses linear algebra and matrices to achieve precisely this goal.

The Direct Approach

Polynomial with Undetermined Coefficients

Since we have established in (2) that S_p(N) is a polynomial of order p+1 with no constant term, we may write down the closed form solution of S_p(N) in the form:

\displaystyle S_p(N) := \sum_{j=1}^{p+1} a_j N^j\ \ \ \ \ (3)

where a_j, \ (j = 1,\ 2,\ \ldots,\ p+1), are coefficients that have yet to be determined. As a first step towards determining these coefficients, observe that every finite summation can be written as a first order recurrence in N by peeling off the last summand:

\displaystyle S_p(N+1) = S_p(N) + (N+1)^p.\ \ \ \ \ (4)

Substituting (3) into (4) gives:

\displaystyle \sum_{j=1}^{p+1}a_j (N+1)^j = \sum_{j=1}^{p+1} a_j N^j + (N+1)^p.\ \ \ \ \ (5)

Summation Manipulations

What follows is a sequence of summation manipulations3 and simplifications of (5) aimed at obtaining an expression from which a closed-form solution for the coefficicents a_{j} becomes transparent:

  • Expand the binomial powers using the binomial formula on both the left-hand side (LHS) and right-hand side (RHS) of (5):

    \displaystyle \sum_{j=1}^{p+1}a_j \sum_{i=0}^j \binom{j}{i}N^i = \sum_{j=1}^{p+1} a_j N^j + \sum_{j=0}^p \binom{p}{j}N^j. \ \ \ \ \ (6)

  • Modify the LHS and RHS of (6) as follows, exploiting the fact that \binom{n}{k} = 0 \mbox{ when } k>n:
    • LHS: Extend upper range of inner sum to match upper range of outer sum.
    • RHS: Peel off 0-index term and add a vacuous p+1-index term to last sum.
    • Result:

      \displaystyle \sum_{j=1}^{p+1}a_j \sum_{i=0}^{p+1} \binom{j}{i}N^i = \sum_{j=1}^{p+1} a_j N^j + 1 + \sum_{j=1}^{p+1} \binom{p}{j}N^j

    • LHS: Interchange order of summations.
    • RHS: Combine terms.
    • Result:

      \displaystyle \sum_{i=0}^{p+1} N^i \sum_{j=1}^{p+1}a_j \binom{j}{i} = 1 + \sum_{j=1}^{p+1} \left[a_j + \binom{p}{j}N^j\right]

    • LHS: Peel off 0-index term from the (now outside) summation on i.
    • RHS: Change dummy index variable.
    • Result:

      \displaystyle \sum_{j=1}^{p+1}a_j + \sum_{i=1}^{p+1} N^i \sum_{j=1}^{p+1}a_j \binom{j}{i} = 1 + \sum_{i=1}^{p+1} \left[a_i + \binom{p}{i}N^i\right]

  • Bring everything over to the left-hand side to obtain a homogeneous linear combination of the monomials \{N^i\}_{i=0}^{p+1}:

    \displaystyle -1 + \sum_{j=1}^{p+1}a_j + \sum_{i=1}^{p+1} \left[\left(\sum_{j=1}^{p+1}a_j \binom{j}{i} \right) - a_i - \binom{p}{i}\right] N^i = 0 \ \ \ \ \ (7)

Linear Independence Reveals Equations for Unknown Coefficients

Since the monomials \{N^i\}_{i=0}^{p+1} are a basis for the vector space \mathbb{P}_{N+1} of polynomials of degree less than or equal to N+1, any linear combination of the \{N^i\}_{i=0}^{p+1} will equal zero if and only if each coefficient equals zero.

Therefore, setting the coefficients from (7) equal to zero gives us a system of p+2 simultaneous linear equations in the p+1 unknown coefficients a_j, one equation for each of the coefficients of N^i:

\displaystyle (i=0)\ \ \ \ \sum_{j=1}^{p+1} a_j =1\ \ \ \ \ (8)

\displaystyle (i=1,\ldots,p)\ \ \ \ \sum_{j=1, j\neq i}^{p+1}\binom{j}{i}a_j = \binom{p}{i} \ \ \ \ \ (9)

\displaystyle (i=p+1)\ \ \ \ \sum_{j=1, j\neq p+1}^{p+1}\binom{j}{i}a_j = \binom{p}{i}\ \ \ \ \ (*)
Observe that (*) reduces to 0=0, since j<p+1 for all j, and \binom{n}{k}=0 when n<k. Thus we are left with a square system of p+1 equations (i=0,\ldots,p) in p+1 unknowns a_j, (j = 1,\ldots,p+1).

By observing that \binom{j}{i}=1 for i=0, we can combine (8) and (9) into the single set of linear equations:

\displaystyle (i=0, \ldots, p)\ \ \ \ \sum_{j=1, j\neq i}^{p+1}\binom{j}{i}a_j = \binom{p}{i} \ \ \ \ \ (10)

Taking into account the inequality j\neq i, and noting once again that \binom{n}{k}=0 when n<k, (10) simplifies to a triangular linear system:

\displaystyle (i=0, \ldots, p)\ \ \ \ \sum_{j=i+1}^{p+1}\binom{j}{i}a_j = \binom{p}{i} \ \ \ \ \ (11)

and hence all p+1 equations in the system are linearly independent.

From Linear System to Matrix Equation

The triangular linear system (11) may be readily solved by combining the p+1 summations into a single matrix equation and using any one of a number of matrix solver packages. But this requires that the system (11) be re-indexed4 so that the start value of index i matches that of index j. Following the usual matrix convention, referred to as 1-indexing, we choose the start indices to be i=j=1.

The 1-indexed system is:

\displaystyle (i=1, \ldots, p+1) \ \ \ \ \sum_{j=i}^{p+1}\binom{j}{i-1}a_j = \binom{p}{i-1} \ \ \ \ \ (12)

Maxima code for solving the 0-indexed5 triangular linear system (13) is given in Appendix B. Octave/Matlab code for solving the 1-indexed triangular linear system (12) is given in Appendix C.

The matrix equation equivalent to (12) is: \mathbf{M} \cdot \mathbf{a} = \mathbf{b}, where

\displaystyle \mathbf{M} := \begin{bmatrix} \binom{1}{0}&\binom{2}{0}&\ldots&\binom{p+1}{0}\rule[-.1cm]{0cm}{.1cm}\\ &\binom{2}{1}&\ldots&\binom{p+1}{1}\\ &\ldots&\binom{j}{i-1}\ldots\\ & & &\binom{p+1}{p} \end{bmatrix},\ \ \ \mathbf{a} := \begin{bmatrix}a_1\\a_2\\ \ldots\\ a_{p+1}\end{bmatrix},\ \ \mbox{and} \ \ \mathbf{b} := \begin{bmatrix}\binom{p}{0}\\\binom{p}{1}\\\ldots\\\binom{p}{p}\end{bmatrix}. \ \ \ \ \ (14)

Since \mathbf{M} is an upper triangular matrix with all non-zero upper triangular entries, we know we can readily solve this for any given values of p and N using a matrix solver system such as Maxima/Mathematica, Octave/Matlab, or Maple, among others. Source code for solutions in Maxima and Octave/Matlab are given in Appendix B and Appendix C, respectively.

Computing the General, Closed-Form Solution

We know from (2) that the general, closed-form solution to (1) is a p+1-degree polynomial (3) in N with rational coefficients a_j given by solving the triangular linear system (12), or equivalently (14).

The highest few coefficients a_j can be readily computed by hand. For all p, we have:

  1. a_{p+1} = \frac{1}{p+1}
  2. a_p = \frac{1}{2}\binom{p}{0} = \frac{1}{2}
  3. a_{p-1} = \frac{1}{12}\binom{p}{1} = \frac{p}{12}
  4. a_{p-3} = -\frac{1}{120}\binom{p}{3}
  5. a_{p-5} = \frac{1}{252}\binom{p}{5}
  6. a_{p-7} = -\frac{1}{240}\binom{p}{7}
  7. a_{p-9} = \frac{1}{132}\binom{p}{9}.
  8. In particular, we claim (without proof) that

    \displaystyle a_{p-i}=c_{i}\binom{p}{i}, \ \ \ \ \ (15)

    where the {c_{i}} are rational coefficients independent of p, and c_{i}=0 for all even positive i. Hence we have:

    \displaystyle a_{p-2} = a_{p-4} = \ldots = 0.

Substituting the above coefficients into the general closed-form polynomial solution gives us:

\displaystyle S_p(N) = \frac{1}{p+1} N^{p+1} + \frac{1}{2} N^p + \frac{p}{12}N^{p-1} -\frac{1}{120}\binom{p}{3}N^{p-3} + \frac{1}{252}\binom{p}{5}N^{p-5}

\displaystyle \ \ \ \ \ \ \ \ \ \ - \frac{1}{240}\binom{p}{7}N^{p-7} + \frac{1}{132}\binom{p}{9}N^{p-9} - \ldots + p c_{p-1} N. \ \ \ \ \ (16)

A listing of the solution formulas for the first ten p values is given in Appendix A.

Conclusion

Reviewing the path we have taken in this 3-part paper: Part 1 motivated the problem and illustrated a recurrence relation method for obtaining the solution for small p. Part 2 generalized this method and found a p-th order recurrence relation for the solution. Part 2 also illustrated, by an induction argument, that the closed form solutions for all p are polynomials of degree p+1 with rational coefficients and no constant term.

This paper used the closed-form polynomial expression motivated in Part 2 to obtain a direct solution. By writing out the polynomial form with undetermined coefficients, the recurrence relation (4) was manipulated into a linear combination of polynomials (7). Using a linear independence argument, (7) was reduced to a triangular linear system (12) and associated matrix equation (14).The general closed-form solution is (16).

Solving the finite-summation-of-integer-powers problem \sum_{k=1}^N k^p, for arbitrary positive integers p and N has provided a natural setting to use a variety of techniques from discrete mathematics and linear algebra, and poses additional interesting questions (characterization of the denominators in the coefficients of (16) (the c_i in (15)) and divisibility properties of S_p for various p). In particular, we have obtained a direct method for solving the finite-summation-of-integer-powers problem that invokes a matrix solution and does not require iteration.

Appendix A: Listing of Solutions to (1)

A listing of the solution formulas for the first ten p values is as follows:

\displaystyle S_1(N) = \sum_{k=1}^{N}k =\frac{{N}^{2}}{2}+\frac{N}{2}

\displaystyle S_2(N) = \sum_{k=1}^{N}{k}^{2}=\frac{{N}^{3}}{3}+\frac{{N}^{2}}{2}+\frac{N}{6}

\displaystyle S_3(N) = \sum_{k=1}^{N}{k}^{3}=\frac{{N}^{4}}{4}+\frac{{N}^{3}}{2}+\frac{{N}^{2}}{4}

\displaystyle S_4(N) = \sum_{k=1}^{N}{k}^{4}=\frac{{N}^{5}}{5}+\frac{{N}^{4}}{2}+\frac{{N}^{3}}{3}-\frac{N}{30}

\displaystyle S_5(N) = \sum_{k=1}^{N}{k}^{5}=\frac{{N}^{6}}{6}+\frac{{N}^{5}}{2}+\frac{5\,{N}^{4}}{12}-\frac{{N}^{2}}{12}

\displaystyle S_6(N) = \sum_{k=1}^{N}{k}^{6}=\frac{{N}^{7}}{7}+\frac{{N}^{6}}{2}+\frac{{N}^{5}}{2}-\frac{{N}^{3}}{6}+\frac{N}{42}

\displaystyle S_7(N) = \sum_{k=1}^{N}{k}^{7}=\frac{{N}^{8}}{8}+\frac{{N}^{7}}{2}+\frac{7\,{N}^{6}}{12}-\frac{7\,{N}^{4}}{24}+\frac{{N}^{2}}{12}

\displaystyle S_8(N) = \sum_{k=1}^{N}{k}^{8}=\frac{{N}^{9}}{9}+\frac{{N}^{8}}{2}+\frac{2\,{N}^{7}}{3}-\frac{7\,{N}^{5}}{15}+\frac{2\,{N}^{3}}{9}-\frac{N}{30}

\displaystyle S_9(N) = \sum_{k=1}^{N}{k}^{9}=\frac{{N}^{10}}{10}+\frac{{N}^{9}}{2}+\frac{3\,{N}^{8}}{4}-\frac{7\,{N}^{6}}{10}+\frac{{N}^{4}}{2}-\frac{3\,{N}^{2}}{20}

\displaystyle S_{10}(N) = \sum_{k=1}^{N}{k}^{10}=\frac{{N}^{11}}{11}+\frac{{N}^{10}}{2}+\frac{5\,{N}^{9}}{6}-{N}^{7}+{N}^{5}-\frac{{N}^{3}}{2}+\frac{5\,N}{66}

Appendix B: Maxima Source Code for Solving (14)

For particular p, the coefficients a_j for the closed form solutions S_p(N) are immediate. Maxima code, is given below, for determining this solution using the 0-indexed triangular linear system (13):

     solution(p):= block([a, eq],	/* give subroutine variables local scope */
     v : makelist(a[i], i, 0, p),	/* create list of unknowns (0-indexed) */
					/* create list of equations (0-indexed) */
     eq : makelist(sum(binom(j+1,i)*a[j],j,i,p) = binom(p,i), i, 0, p),
     linsolve(eq, v)
  )$

A more elaborate function that computes the series form and factored forms of the solution is given below. This was used to generate the first ten solution listings.

	SpN_mat(p):= block([a, eq],	/* give subroutine variables local scope */
	v : makelist(a[i], i, 0, p),	/* create list of unknowns (0-indexed)   */
					/* create list of equations (0-indexed)  */
	eq : makelist(sum(binom(j+1,i)*a[j],j,i,p) = binom(p,i), i, 0, p),  
     
	/* find coefficients of solution polynomial by solving linear system */
	sol : linsolve(eq, v),
     
	/* create polynomial: inner product of {N^i} with coefficients */
	pol : makelist(N^(i+1),i,0,p),	/* monomials {N^i} */
     
	/* closed form formula */
	cff : rhs(sol.pol), 	/* inner product of {N^i} with coefficients */     
	cff 				/* return closed form formula in series */
  )$

Download the complete Maxima source code listing here.

Appendix C: Octave/Matlab Source Code for Solving (14)

For particular p, the coefficients a_j for the closed form solutions S_p(N) are immediate. Octave/Matlab code, is given below, for determining this solution using the 1-indexed triangular linear system (12):

% Functions listing coefficients from a_p+1 to a_1 as fractions

function c = sumkp_matrix(p)
    M=zeros(p+1); 		% initialize matrix M
    for i=1:p+1
        for j=i:p+1		% set upper triangular elements
            M(i,j)=nchoosek(j,i-1);  %  Equation (12)
        end;
    end;
    
    b=zeros(p+1,1); % initialize column matrix b
    for i=1:p+1
        b(i)=nchoosek(p,i-1);  %  Equation (12)
    end
    
    c=inv(M)*b;  % Solve to obtain coefficients
       
    c=fliplr(c');		% List from highest index to lowest
    sml=abs(c)<1e-8;	% Numerical correction: set tiny coeffients to exactly zero
    c(find(sml==1))=0;
    
    disp("Solution Coefficients (High to Low order):")
    c=rats(c);  		% return coefficients as fractions

end;

Download the complete Octave source code listing here.

This article is available for download as a PDF here.

The Maxima and Octave complete source code listings implementing the direct matrix solution are available for download here.

References

[Ebr10a]
Assad Ebrahim.
Finite summation of integer powers xp, part 1.
January 2010.

[EO10b]
Assad Ebrahim.
Finite summation of integer powers xp, part 2.
February 2010.

[GKP]
Ron Graham, Donald Knuth, and Oren Patashnik.
Concrete Mathematics: A Foundation for Computer Science.
Addison Wesley.

[Knu]
Donald Knuth.
The Art of Computer Programming (3 Volumes).

Footnotes

  1. The notation \binom{n}{k} used above denotes binomial coefficients, often expressed verbally as "n choose k". Other representations of these coefficients include C(n,k) and C_{k}^{n}. ↩
  2. Maxima code for automatically computing solutions to (1) using solution formula (2) is given in Part 2 of this paper . ↩
  3. Such manipulations are taught in material such as (GKP) and (Knu). ↩
  4. Index manipulation is covered in (GKP) and (Knu). ↩
  5. 0-indexing is more convenient for various computational software packages and programming langauges, including Maxima. The 0-indexed system is:

    \displaystyle (i=0,\ldots,p)\ \ \ \ \sum_{j=i}^{p}\binom{j+1}{i}a_j =\binom{p}{i} \ \ \ \ \ (13)

    ↩

If you enjoyed this article: :
 (1) Leave a comment below!
(2) Mention and reshare the article!

(3) Subscribe to our RSS feed! (Don't miss the next article!)

2 comments to Sum of Integer Powers (Part 3)

Leave a Reply

You can use these HTML tags

<a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

  

  

  

Your comments are valued! (Please indulge the gatekeeping question as spam-bots cannot (yet) do simple arithmetic...) - required

Optionally add an image (JPEG only)

 

«  
Home
  »

Stats: 1,066,417 article views since 2010 (March update)

Dear Readers:

Welcome to the conversation!  We publish long-form pieces as well as a curated collection of spotlighted articles covering a broader range of topics.   Notifications for new long-form articles are through the feeds (you can join below).  We love hearing from you.  Feel free to leave your thoughts in comments, or use the contact information to reach us!

Reading List…

Looking for the best long-form articles on this site? Below is a curated list by the main topics covered.

Mathematics-History & Philosophy

  1. What is Mathematics?
  2. Prehistoric Origins of Mathematics
  3. The Mathematics of Uruk & Susa (3500-3000 BCE)
  4. How Algebra Became Abstract: George Peacock & the Birth of Modern Algebra (England, 1830)
  5. The Rise of Mathematical Logic: from Laws of Thoughts to Foundations for Mathematics
  6. Mathematical Finance and The Rise of the Modern Financial Marketplace
  7. A Course in the Philosophy and Foundations of Mathematics
  8. The Development of Mathematics
  9. Catalysts in the Development of Mathematics
  10. Characteristics of Modern Mathematics

Electronic & Software Engineering

  1. Electronics in the Junior School - Gateway to Technology
  2. Coding for Pre-Schoolers - A Turtle Logo in Forth
  3. Experimenting with Microcontrollers - an Arduino development kit for under £12
  4. Making Sensors Talk for under £5, and Voice Controlled Hardware
  5. Computer Programming: A brief survey from the 1940s to the present
  6. Forth, Lisp, & Ruby: languages that make it easy to write your own domain specific language (DSL)
  7. Programming Microcontrollers: Low Power, Small Footprints & Fast Prototypes
  8. Building a 13-key pure analog electronic piano.
  9. TinyPhoto: Embedded Graphics and Low-Fat Computing
  10. Computing / Software Toolkits
  11. Assembly Language programming (Part 1 | Part 2 | Part 3)
  12. Bare Bones Programming: The C Language

Pure & Applied Mathematics

  1. Fuzzy Classifiers & Quantile Statistics Techniques in Continuous Data Monitoring
  2. LOGIC in a Nutshell: Theory & Applications (including a FORTH simulator and digital circuit design)
  3. Finite Summation of Integer Powers: (Part 1 | Part 2 | Part 3)
  4. The Mathematics of Duelling
  5. A Radar Tracking Approach to Data Mining
  6. Analysis of Visitor Statistics: Data Mining in-the-Small
  7. Why Zero Raised to the Zero Power IS One

Technology: Sensors & Intelligent Systems

  1. Knowledge Engineering & the Emerging Technologies of the Next Decade
  2. Sensors and Systems
  3. Unmanned Autonomous Systems & Networks of Sensors
  4. The Advance of Marine Micro-ROVs

Math Education

  1. Teaching Enriched Mathematics, Part 1
  2. Teaching Enriched Mathematics, Part 2: Levelling Student Success Factors
  3. A Course in the Philosophy and Foundations of Mathematics
  4. Logic, Proof, and Professional Communication: five reflections
  5. Good mathematical technique and the case for mathematical insight

Explore…

Timeline