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 \left( \begin{array}{c} n \\ k \end{array} \right) = \frac{n!}{(n-k)!k!} is the binomial coefficient:

S_m(n) = \frac{1}{m+1} \sum_{k=1}^n \left( \begin{array}{c} m+1 \\ k \end{array} \right) B_k n^{m+1-k}

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.

About these ads

Pages: 1 2

4 Responses to “Sums Of Powers”

  1. Graham said

    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)
    
  2. Graham said

    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.

  3. 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]]
    
  4. [...] today’s Programming Praxis exercise, our goal is to implement an algorithm that calculates the bernoulli [...]

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 627 other followers

%d bloggers like this: