Newly Discovered Method for Summing nth Powers of Integers

[ktoz]

While responding to another thread about summing nth powers of integers, I came up with what might be a new method. There is an old one (Faulhabers formula) http://mathworld.wolfram.com/FaulhabersFormula.html" , but mine seems to be considerably simpler.

Most likely, someone has already discovered this, but I wouldn't know, no math training beyond high school...

$$\sum_{j=0}^m j^p \quad = \quad p! \sum_{j=0}^m Q_{pj} R_{pmj}$$

$$Q_{pk} \quad = \quad \left\{ \begin{array}{ll} \sum_{j=0}^k -1^j \frac{ (k - j + 1)^p}{j!(p - j)!} \quad & k < p \\ 1 & k \geq p \end{array}$$

$$R_{pmk} \quad = \quad \prod_{l=1}^p 1 + \frac{m - l}{k}$$

Ken

Added later: With a slight variation, equation also yeilds the following interesting relation

$$x^n \quad = \quad p! \sum_{j=0}^x Q_{nj} R_{n-1xj}$$

 

[shmoe]

What is $$R_{pmk}$$ when k=0? As it is you are dividing by zero. Ignoring this division by zero problem for the moment, you'd end up with a polynomial in m of degree p. This can't be correct. You'll notice that Faulhaber's gives a polynomial of degree p+1 for the sum of pth powers, and if yours is giving a polynomial then they must be equal.

In any case Faulber's still seems simpler as it tells you the coefficients of this polynomial directly (in terms of the Bernoulli coefficients), rather than a sum over m polynomials (as yours is summing over the R's)

 

[ktoz]

What is $$R_{pmk}$$ when k=0? As it is you are dividing by zero.

Right you are. I worked out all the details/got it working in AppleScript and got mixed up translating AppleScript to mathematical notation. It should be:

$$R_{pmk} \quad = \quad \prod_{l=1}^p 1 + \frac{m - k}{l}$$

Ignoring this division by zero problem for the moment, you'd end up with a polynomial in m of degree p. This can't be correct. You'll notice that Faulhaber's gives a polynomial of degree p+1 for the sum of pth powers, and if yours is giving a polynomial then they must be equal.

Does the above change solve the "p" vs "p+1" problem? The notation may be off, but the testing program works fine.

In any case Faulber's still seems simpler as it tells you the coefficients of this polynomial directly (in terms of the Bernoulli coefficients), rather than a sum over m polynomials (as yours is summing over the R's)

I checked out Bernoulli numbers http://mathworld.wolfram.com/BernoulliNumber.html" and, speaking in programming terms, $$\frac{x}{e^x - 1}$$ is an illusion of sorts in that you still have to perform a loop to calculate $$e^x$$. The real representation would be

$$\frac{x}{-1 + \prod_{j=0}^x e}$$

Since "e" is a trancendental, the calculations would be more intensive than the simple integers I'm working with in $$R_{pmk}$$. So Bernoullis may look cleaner symbolically, but in practice, wouldn't $$\prod_{j=1}^l 1 + \frac{m - j}{l}$$ actually run faster?

 

[DeadWolfe]

If that actually works, I would imagine the fact that it doesn't need the Bernoulli numbers would cause it run faster for large values of p.

 

[ktoz]

It works nicely (at least the script does). I'm a mathematical notation newbie, but I think it's correct now with the fix:

$$R_{pmk} \quad = \quad \prod_{l=1}^p 1 + \frac{m - k}{l}$$

 

[shmoe]

Latex doesn't seem to be turning up when I rty to preview, so here's hoping it works and I didn't make any typos.

Right you are. I worked out all the details/got it working in AppleScript and got mixed up translating AppleScript to mathematical notation. It should be:
$$R_{pmk} \quad = \quad \prod_{l=1}^p 1 + \frac{m - k}{l}$$

edit-the latex I had here still wasn't displaying properly, but it's redundant. I figured out what seemed to be wrong, what you've got gives the sum up to (m+1)^p, not just m^p. No biggie, but otherwise it seems to work.

Does the above change solve the "p" vs "p+1" problem? The notation may be off, but the testing program works fine.

Sorry, forget about that. Yours isn't going to be a polynomial at all (not in the way I was implying), I'm not sure what I was thinking.

I checked out Bernoulli numbers http://mathworld.wolfram.com/BernoulliNumber.html" and, speaking in programming terms, $$\frac{x}{e^x - 1}$$ is an illusion of sorts in that you still have to perform a loop to calculate $$e^x$$. The real representation would be
$$\frac{x}{-1 + \prod_{j=0}^x e}$$
Since "e" is a trancendental, the calculations would be more intensive than the simple integers I'm working with in $$R_{pmk}$$. So Bernoullis may look cleaner symbolically, but in practice, wouldn't $$\prod_{j=1}^l 1 + \frac{m - j}{l}$$ actually run faster?

e being transcendental doesn't matter (and I don't understand your "real representation"). I don't know the most efficient way to calculate bernoulli numbers, but if you need up to the pth the naive way would involve multiplying two polynomials of degree p (one of which has 0 for half it's coefficients) and recursively solving the coefficients (which isn't as bad as it sounds). What's really nice is once you've found these coefficients you can vary m as much as you like and finding the new value of your power sum just involves evaluating a p+1 degree polynomial.

If you look at your proposed method, you have to find all these $$R_{pmk}$$ one for each value of k from 1 to m. Each of these involves p multiplications. On the surface this is the same work it would take to find this sum in the obvious way, multiply each number from 1 to m by itself p times then add (though in practice finding a pth power is considerably faster than multiplying p different numbers together). If you then wanted to vary m, it looks like you'd have to go back and recalculate all these R's again. This isn't even taking into account the Q's you have to find using their sum (up to p of these are non-trivial).