## Sums Of Powers

### February 11, 2011

We saw the calculation of the Bernoulli numbers in the previous exercise. In today’s exercise we will study the Bernoulli numbers in their original context in computing the sums of powers. For instance, *S*_{10}(1000) = 1^{10} + 2^{10} + … + 1000^{10} = 91409924241424243424241924242500, a result calculated by the Swiss mathematician Jakob Bernoulli in “less than half of a quarter of an hour” in the late seventeenth century.

Any discussion of Bernoulli numbers is complicated by the fact that there are several similar sequences that go by that name. The oldest is -1/2, 1/6, 0, -1/30, 0, 1/42, 0, -1/30, 0, 5/66, …; that’s the sequence we calculated previously. The modern version of the sequence prepends a 1, so it’s 1, -1/2, 1/6, 0, -1/30, 0, 1/42, 0, -1/30, 0, 5/66, …. And sometimes the -1/2 is given as +1/2, so the sequence is 1, 1/2, 1/6, 0, -1/30, 0, 1/42, 0, -1/30, 0, 5/66, …. Even worse, sometimes the Bernoulli numbers are numbered counting from zero, though mostly they are numbered counting from one, so you have to be careful when you talk about the Bernoulli numbers.

The particular sequence that we want to look at is the third sequence that starts with 1, 1/2, 1/6, 0, -1/30, …. In 1999, S. Akiyama and Y. Tanigawa found a simple way to compute that sequence in linear time using a dynamic program represented in a triangular matrix:

1 | 1/2 | 1/3 | 1/4 | 1/5 |

1/2 | 1/3 | 1/4 | 1/5 | … |

1/6 | 1/6 | 3/20 | … | … |

0 | 1/30 | … | … | … |

-1/30 | … | … | … | … |

In the first row, A_{0,j} is computed as 1 / (*j*+1). In subsequent rows, A_{i,j} is computed as (*j*+1) × (A_{i−1,j} − A_{i−1,j+1}). The Bernoulli number B_{n} is located in the first column at A_{n,0}.

Given Bernoulli numbers, sums of powers can be computed by Bernoulli’s formula, where is the binomial coefficient:

Your task is to write functions that return a list of Bernoulli numbers up to a requested limit and that compute the sum of powers S_{m}(*n*). When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

I used this as an opportunity to play with NumPy’s (http://numpy.scipy.org/)

n-dimensional arrays. Also uses the fractions module for rational numbers.

Small typo that’s bothering me. It’s equivalent, but line 17 should read

`return a_t(n)[n, 0]`

instead to be consistent with the rest.

I believe the formula given for s in the exercise is wrong: it should be sum of k = 1 to m + 1 (or 0 to m), not 1 to n. At least that’s the version I found on wikipedia that worked for me after the 1 to n one didn’t.

My Haskell solution (see http://bonsaicode.wordpress.com/2011/02/11/programming-praxis-sums-of-powers/ for a version with comments):

[…] today’s Programming Praxis exercise, our goal is to implement an algorithm that calculates the bernoulli […]