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, S10(1000) = 110 + 210 + … + 100010 = 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, A0,j is computed as 1 / (j+1). In subsequent rows, Ai,j is computed as (j+1) × (Ai−1,j − Ai−1,j+1). The Bernoulli number Bn is located in the first column at An,0.
Given Bernoulli numbers, sums of powers can be computed by Bernoulli’s formula, where 
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 Sm(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.
Programming Praxis 
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.
#!/usr/bin/env python import numpy as np from fractions import Fraction from operator import mul def a_t(n): A = np.zeros(shape=(n + 1, n + 1), dtype=Fraction) for j in xrange(n + 1): A[0, j] = Fraction(1, j + 1) for i in xrange(1, n + 1): for j in xrange(n - i + 1): A[i, j] = (j + 1) * (A[i - 1, j] - A[i - 1, j + 1]) return A def b(n): return a_t(n)[n][0] def c(n, k): num = reduce(mul, xrange(n - k + 1, n + 1), 1) denom = reduce(mul, xrange(1, k + 1), 1) return num / denom def s_naive(m, n): return sum(pow(k, m) for k in xrange(1, n + 1)) def s(m, n): s = sum(c(m + 1, k) * b(k) * pow(n, m + 1 - k) for k in xrange(m + 1)) return s // (m + 1) if __name__ == "__main__": print b(18) print a_t(18)[:, 0] print s_naive(10, 3) print s(10, 3) print s_naive(10, 1000) print s(10, 1000)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):
import Data.Ratio a :: (Integral a, Integral b) => a -> a -> Ratio b a i j = iterate (\xs -> zipWith (*) [1..] $ zipWith (-) xs (tail xs)) (map (1 %) [1..]) !! fromIntegral i !! fromIntegral j bernoullis :: (Integral a, Integral b) => a -> [Ratio b] bernoullis upto = map (flip a 0) [0..upto] choose :: Integral a => a -> a -> Ratio a choose n k = product [1..n] % (product [1..n-k] * product [1..k]) s :: Integral a => a -> a -> Ratio a s m n = 1 % (m+1) * sum [choose (m+1) k * a k 0 * (n%1)^(m+1-k) | k <- [0..m]][…] today’s Programming Praxis exercise, our goal is to implement an algorithm that calculates the bernoulli […]